Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Posting-Version: version B 2.10.2 9/3/84; site panda.UUCP
Path: utzoo!decvax!genrad!panda!sources-request
From: sources-request@panda.UUCP
Newsgroups: mod.sources
Subject: IEEE Calculator (part 4 of 6)
Message-ID: <863@panda.UUCP>
Date: Tue, 3-Sep-85 22:11:44 EDT
Article-I.D.: panda.863
Posted: Tue Sep  3 22:11:44 1985
Date-Received: Wed, 4-Sep-85 23:10:30 EDT
Sender: jpn@panda.UUCP
Lines: 1391
Approved: jpn@panda.UUCP

Mod.sources:  Volume 3, Issue 6
Submitted by: decvax!decwrl!sun!dgh!dgh (David Hough)

#! /bin/sh
: make a directory, cd to it, and run this through sh
echo If this kit is complete, "End of Kit" will echo at the end
echo Extracting base.i
cat >base.i <<'End-Of-File'
(* File base.i, Version 8 October 1984.  *)

procedure mpyten ( var x : internal ; n : integer  ) ;

        (* Multiplies x by 10**n, using table of powers of ten.  *)
        
var
n1, n2 : integer ;

begin
n1 := abs(n) div 32 ;
n2 := abs(n) mod 32 ;
if n1 < 32 then begin
if n > 0 then begin
if n2 > 0 then multiply( x, tensmall[n2], x ) ;
if n1 > 0 then multiply( x, tenbig[n1], x) ;
end
else if n < 0 then begin
if n2 > 0 then divide ( x, tensmall[n2], x ) ;
if n1 > 0 then divide ( x, tenbig[n1], x ) ;
end ;
end
else begin (* n is too big.  *)
if n > 0 then begin
makeinf(x) ;
setex ( overfl ) ;
end
else begin
makezero(x) ;
setex ( underfl ) ;
end ;
end ;
end ;


procedure xtodec ( x : internal ; r : roundtype ;
                var s : strng ) ;
                
                (* Converts abs(x) to an integral value, then that
                value is converted to a strng of ASCII digits.  *)
                

var
tf : excepset ;
acc : -20..+20  ; (* divide accumulator *)
i, j : integer ;
carry : boolean ;
nib : nibarray ;
last : integer ;

begin
roundint ( x, r, xprec ) ;
fpstatus.curexcep := fpstatus.curexcep - [inxact] ; 
        (* Don't care about inxact.  *)
if kind(x) = zerokind then makeucsdstring('0 ',s) else begin
s[0] := chr(0) ;
last := x.exponent - 1 ; (* Last is the last active bit of the accumulator.  *)

repeat (* Get one digit per cycle until there's nothing left.  *)

(* For each digit, do a NON RESTORING divide by TEN.  *)

acc := - 10 ;

for j := 0 to last do begin 
        (* Do one divide minicycle for each bit of dividend, plus one extra.  *)
acc := acc + acc ; (* Double remainder.  *)
if x.significand[j] then acc := acc + 1 ; (* Shift in next bit of dividend.  *)

if acc < 0 then begin (* Do add-ten cycle.  *)
acc := acc + 10 ;
end
else begin (* Do subtract-ten cycle.    *)
acc := acc - 10 ;
end ;
x.significand[j] := acc >= 0  ; (* Sign of this step determines quotient bit
                        and whether add or subtract next time.  *)
(* End around complement rotate for quotient bit.  *)
end (* Divide mini-cycle.  *) ;

if acc < 0  then begin 
        (* Remainder is negative so add 10.  *)
acc := acc + 10 ;
end ;

precatchar(' ',s) ;
s[1] := chr(ord('0') + acc ) ;

j := firstbit ( x, 0, last ) ;
if j <= last then 
        begin
        left(x, j) ;
        end ;
last := last - j ;
until last < 0 ;
        (* Keep doing divide cycles until the quotient is zero.  *)

end ;

end ;

procedure subdec (* x : internal ; p1, p2 : integer ; var s: strng *) ;
        (* s receives a strng of decimal digits representing the integer in
        x.significand[p1]..x.significand[p2], right justified.  *)
        
var
j, i : integer ;
nib : nibarray ;

begin

if p2 = stickybit then begin (* Avoid trying to donormalize sticky bit.  *)
for i := p1 to p2 do
x.significand[i-1] := x.significand[i] ;
p1 := p1 - 1 ;
p2 := p2 - 1 ;
end ;

for i := 0 to (p1-1) do 
x.significand[i] := false ; (* Clear bits outside field.  *)
for i := (p2+1) to stickybit do
x.significand[i] := false ;

x.exponent := p2 + 1 ;
donormalize(x) ;
xtodec( x, rnear, s) ;
end ;


procedure findbinten ( x : internal ; n : integer ;
                var s : strng ; var p : integer  ) ;
                
        (* Converts x into s and p, s a strng of exactly n significant
        digits, so
                x .=. s * 10^p *)

var
e1, e2, dp : integer ;
tf : excepset ;
t : internal ;
cc : conditioncode ;
i : integer ;
norm : boolean ;
et, dt : integer ;
fraction : integer ;
spill : boolean ;

begin
        (* First ESTIMATE p.
                We want a p such that
                10^(n-1) <= int(abs(x)*10^p) < 10^n
                
                For a first guess, use
                n - log10( 2**e * 1+f ) 
                which we approximate by
                n - ((77+(1/16))/256) * (e + f )
                
                
                e is broken into two pieces to get benefit of 
                22 bit product.  *)
                
norm := x.significand[0] ; (* If x is not normalized then we don't force
        n significant digits in E format.  *)
e2 := (x.exponent-1) div 256 ; (* e = e1 + 256*e2 *)
e1 := (x.exponent-1) - 256 * e2 ;
fraction := xbyte(x,1,8) ;
e1 := 77 * e1 + 16 * e2 ; (* First order contribution from e1 and second
        order from e2. *)
if norm then e1 := e1 + ( 77 * fraction ) div 256 ;
        (* If normalized add in a contribution for the fraction.
        If not normalized, assume significand of 1.00000...  *)
et := e1 div 256 ; dt := e1 - et * 256 ;
        (* We are never sure how Pascal div and mod work on negative
        numbers.  *)
if dt > 0  then et := et + 1 ;
        (* but we try to get et as close as possible anyway to
        the ceiling of e1/256.  *)
e2 := 77 * e2 ;
p := n  - (et + e2) ;

        (* Now remedy flaws in approximation by altering p as necessary.  *)

tf := fpstatus.curexcep  ; (* Save exceptions.  *)
dp := 0 ; (* Assume no correction required.  *)
repeat
fpstatus.curexcep  := tf ; 
(* Restore exception status.  Don't want inexact flag set
                for inappropriate p.  *)
p := p + dp ; (* Correct p.  *)
dp := 0 ; (* Assume no correction required.  *)
t := x ;
mpyten(t, p) ; (* Multiply t by appropriate power.  *)
roundint( t, fpstatus.mode.round, xprec ) ;

t.sign := false ; (* Use absolute values for comparison.  *)
compare ( t, tensmall[n], cc ) ;  (* t must not exceed n sig digits.  *)
case cc of
otherwise ;
greater : dp := -1 ; (* If too big, correct and repeat.  *)
equal : begin
        if norm or (n=1) then begin (* We want only n digits.  *)
        p := p - 1 ; (* If almost, avoid full repeat of process.  *)
        t := tensmall[n-1] ;
        end 
        else begin (* If not normalized we want n-1 digits.  *)
        p := p - 2 ;
        t := tensmall[n-2] ;
        end ;
        end ;
lesser : begin
        compare(t,tensmall[n-1],cc) ;
         if norm then begin (* If normalized, insure enough sig digits. *)
        if cc = lesser then dp :=  + 1 ; (* If not enough digits, correct *)
        end  (* Need exactly n digits.  *)
        else begin (* If unnormalized, want no more than n-1 digits.  *)
        if cc <> lesser then dp := -1 ; (* Try again for less than n.  *)
        end ;
        end ;
end ;
t.sign := x.sign ;
{
if dp <> 0 then writeln(' Power of ten correction: ',dp) ;
}
spill := ([underfl, overfl] * fpstatus.curexcep ) <> [] ;
until (dp = 0) or (kind(t)=zerokind) or spill ;
        (* Repeat until no correction necessary or over/underfl occurs.  *)
        (* or the number rounds to normalized zero.  *)
if spill then
        begin
        (* String of asterisks if overfl.  *)
        s[0] := chr(0) ;
        while length(s) < n do concatchar(s,'*') ;
        end 
else 
        begin
        xtodec(*2*) ( t, fpstatus.mode.round, s(*, n*) ) ; (* Get strng.  *)
        end ;
while length(s) < n do precatchar( '0', s ) ;
        (* Add unnormalizing digits.  *)
p := - p ;
end ;

procedure decint ( s : strng ; var x : internal ; var e : integer ) ;

        (* DECINT converts a strng s of decimal digits into x and e
        such that
                s = x * 10^e
        If s >= 2^(leastsigbit+1) then the sticky bit may be set. *)

const
last = 69 ; (* last bit of decint accumulator *)
var
i,j : integer ;
acc : array[0..last] of boolean ; (* Accumulator long enough to hold 20
        significant digits and a sticky bit *)
carry, zero : boolean ;
n : nibarray ;

procedure tenmult ; (* multiplies acc by 10 *)
var
i : integer ;
carry : boolean ;

begin
for i := 0 to (last-3) do acc[i] := acc[i+3] ; (* multiply acc by 8 *)
for i := (last-2) to last do acc[i] := false ;
carry := false ;
for i := (last-1) downto 2 do
adder( acc[i], acc[i-2], acc[i], carry) ;
for i := 1 downto 0 do
adder( acc[i], false, acc[i], carry) ;
end ;

begin (* decint *)
for i := 0 to last do acc[i] := false ;
e := 0 ;
zero := true ;
for j := 1 to length(s) do begin
if s[j] = '0' then begin
if not zero then e := e + 1 ;
end 
else begin
if not zero then begin
e := e + 1 ;
while (e>0) and (not acc[0]) and (not acc[1]) and (not acc[2]) and 
(not acc[3]) do begin (* multiply by ten *)
tenmult ;
e := e - 1 ;
end ;
end ;
zero := false ;
if e > 0 then acc[last] := true (* set sticky bit on *)
else begin (* add digit *)
hexnibble( s[j], n) ;
carry := false ;
for i := (last-1) downto (last-4) do
adder( acc[i], n[i+4-last], acc[i], carry) ;
for i := (last-5) downto 0 do
adder( acc[i], false, acc[i], carry) ;
end ;
end ;
end ;

if zero then makezero(x) else begin
        (* Now acc[0..(last-1)] represents an integer value,
        acc[last] a sticky bit, to be multiplied by 10^e *)
        
        i := 0 ;
        while ( i < last ) and not acc[i] do i := i + 1 ;
                (* search for first nonzero bit *)
        x.exponent := last - i ; (* Set exponent of result *)
        for j := 0 to (last-i-1) do acc[j] := acc[j+i] ;
                (* normalize *)
        for j := (last-i) to (last-1) do acc[j] := false  ; 
        for i := 0 to (stickybit-1) do x.significand[i] := acc[i] ;
        x.significand[stickybit] := acc[last] ;
        for i := stickybit to (last-1) do
        x.significand[stickybit] := x.significand[stickybit] or acc[i] ;
        end ;
end ;

procedure putdec (* s : strng ; p1, p2 : integer ; 
                var x : internal ; var error : boolean *) ;
                
                (* Interprets s as a decimal integer, puts value in bits
                p1..p2 of x.significand.
                Sets Error if any significant bits don't fit in field.  *)

var
i, j : integer ;
y : internal ;
e : integer ;
f : excepset ;

begin
decint( s, y, e ) ;
f := fpstatus.curexcep ;
mpyten ( y, e ) ;
fpstatus.curexcep := f ; (* Ignore exceptions that may arise.  *)
e := y.exponent ;
error := (e > (leastsigbit+1)) ;
        (* Bad news if s too big to fit in 64 bits.  *)
if not error then begin
if kind(y) = zerokind then begin
for i := p1 to p2 do
x.significand[i] := false ; (* Set up zero field.  *)
end
else begin
if (p2-p1+1) >= e then begin (* y fits in field.  *)
for i := p2 downto (p2+1-e) do
x.significand[i] := y.significand[i+e-1-p2] ;
for i := (p2-e) downto p1 do
x.significand[i] := false ;
end
else begin
for i := p2 downto p1 do
x.significand[i] := y.significand[i+e-1-p2] ;
for i := (p1-p2+e-2) downto 0 do
error := error or y.significand[i] ; (* Check for lost significant bits.  *)
end end end
end ;

procedure todecint ( x : internal ; var s : strng ) ;

(* if x is an integer less than 2**15,
then s receives the decimal digits representing x.
Otherwise s is set to empty. *)

var
i, k : integer ;

begin
s[0] := chr(0) ;
if kind(x) = zerokind then makeucsdstring('0 ',s) else
if (abs(kind(x)) = normkind) and (x.exponent <= 15) and (x.exponent >= 1)
then begin
if zerofield ( x, x.exponent, stickybit ) then begin (* it's all integer *)
k := 0 ;
for i := 0 to (x.exponent-1) do begin (* Accumulate k.  *)
k := k + k ;
if x.significand[i] then k := k + 1 ;
end ;
while k > 0 do begin
precatchar( chr(ord('0') + (k mod 10)), s) ;
k := k div 10 ;
end ;
if x.sign then precatchar( '-', s )  ;
end end
end ;

procedure bindec (* x : internal ; var s  : strng *)  ;
(* converts x to decimal format *)

var
e, i, j, k : integer ;
nib : nibarray ;
t : strng ;
tf : excepset ;
ns : integer ;

begin
case abs(kind(x)) of
zerokind : if x.sign 
	then makeucsdstring('-0',s) else  makeucsdstring('0 ',s) ;

unnormkind, normkind : begin
todecint ( x, s ) ;
if length(s) < 1 then begin (* Can't represent as integer; too bad.  *)
ns := 19 ;
case storagemode of (* Set number of significant digits output.  *)
flt32 : ns := 9 ;
f64 : ns := 17 ;
otherwise ;
end ;
findbinten( x, ns, s, e ) ;
tf := fpstatus.curexcep ;
fpstatus.curexcep := [] ;
if not (overfl in tf) then begin
if e <> 0 then begin (* x not an integer so write it in E format.  *)
e := e + ns-1 ;
for i := length(s) downto 2 do s[i+1] := s[i] ;
s[2] := '.' ;
s[0] := chr(length(s)+1) ;
if e <> 0 then begin
concatchar(s, 'E') ;
if e > 0 then concatchar(s, '+') ;
intdec(e, t) ;
s := concat(s,t) ;
end ;
end ;
end ;
if x.sign then precatchar('-',s) ;
end ;
end ;

infkind, nankind : nanascii ( x, false, s ) ;

otherwise
end ;
end ;



procedure decbin (* s : strng ; var x : internal ; var error : boolean *) ;
(* converts decimal strng s to internal format *)
(* error is set true if bad format *)

type
stringclass = (nonnumeric, truezero, nonzero) ; (* types of strng *)

var
class : stringclass ;
i, k,  min : integer ;
e1, e2 : integer ;
sub : strng ;
t : strng ;
esign : boolean ;
nib : nibarray ;
ee, ie : integer ;


procedure checkadd ( x, y : integer ; var z : integer ; var error : boolean ) ;
        (* Computes z := x + y except if z overflows error is set to true
                and z is set to maxexp-1 or minexp+1 *)
begin
error := false ;
z := x + y ;
if (x>0) and (y>0) and (z<=0) then begin
z := maxexp - 1 ;
error := true ;
end
else if (x<0) and (y<0) and (z>=0) then begin
z := minexp + 1 ;
error := true ;
end ;
end ;

procedure bump ; (* removes first character from strng t *)
begin
delete (t,1,1) 
end ;


begin
class := nonnumeric ;
error := false ;
esign := false ;
x.sign := false ;
x.exponent := 0 ;
ee := 0 ; ie := 0 ;
for i := 0 to stickybit do x.significand[i] := false ;
sub[0] := chr(0) ; (* substring for accumulating significant digits *)
t[0] := chr(0) ;
for i := 1 to length(s) do if s[i] <> ' ' then concatchar(t,upcase(s[i])) ;
concatchar(t,'!') ; (* this marks the end of the input strng *)

if t[1] = '+' then bump else if t[1] = '-' then begin (* handle negative *)
x.sign := true ;
bump
end ;
while t[1] = '0' do begin
class := truezero ;
bump ; (* delete leading zeros *)
end ;
while t[1] in digitset do begin (* digits before point *)
class := nonzero ;
concatchar(sub, t[1]) ;
bump
end ;
if t[1] = '.' then begin (* check for point *)
bump ;
while t[1] in digitset do begin (* process digits after point *)
if (t[1] <> '0') or (class = nonzero) then class := nonzero 
else class := truezero ;
concatchar( sub, t[1]) ;
ie := ie - 1 ;
bump ; 
end ;  
end ;
ee := 0 ;
if t[1] = 'E' then bump ; (* handle E for exponent *)
if t[1] = '+' then bump else if t[1]='-' then begin (* exponent sign *)
esign := true ;
bump
end ;
while t[1] in digitset do begin (* exponent digits *)
if ee > ((maxexp - (ord(t[1])-ord('0'))) div 10 ) then begin
error := true ;
ee := maxexp - 1 ;
end else
begin
ee := 10 * ee + ord(t[1]) - ord('0') ;
end ;  bump  end ;
if class = truezero then x.exponent := minexp  else begin
if esign then ee := -ee ;
checkadd(ee,ie,ee,error) ; (* ee := ee + ie *)
if not error then begin 
                (* Minimize ee if possible by adding zeros to sub *)
ie := 19 - length(sub) ; (* Maximum number of zeros to add.  *)
if (ee>0) and (ie>0) then begin (* Go ahead and add.  *)
if ee < ie then ie := ee ; (* Only add enough to reduce ee to 0.  *)
ee := ee - ie ;
for i := 1 to ie do concatchar( sub, '0') ;
end ;
decint ( sub, x, ie ) ; (* Convert substring to x and ie *)
checkadd( ee, ie, ee, error ) ; (* Add in ie to exponent.  *)
end ;
if not error then 
mpyten ( x, ee ) ; (* Adjust x by appropriate power of ten.  *)
end ;
if class = nonnumeric  then 
        (* the following code checks for INFs and NANs *)
begin
NANer ( s, false, x, error ) 
end 
else
if length(t) > 1 then error := true ;
if error  then 
        begin
        makenan(nanascbin,x) ;
        end ;
end ;

procedure display (* x : internal *) ;
(* displays x in decimal and binary *)

begin
write(' Hex: ') ; displayhex(x) ;
write(' Dec: ') ; displaydec(x) ;
end ;
End-Of-File
echo Extracting utility.i
cat >utility.i <<'End-Of-File'
(* File utility.i, Version 8 October 1984. *)

function length ( var x : strng ) : integer ;

begin (* concat *)
length := ord(x[0]) ;
end   (* concat *) ;

procedure displayhex ( x : internal ) ;

var s : strng ;
i : integer ;

begin
binhex(x,s) ;
for i := 1 to length(s) do write(s[i]) ;
writeln ;
end ;

procedure displaydec ( x : internal ) ;

var s : strng ;
i : integer ;

begin
bindec(x,s) ;
for i := 1 to length(s) do write(s[i]) ;
writeln ;
end ;

procedure concatchar ( var s : strng ; c : char ) ;
(* concatenates the character c onto the end of s *)
var
ls : integer ;
begin
ls := ord(s[0]) + 1 ;
s[ls] := c ;
s[0] := chr(ls) ;
end ;

function upcase ( c : char ) : char ;
begin
if ('a' <= c) and (c <= 'z') then upcase := chr(ord(c)-32) else upcase := c
end  ;

function stackspace : integer ;

        (* Computes number of stack entries left.
        As this number approaches zero, operation becomes dangerous.  *)

var space : integer ;

begin
stackspace := 10000 ;
end ;

procedure hexnibble ( h : char ; var n : nibarray ) ;
(* Converts ASCII hexit h into a nibarray *)
var
i, w : integer ;
begin
if h in digitset then w := ord(h)-ord('0') else w := ord(h) - ord('A') + 10 ;
for i := 3 downto 0 do begin
n[i] := odd(w) ;
w := w div 2 ;
end ;
end ;

function nibblehex (* n : nibarray ) : char *)  ;
(* converts a nibarray into a hexit ASCII character *)
var
i, w : integer ;
c : char ;

begin
w := 0 ;
for i := 0 to 3 do begin
w := w + w ;
if n[i] then w := w + 1 ;
end ;
if w < 10 then c := chr(ord('0') + w) else c := chr(ord('A') + w - 10 ) ;
nibblehex := c ;
end ;

procedure displayexcep ( es : excepset ) ;
        (* Displays exception names that are enabled.  *)
var i : xcpn ;
begin
for i := invop to inexact do
if i in es then write(' ',xcpnname[i],' ') ;
end ;

procedure displaystatus ;
        (* Displays current mode, trap, exception flags.  *)

begin
write(' Modes: ') ;
case fpstatus.mode.round of
rneg : write(' RM ') ;
rpos : write(' RP ') ;
rzero : write(' RZ ') ;
otherwise
end ;
case fpstatus.mode.precision of
sprec: write(' R24 ') ;
dprec: write(' R53 ') ;
otherwise
end ;
if fpstatus.mode.clos = proj then write(' PROJ ') ;
if fpstatus.mode.norm = warning then write(' WARN ') ;
case storagemode of
i16 : write(' I16 ') ;
i32 : write(' I32 ') ;
i64 : write(' I64 ') ;
flt32 : write(' F32 ') ;
f64 : write(' F64 ') ;
ext80 : write(' X80 ') ;
otherwise
end ;
writeln ;
if fpstatus.trap <> [] then begin (* Write out trap line.  *)
write(' Traps: ') ;
displayexcep( fpstatus.trap ) ;
writeln ;
end ;
if fpstatus.excep <> [] then begin (* Write out exceptions.  *)
write(' Exceptions: ') ;
displayexcep( fpstatus.excep ) ;
writeln ;
end ;
end ;

procedure trapmessage ;

        (* Announces name of trap that would occur now.  *)
        
var
tset : excepset ;
f : xcpn ;

begin
tset := fpstatus.trap * fpstatus.curexcep ;
if tset <> [] then begin
f := invop ; (* Start with highest priority exception.  *)
while (f <> cvtovfl) and not (f in tset) do f := succ(f) ;
if f in tset then writeln( xcpnname[f],' TRAP occurs.  ') ;
end ;
end ;

procedure setex (* e : xcpn *) ;
        (* Turns on exception flag in curexcep.  *)
begin
fpstatus.curexcep := fpstatus.curexcep  + [e] ;
end ;

function zerofield (* x : internal ; p1, p2 : integer ) : boolean *) ;

        (* Returns true if x.significand[p1..p2] is all zeros.  *)
        
var i : integer ;

begin
i := p1 ;
while ( i < p2 ) and not x.significand[i] do i := i + 1 ;
zerofield := ( i >= p2 ) and not x.significand[p2]  ;
        (* Can't test bit p2 in main loop ; would cause range error if
        p2 were stickybit, on subsequent test.  *)
end ;

function firstbit (* x : internal ; p1, p2 : integer ) : integer *)  ;

        (* Returns index of leftmost onebit in field
        x.significand[p1..p2].
        If field is zero, returns p2+1.  *)
        
var i : integer ;

begin
i := p1 ;
while ( i < p2 ) and not x.significand[i] do i := i + 1 ;
if ( i >= p2 ) and not x.significand[p2] then i := p2+1   ;
        (* Can't test bit p2 in main loop ; would cause range error if
        p2 were stickybit, on subsequent test.  *)
firstbit := i ;
end ;

function lastbit (* x : internal ; p1, p2 : integer ) : integer *)  ;

        (* Returns index of rightmost nonzero bit in field
        x.significand[p1..p2].
        If field is zero, returns p1-1.  *)
        
var i : integer ;

begin
i := p2 ;
while ( i > p1 ) and not x.significand[i] do i := i - 1 ;
if ( i <= p1 ) and not x.significand[p1] then i := p1 - 1  ;
        (* Can't test bit p1 in main loop ; would cause range error if
        p1 were zero, on subsequent test.  *)
lastbit := i ;
end ;

function kind (* x : internal ) : integer *)  ;
(* returns kind(x) but all NANs have kind=4 in order to fit int16 *)
var
i, tkind : integer ;
begin
if x.exponent = maxexp then begin (* inf or nan *)
if zerofield ( x, 1, stickybit )  then tkind := ord(infkind) 
else tkind := ord(nankind) ;
end
else 
if (x.exponent <= minexp) and zerofield(x, 0, stickybit)
then tkind := ord(zerokind)
else if x.significand[0] = true then tkind := ord(normkind)
else tkind := ord(unnormkind) ;
if x.sign then tkind := -tkind ;
kind := tkind ;
end ;

procedure makezero (* var x : internal *) ;

        (* makes x into a zero.  Does not change sign of x *)

var i : integer ;

begin
x.exponent := minexp ;
for i := 0 to stickybit do x.significand[i] := false ;
end ;


procedure makeinf (* var x : internal *) ;

        (* makes x into a infinity.  Does not change sign of x *)

var i : integer ;

begin
x.exponent := maxexp ;
for i := 0 to stickybit do x.significand[i] := false ;
end ;

procedure makenan (* n : integer ; var x : internal  *) ;

        (* makes x a NAN and inserts n in its more significant field.
        Sets NV Operand flag in curexcep. *)

var
i : integer ;

begin
x.exponent := maxexp ;
for i := 0 to stickybit do x.significand[i] := false ;
i := 15 ;
while n <> 0 do begin
x.significand[i] := odd(n) ;
n := n div 2 ;
i := i - 1 ;
end ;
setex( invop ) ;
end ;

function unzero (* x : internal ) : boolean *)  ;

        (* returns TRUE if x is an unnormalized zero,
        FALSE otherwise *)


begin
unzero := zerofield( x, 0, stickybit ) and (x.exponent > minexp)  ;
end ;

procedure pushstack (* x : internal *) ;
(* pushes x on stack *)
        (* In case of NV exception and trapping NAN, makes the NAN
        non-trapping.  *)

var
p : pstack  ;

begin
if stackspace >= 3 then begin
new(p) ;
if (invop in fpstatus.curexcep) and (abs(kind(x))=nankind) then
x.significand[1] := false ; (* By convention bit 1 determines trapping/
                                non-trapping.  *)
p^.x := x ;
p^.next := stack ;
stack := p ;
end else 
writeln(' ERROR: not enough space for push! ') ;
end ;

procedure popstack (* var x : internal *) ;
        (* pops stack to x, or sets x to 0 if stack is empty *)
begin
if stack=nil then begin
x.sign := false ;
makezero(x) ;
end
else begin
x := stack^.x ;
stack := stack^.next ;
if (abs(kind(x))=nankind) and x.significand[1] then (* It's a Trapping NAN. *)
setex( invop ) ;
end
end ;

procedure donormalize (* var x : internal *) ;
        (* normalizes x *)
        (* Unnormalized zeros are set to normalized zeros.
a           INFs and NANs are not changed *)

var
i, j : integer ;

begin
if x.exponent < maxexp then begin
i := firstbit( x, 0, stickybit )  ;
if i > stickybit then x.exponent := minexp (* zero *) else
if i > 0 then begin
x.exponent := x.exponent - i ;
for j := i to stickybit do 
x.significand[j-i] := x.significand[j] ;
for j := ((stickybit+1)-i) to (stickybit-1) do
x.significand[j] := x.significand[stickybit] ;
{
if (x.significand[stickybit]) and (i>1)  then 
writeln(i,' ERROR: Normalizing Sticky Bit ') ;
        (* It's OK to shift sticky bit into next to last position because
        during rounding the last two positions are stuck together.
        It is definitely not OK to shift a sticky bit any further left.  *)
}
end ;
end ;
end ;

procedure right (* var x : internal ; n : integer *) ;
        (* does sticky right shift of internal x *)
        
var
i : integer ;

begin
{
if (0 > n) or (n > stickybit) then
writeln(' Funny Right ',n) ;
}
if n > stickybit then n := stickybit ; (* It's all the same for large n.  *)
x.significand[stickybit] := not zerofield( x, (stickybit-n), stickybit) ;
for i := (stickybit-1) downto n do
x.significand[i] := x.significand[i-n] ;
for i := (n-1) downto 0 do
x.significand[i] := false ;
end ;

procedure left (* var x : internal ; n : integer *) ;

        (* Lefts shifts significand of x, n times *)
        
var i : integer ;

begin
{
if (0 > n) or ( n > stickybit) then
writeln(' Funny left ',n) ;
}
if n > stickybit then n := stickybit ;  (* All the same for large n.  *)
for i := 0 to (stickybit-1-n) do
x.significand[i] := x.significand[i+n] ;
{
if x.significand[stickybit] and (n>1)  then
writeln(n,' Error: LEFT shift of STICKY bit ') ;
}
for i := (stickybit-n) to (stickybit-1) do
x.significand[i] := x.significand[stickybit] ;
end ;

procedure roundkcs (* var x : internal ; r : roundtype ;
        p : precisetype  *) ;
        
        (* Rounds x according to rounding mode and rounding precision.
        Sets Inexact flag in curexcep if appropriate. *)
        
var i : integer ;
akx : integer ;

procedure dorn ; (* round to nearest *)

var
i : integer ;
carry : boolean ;

begin
carry := true ;
i := (leastsigbit+1) ;
while (i>=0) and carry do begin
adder(x.significand[i], false, x.significand[i], carry) ;
i := i-1 ;
end ;
if carry then begin (* carry out of most significant bit occurred. *)
x.significand[0] := true ;
x.exponent := x.exponent + 1 ;
end
else begin (* Check for ambiguous case *)
if zerofield( x, leastsigbit+1, stickybit )  
then (* It is the ambiguous case. *)
x.significand[leastsigbit] := false ; (* force round to even. *)
end ;
end ;

procedure doro ; (* Round away from zero (Outward Round) *)

var
i : integer ;
carry : boolean ;

begin
carry := false ;
for i := (leastsigbit+1) to stickybit do carry := carry or x.significand[i] ;
if carry then begin (* propagate a carry *)
i := leastsigbit ;
while (i >= 0) and carry do begin
adder(x.significand[i], false, x.significand[i], carry) ;
i := i - 1 ;
end ;
if carry then begin (* Carry out occurred, so renormalize. *)
x.significand[0] := true ;
x.exponent := x.exponent + 1 ;
end ;
end ;
end ;

begin (* round *)
akx := abs(kind(x)) ;
if akx in [unnormkind, normkind] then begin
case p of 
sprec: begin
right(x, 40) ;
x.exponent := x.exponent + 40 ;
end ;

dprec : begin
right(x, 11) ;
x.exponent := x.exponent + 11 ;
end ;

otherwise
end ;
for i := (leastsigbit+1) to stickybit do if x.significand[i] then 
begin
setex( inxact ) ;
end ;
case r of
rnear : begin 
dorn ;
end ;
rneg : if x.sign then doro  ;
rpos : if not x.sign then doro ;
otherwise
end ;

for i := (leastsigbit+1) to stickybit do x.significand[i] := false ;
        (* Eliminate G, R, and S bits. *)
case p of
 sprec: begin
left(x,39) ;
x.exponent := x.exponent - 39 ; 
if not  x.significand[0] then 
begin
left( x, 1) ;
x.exponent := x.exponent - 1 ;
end ;
end ;

dprec: begin
left(x,10) ;
x.exponent := x.exponent - 10 ; 
if not  x.significand[0] then 
begin
left(x,1) ;
x.exponent := x.exponent - 1 ;
end ;
end ;

otherwise
end ;
end ;
end ;

procedure roundint (* var x : internal ; r : roundtype ; p : precisetype *) ;
        
        (* Rounds x to an integral value in accordance with modes.  *)

var
akx, i, count : integer ;

begin
akx := abs(kind(x)) ;
if akx in [unnormkind, normkind]  then begin
if (x.exponent >= (leastsigbit+1)) then count := 0
else if x.exponent <= (leastsigbit+1-stickybit) then count := stickybit
else count := (leastsigbit+1) - x.exponent ; 
        (* Compute shift count of bits to get rid of.  *)
case p of (* But allow for rounding to shorter precisions, too.  *)
 sprec: if count < 40 then count := 40 ;
 dprec: if count < 11 then count := 11 ;
otherwise
end ;
if count > 0 then right ( x, count) ;
roundkcs( x, r, xprec) ; (* Do rounding.  *)
if count > leastsigbit then begin (* Limit left shifts
        for 0 < x < 1 which must be rounded either to 0 or 1.  *)
        count := leastsigbit ;
        x.exponent := 1 ;
        end ;
if count > 0 then begin
left(x, count-1 ) ;
if x.significand[0] then x.exponent := x.exponent + 1 (* Rounding carry out.  *)
else left(x, 1) ;
end ;
if zerofield ( x, 0, stickybit ) then x.exponent := minexp ;
        (* No significant bits left so make it a true zero.  *)
end end ;


procedure picknan (* x, y : internal ; var z : internal *) ;

        (* Sets z to whichever of x or y is a NAN.
        If both are NANs, sets z to the one with the 
        greatest significand.  *)
        
var i : integer ;

begin
if abs(kind(x)) = nankind then 
if abs(kind(y)) = nankind then begin
i := 0 ;
while (i <= leastsigbit) and (x.significand[i] = y.significand[i]) do 
        i := i + 1 ;
if x.significand[i] then z := x else z := y ;
end 
else z := x else z := y
end ;

function equalinternal ( x1, x2 : internal ) : boolean ;
        (* Returns true if x1 = x2 *)
var
t : boolean ;
i : integer ;

begin
t := (x1.sign=x2.sign) and (x1.exponent=x2.exponent) ;
if t then for i := 0 to stickybit do
t := t and (x1.significand[i]=x2.significand[i]) ;
equalinternal := t ;
end ;

function concat ( x, y : strng ) : strng ;

var t : strng ;
i, lx, ly : integer ;

begin (* concat *)
t := x ;
lx := ord(x[0]) ; ly := ord(y[0]) ;
for i := 1 to ly do t[lx+i] := y[i] ;
t[0] := chr(lx+ly) ;
concat := t ;
end   (* concat *) ;

procedure precatchar ( c : char ; var x : strng ) ;

var i, ls : integer ;

begin (* precatchar *)
ls := ord(x[0]) ;
for i := ls downto 1 do x[i+1] := x[i] ;
x[1] := c ;
x[0] := chr(ls+1) ;
end   (* precatchar *) ;

procedure delete ( var x : strng ; index, count : integer ) ;

var i : integer ;

begin (* delete *)
for i := index+count to length(x) do
	x[i-count] := x[i] ;
x[0] := chr(length(x)-count) ;
end   (* delete *) ;

procedure makeucsdstring( x : strng; var t : strng );
	(* Converts a constant string to UCSD form. *)
var i, l : integer ;
begin
l := 0 ;
while (33 <= ord(x[l])) and (ord(x[l]) <= 126) do l := l + 1 ;
for i := l downto 1 do t[i] := x[i-1] ;
t[0] := chr(l) ;
end ;

procedure copy ( s : strng; index, count : integer ; var t : strng ) ;

var i : integer ; l : integer ;

begin (* copy *)
t[0] := chr(count) ;
for i := 1 to count do t[i] := s[index+i-1] ;
end   (* copy *) ;

procedure insert( s : strng ; var d : strng ; index : integer ) ;

var i,ld,ls : integer ;

begin (* insert *)
ls := ord(s[0]) ; ld := ord(d[0]) ;
for i := ld downto index do d[i+ls] := d[i] ;
for i := ls downto 1 do d[index+i-1] := s[i] ;
d[0] := chr(ls+ld) ;
end   (* insert *) ;

function pos ( c : char ; s : strng ) : integer ;

var i, l : integer ;

begin (* pos *)
l := ord(s[0]) ;
i := 1 ;
while (i <= l) and (s[i] <> c) do i := i + 1 ;
if i <= l then pos := i else pos := 0 ;
end   (* pos *) ;

function sequal ( s, c : strng ) : boolean ;
      (* Compares UCSD string s to C string c and returns true if equal. *)
var
i, ls, lc : integer ;

begin (* sequal *)
lc := 0;
ls := ord(s[0]) ;
while (33 <= ord(c[lc])) and (ord(c[lc]) <= 126) do lc := lc + 1 ;
if lc <> ls then 
	begin
	sequal := false ;
	end
 else
	begin
	i := 1 ;
	while (i <= ls) and (s[i] = c[i-1]) do i := i + 1 ;
	sequal := i > ls ;
	end ;
end   (* sequal *) ;


End-Of-File
echo Extracting init.i
cat >init.i <<'End-Of-File'

(* File init.i, Version 4 February 1985. *)


PROCEDURE initialize ;

        (* does all the initializing of tables and variables.  *)

var error : boolean ;

procedure inittensmall ;

        (* procedure to initialize the table of small  powers of ten *)
        
var
i : integer ;
j : integer ;
x : internal ;
carry : boolean ;
last : integer ;
error : boolean ;

begin
(* Make 10^0=1 *)
for i := 1 to stickybit do x.significand[i] := false ;
x.sign := false ;
x.exponent := 1 ;
x.significand[0] := true ;
tensmall[0] := x ;

(* Make other exact powers of ten.  *)
last := 0 ; (* Last non-zero bit. *)
for j := 1 to 28 do begin
x.exponent := x.exponent + 3 ; (* Multiply by 8 first.  *)
last := last + 2 ; (* At least 2 more significant bits.  *)
carry := false ;
for i := last downto 2 do
adder(x.significand[i-2], x.significand[i],x.significand[i],carry) ;
for i := 1 downto 0 do
adder(false, x.significand[i],x.significand[i],carry) ;
if carry then begin (* Overflowed slightly.  *)
x.exponent := x.exponent + 1 ;
last := last + 1 ;
for i := last downto 1 do
x.significand[i] := x.significand[i-1] ;
x.significand[0] := true ;
end ;
tensmall[j] := x ;
end ;

hexbin(' .a18f 07d7 36b9 0be5 4 h  97', tensmall[29], error ) ;
hexbin(' .c9f2 c9cd 0467 4ede c h 100', tensmall[30], error ) ;
hexbin(' .fc6f 7c40 4581 2296 4 h 103', tensmall[31], error ) ;
end ;

procedure inittenbig ;

        (* procedure to initalize the table of large powers of ten *)
        
var
error : boolean ;

begin
tenbig[0] := tensmall[0] ;
hexbin(' .9dc5 ada8 2b70 b59d c h 107',tenbig[1], error ) ;
hexbin(' .c278 1f49 ffcf a6d5 4 h 213',tenbig[2], error ) ;
hexbin(' .efb3 ab16 c59b 14a2 c h 319',tenbig[3], error ) ;
hexbin(' .93ba 47c9 80e9 8cdf c h 426',tenbig[4], error ) ;
hexbin(' .b616 a12b 7fe6 17aa 4 h 532',tenbig[5], error ) ;
hexbin(' .e070 f78d 3927 556a c h 638',tenbig[6], error ) ;
hexbin(' .8a52 96ff e33c c92f c h 745',tenbig[7], error ) ;
hexbin(' .aa7e ebfb 9df9 de8d c h 851',tenbig[8], error ) ;
hexbin(' .d226 fc19 5c6a 2f8c 4 h 957',tenbig[9], error ) ;
hexbin(' .8184 2f29 f2cc e375 c h 1064',tenbig[10], error ) ;
hexbin(' .9fa4 2700 db90 0ad2 4 h 1170',tenbig[11], error ) ;
hexbin(' .c4c5 e310 aef8 aa17 4 h 1276',tenbig[12], error ) ;
hexbin(' .f28a 9c07 e9b0 9c58 c h 1382',tenbig[13], error ) ;
hexbin(' .957a 4ae1 ebf7 f3d3 c h 1489',tenbig[14],error) ;
hexbin(' .b83e d8dc 0795 a262 4 h 1595',tenbig[15], error) ;
end ;
 
procedure inittenhuge ;

var
error : boolean ;

begin
hexbin(' .e319 a0ae a60e 91c6 c h 1701',tenbig[16], error) ;
hexbin(' .8bf6 1451 432d 7bc2 c h 1808', tenbig[17], error) ;
hexbin(' .ac83 fb89 6b67 95fc c h 1914', tenbig[18], error) ;
hexbin(' .d4a4 4fb4 b8fa 79af c h 2020', tenbig[19], error) ;
hexbin(' .830c f791 e54a 9d1c c h 2127', tenbig[20], error) ;
hexbin(' .a188 4b69 ade2 4964 4 h 2233', tenbig[21], error) ;
hexbin(' .c71a a36a 1f8f 01cb c h 2339', tenbig[22], error) ;
hexbin(' .f56a 298f 4370 28f3 4 h 2445', tenbig[23], error) ;
hexbin(' .973f 9ca8 cd00 a68c 4 h 2552', tenbig[24], error) ;
hexbin(' .ba6d 9b40 d7cc 9ecc c h 2658', tenbig[25], error) ;
hexbin(' .e5ca 5a0b 8d73 7f0e 4 h 2764', tenbig[26], error) ;
hexbin(' .8d9e 89d1 1346 bda5 4 h 2871', tenbig[27], error) ;
hexbin(' .ae8f 2b2c e3d5 dbe9 c h 2977', tenbig[28], error) ;
hexbin(' .d729 3020 5a0c 1b2f c h 3083', tenbig[29], error) ;
hexbin(' .849a 672a 0d2e cfd1 c h 3190', tenbig[30], error) ;
hexbin(' .a372 2c13 41fa 93de 4 h 3296', tenbig[31], error) ;
end ;



begin
digitset := [ '0' .. '9' ] ;
hexset := digitset + [ 'A' .. 'F' ] ;
stack := nil ;
storagemode := unrounded ;
testflag := false ;
fpstatus.mode.round := rnear ;
fpstatus.mode.precision := xprec ;
fpstatus.mode.clos := affine ;
fpstatus.mode.norm := normalizing ;
fpstatus.curexcep := [] ;
fpstatus.excep := [] ;
fpstatus.trap := [] ;
leftnan[1] := 0 ;
leftnan[2] := 24 ;
leftnan[3] := 53 ;
leftnan[4] := leastsigbit + 1 ;
rightnan[1] := leftnan[2] - 1 ;
rightnan[2] := leftnan[3] - 1 ;
rightnan[3] := leftnan[4] - 1 ;
 rightnan[4] := stickybit ;
xcpnname[cvtovfl] := 'IV' ;
xcpnname[overfl] := 'OV' ;
xcpnname[underfl] := 'UN' ;
xcpnname[div0] := 'D0' ;
xcpnname[invop] := 'NO' ;
xcpnname[inxact] := 'NX' ;
inittensmall ;
inittenbig ;
inittenhuge ;
(*
decbin ( ' 3.1415926535 89793 23846 26433', pi, error ) ;
decbin ( ' 2.7182818284 59045 23536 02874', e, error ) ;
*)
hexbin ( ' .c90fdaa22168c234c h 2', pi, error ) ;
decbin ( ' 2.7182818284 59045 23536 02874', e, error ) ;
end ;


End-Of-File
echo ""
echo "End of Kit"
exit