Usenet 1994 January/usenetsourcesnewsgroupsinfomagicjanuary1994.iso/sources/unix/volume3/ieee/part6

Internet Message Format | 1986-11-30 | 26.2 KB

From: genrad!decvax!decwrl!sun!dgh!dgh (David Hough)
Subject: IEEE Calculator (part 6 of 6)
Newsgroups: mod.sources
Approved: [e-mail]

Mod.sources:  Volume 3, Issue 8
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 forward.i
cat >forward.i <<'End-Of-File'
(* File forward.i, Version 8 October 1984.  *)

                        (* FORWARDs for CALC   *)

procedure adder ( x, y : boolean ; var z : boolean ; var carry : boolean ) ;
forward ;
procedure suber ( x, y : boolean ; var z : boolean ; var carry : boolean ) ;
forward ;

procedure bindec ( x : internal ; var s : strng ) ;
forward ;

procedure binhex ( x : internal ; var s : strng ) ;
forward ;

procedure compare ( x, y : internal ; var cc : conditioncode ) ;
forward ;

procedure add ( x, y : internal ; var z : internal ) ; forward ;

procedure decbin ( s : strng ; var x : internal ; var error : boolean ) ;
forward ;

procedure hexbin ( s : strng ; var x : internal ; var error : boolean ) ;
forward ;

procedure putdec ( s : strng ; p1, p2 : integer ; var x : internal ;
                        var error : boolean ) ;
forward ;

procedure subdec ( x : internal ; p1, p2 : integer ; var s : strng ) ;
forward ;

function nibblehex ( n : nibarray ) : char ;
forward ;

procedure setex ( e : xcpn ) ; forward ;

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

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

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

function kind ( x : internal ) : integer ;
forward ;

procedure makezero ( var x : internal ) ;
forward ;

procedure makeinf ( var x : internal ) ;
forward ;

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

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

procedure donormalize ( var x : internal ) ;
forward ;

procedure right ( var x: internal ; n : integer ) ;
forward ;

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

procedure roundkcs ( var x: internal ; r : roundtype ; p : extprec ) ;
forward ;

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

procedure roundint ( var x : internal ; r : roundtype ; p : extprec ) ;
forward ;

procedure unpackinteger ( y : cint64 ; var x : internal ; itype : inttype ) ;
forward ;

procedure store ( var x : internal ) ;
forward ;

procedure display ( x : internal ) ; forward ;

procedure popstack ( var x : internal ) ; forward ;

procedure pushstack ( x : internal ) ; forward ;

procedure dotest ( s : strng ; var found : boolean ; x, y : internal ) ;
forward ;


End-Of-File
echo Extracting function.i
cat >function.i <<'End-Of-File'
(* File Calc:Function, Version 24 May 1984.    *)

procedure compare (* x, y : internal ; var cc : conditioncode *) ;

        (* computes X comp Y and returns result cc *)


procedure docomp ;
        (* does bitwise X comp Y and returns result in cc *)
        (* Works like -INF < -NORM < -0 < +0 < +NORM < +INF
        so don't use to compare two zeros or to compare with INF
        in proj mode *)

var i : integer ;

begin
cc := equal ;
if x.sign <> y.sign then
if x.sign then cc := lesser else cc := greater
else begin (* same signs *)
if x.exponent < y.exponent then cc := lesser
else if x.exponent > y.exponent then cc := greater
else begin (* same sign and same exponent *)
i := 0 ;
while ( i < stickybit) and (x.significand[i]=y.significand[i]) do i := i + 1 ;
if i <> stickybit then
if x.significand[i] then cc := greater else cc := lesser ;
end ;
if x.sign then case cc of (* if negative numbers, reverse conclusion *)
lesser : cc := greater ;
greater : cc := lesser ;
end ;
end ;
end ;

begin (* compare *)
roundkcs ( x, fpstatus.mode.round, xprec ) ; (* Preround; ignore inxact. *)
donormalize(x) ;
roundkcs (y, fpstatus.mode.round, xprec ) ; (* Preround; ignore inxact. *)
donormalize(y) ;
fpstatus.curexcep := fpstatus.curexcep - [inxact] ; (* Ignore.  *)
case abs(kind(x)) of

zerokind: case abs(kind(y)) of
zerokind : cc := equal ;
normkind : docomp ;
infkind : if fpstatus.mode.clos = proj then cc := notord
else docomp ;
nankind : cc := notord ;
end ;

normkind : case abs(kind(y)) of
zerokind , normkind : docomp ;
infkind : if fpstatus.mode.clos = proj then cc := notord
else docomp ;
nankind : cc := notord ;
end ;

infkind : case abs(kind(y)) of
zerokind, normkind : if fpstatus.mode.clos = proj then
        cc := notord else docomp ;
infkind : if fpstatus.mode.clos = proj then cc := equal
else docomp ;
nankind : cc := notord ;
end ;

nankind : cc := notord ;
end ;

end ;

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

        (* Does z := x + y  *)

procedure doadd ;

        (* Does true add z := x + y
        Neither x nor y should be INF or NAN
        x and y should not both be true zero.  *)

var
carry : boolean ;
i : integer ;
shiftcount : integer ;

begin
roundkcs( x, fpstatus.mode.round, xprec) ; (* Pre-round.  *)
roundkcs( y, fpstatus.mode.round, xprec) ;

z.sign := x.sign ;

if x.exponent >= y.exponent then begin (* Align smaller operand.  *)
z.exponent := x.exponent ;
shiftcount := x.exponent - y.exponent ;
if shiftcount < 0 then shiftcount := maxexp ;
right( y, shiftcount ) ;
end
else begin
z.exponent := y.exponent ;
shiftcount := y.exponent - x.exponent ;
if shiftcount < 0 then shiftcount := maxexp ;
right( x, shiftcount ) ;
end ;

carry := false ;
for i := stickybit downto 0 do
adder( x.significand[i], y.significand[i], z.significand[i], carry ) ;

if carry then begin (* Renormalize for carry out.  *)
right( z, 1 ) ;
z.significand[0] := true ;
z.exponent := z.exponent + 1 ;
end ;

end ;

procedure dosub ;

        (* Does true subtract z := x - y,
        neither of which should be INF or NAN,
        only one of which should be true zero.  *)

var
carry : boolean ;
shiftcount, i : integer ;
postnorm : boolean ;
rnd : roundtype ;
redo : boolean ;
xur, yur : internal ; (* Save x and y unrounded operands.  *)

begin
        (* Proper preround is ambiguous in the case when the indicated
        rounding mode is RZ and a true subtract is required.
        x and y should be rounded RM if the result is positive,
        RP if the result is negative.
        Occasionally the result comes out reversed and the operands
        have to be re-pre-rounded.
        All this makes a good argument for not having pre-rounding or at
        least fudging this one annoying case.  *)

y.sign := not y.sign ; (* Reverse sign of y so x and y have same signs.  *)
xur := x ; yur := y ; (* Save unrounded operands.  *)
redo := false ; (* Set true if we go through loop twice.  *)

repeat
x := xur ; y := yur ; (* Restore unrounded state.  *)
rnd := fpstatus.mode.round ; (* This is almost always appropriate.  *)

if x.exponent >= y.exponent then begin
z.sign := x.sign ;
z.exponent := x.exponent ;
if rnd = rzero then begin
if x.sign <> redo then rnd := rpos else rnd := rneg end  ;
roundkcs ( x, rnd, xprec ) ;
roundkcs ( y, rnd, xprec ) ;
shiftcount := x.exponent - y.exponent ;
if shiftcount < 0 then shiftcount := maxexp ;
right ( y, shiftcount ) ;
end
else begin
z.sign := not y.sign ;
z.exponent := y.exponent ;
if rnd = rzero then begin  (*  Bad case.  *)
if y.sign <> redo then rnd := rpos else rnd := rneg end  ;
roundkcs( x, rnd, xprec ) ;
roundkcs( y, rnd, xprec ) ;
shiftcount := y.exponent - x.exponent ;
if shiftcount < 0 then shiftcount := maxexp ;
right ( x, shiftcount ) ;
end ;

postnorm := x.significand[0] or y.significand[0] ;
        (* Post normalization occurs if larger operand is normalized.  *)

carry := false ;
if z.sign = x.sign then for i := stickybit downto 0 do
suber( x.significand[i], y.significand[i], z.significand[i], carry )
else for i := stickybit downto 0 do
suber( y.significand[i], x.significand[i], z.significand[i], carry ) ;

if carry then begin (* Sign reversal.  *)
z.sign := not z.sign ;
carry := true ; (* For end-around carry.  *)
for i := stickybit downto 0 do
adder( not z.significand[i], false, z.significand[i], carry ) ;
(* Complement.  *)
redo := fpstatus.mode.round = rzero ; (* Pre-round was inappropriate.  *)
if redo then writeln(' RE-PRE-ROUND required for SUBTRACT.  ') ;
end ;

until not redo ;

if postnorm then donormalize(z) ; (* Do renormalization.  *)

store(z) ; (* Force round to storage mode to determine if the result will
                be zero; if so, special sign rules apply.  *)
if zerofield( z, 0, leastsigbit )  then (* Correct sign for zero result.  *)
z.sign := fpstatus.mode.round = rneg ;  (* Which is +0 except in RM mode.  *)

end ;

begin (* add *)
if (abs(kind(x))=nankind) or (abs(kind(y))=nankind) then
picknan( x, y, z) else

case abs(kind(x)) of

zerokind : case abs(kind(y)) of

zerokind : begin (* +-0 + +-0 *)
z := x ;
if x.sign <> y.sign then z.sign := fpstatus.mode.round = rneg ;
        (* +0 + -0 usually has sign +0 except in RM mode.  *)
end ;
unnormkind, normkind : if x.sign = y.sign then doadd else dosub ;
infkind : z := y ;
end ;

unnormkind, normkind : if abs(kind(y)) = infkind then z := y else
if x.sign = y.sign then doadd else dosub ;

infkind : if abs(kind(y)) <> infkind then z := x else (* INF +- INF *)
if (fpstatus.mode.clos = proj) or (x.sign <> y.sign)
then makenan ( nanadd, z )
else z := x ;

end ;

end ;

procedure multiply (x , y : internal ; var z : internal  ) ;

        (* routine to multiply two internal numbers and return z := x*y
        with curexcep  containing flags set. *)

var dorout : boolean ;

procedure domult ;

        (* does the multiply of z := abs(x) * abs(y) *)

var
i, j, j0, n0, n1  : integer ;
carry : boolean ;
r : roundtype ;
xlast, ylast : integer ;
xfirst, yfirst : integer ;

                (* The following subprocedures do various cases of
                multiply substeps.  Based on Booth algorithm ideas.  *)

procedure don0 ;

        (* Multiplies z by n0 zeros, i. e. right shifts n0 times,
        and resets n0.  *)

var i : integer ;

begin
z.significand[stickybit] := not zerofield ( z, stickybit-n0, stickybit ) ;
                (* Shift bits into stickybit.  *)
for i := (stickybit-1) downto n0 do
z.significand[i] := z.significand[i-n0] ; (* Really shift bits.  *)
for i := (n0-1) downto 0 do
z.significand[i] := false ; (* Shift in zeros at left.  *)
n0 := 0 ;
end ;

procedure do11 ;

        (* Does add y and shift to z.  *)

var
j : integer ;
carry : boolean ;

begin
if z.significand[stickybit-1] then z.significand[stickybit] := true ;
for j := (stickybit-1) downto (ylast+2) do  (* These bits only shift.  *)
z.significand[j] := z.significand[j-1] ;
carry := false ;
for j := (ylast+1) downto (yfirst+1) do
adder( z.significand[j-1], y.significand[j-1], z.significand[j], carry ) ;
z.significand[yfirst] := carry ;
end ;

procedure do21 ;

        (* Does twice:  add y and shift z.  *)

begin
do11 ; do11 ;
end ;

procedure don1 ;

        (* Does n1 times: add y and shift z, by
        1) subtract y
        2) shift z n1 times
        3) add 2*y
                                        *)

var
j, j0 : integer ;
carry : boolean ;

begin
if n1=1 then do11 else if n1=2 then do21 else begin
                        (* Do subtract.  *)
carry := false ;
for j := (ylast) downto (yfirst) do
suber( z.significand[j], y.significand[j], z.significand[j], carry ) ;
for j := (yfirst-1) downto 0 do (* Ripple carry.  *)
z.significand[j] := true ;
                        (* Do shift.  *)
z.significand[stickybit] := not zerofield( z, stickybit-n1, stickybit ) ;
for j := (stickybit-1) downto n1 do
z.significand[j] := z.significand[j-n1] ;
for j := (n1-1) downto 0 do
z.significand[j] := true  ;
                        (* Do add 2*y.  *)
carry := false ;
for j := ylast downto yfirst do
adder( z.significand[j], y.significand[j], z.significand[j], carry ) ;
end ;
for j := (yfirst-1) downto 0 do
z.significand[j] := false ;
n1 := 0 ;
end ;

begin
        (* Preround. *)
dorout := false ;
case fpstatus.mode.round of
rnear, rzero : r := fpstatus.mode.round ;
rneg : if x.sign = y.sign then r := rzero else dorout := true ;
rpos : if x.sign = y.sign then dorout := true   else r := rzero ;
end ;
if dorout then
        begin (* round out *)
        if x.sign then roundkcs( x, rneg, xprec ) else roundkcs( x, rpos, xprec ) ;
        if y.sign then roundkcs( y, rneg, xprec ) else roundkcs( y, rpos, xprec ) ;
        end   (* round out *)
        else
        begin
        roundkcs( x, r, xprec ) ;
        roundkcs( y, r, xprec ) ;
        end ;

xfirst := firstbit( x, 0, leastsigbit) ;
yfirst := firstbit( y, 0, leastsigbit ) ;
if xfirst <= leastsigbit then
xlast := lastbit( x, xfirst, leastsigbit) else xlast := -1 ;
if yfirst <= leastsigbit then
ylast := lastbit( y, yfirst, leastsigbit) else ylast := -1 ;

if (xfirst > xlast) or (yfirst > ylast)  then begin
(* z is unnormalized zero.  *)
makezero(z) ;
z.exponent := x.exponent + y.exponent - 1 ;
end

else begin (* Both operands are non-zero.  *)

z.exponent := x.exponent + y.exponent ;
for i := 0 to stickybit do z.significand[i] := false ;
n1 := 0 ; n0 := 0 ;
for i := xlast downto xfirst do  (* for each bit of x *)
if x.significand[i] then
begin
        (* Encounter One.  *)
if n0 > 0 then begin
don0 ;
n1 := 1 ;
end
else n1 := n1 + 1
end

else begin
        (* Encounter Zero.  *)
if n1 > 0 then begin
don1 ;
n0 := 1 ;
end
else n0 := n0 + 1
end ;

if n1 > 0 then don1 ; (* Tidy up at end.  *)
n0 := n0 + xfirst ;
(* Additional right shifts necessary if x not normalized.  *)
if n0 > 0 then don0 ;

if not z.significand[0] then begin (* Do one donormalize cycle *)
z.exponent := z.exponent - 1 ;
left(z, 1 ) ;
end ;
end ;
if (x.exponent < 0) and (y.exponent < 0) and (z.exponent > 0) then begin
        (* Gross underfl has caused integer overfl leading to large
        positive exponent.  *)
right(z,stickybit) ;
z.exponent := minexp ;
setex(underfl) ;
end ;
end ;

begin
if (abs(kind(x))=nankind) or (abs(kind(y))=nankind) then picknan(x,y,z) else
case abs(kind(x)) of
zerokind : case abs(kind(y)) of
zerokind, unnormkind,  normkind : z := x ;
infkind : makenan(nanmul, z) ;
end ;

unnormkind : case abs(kind(y)) of
zerokind: z := y ;
unnormkind, normkind : domult ;
infkind : if unzero(y) then makenan( nanmul, z) else z := y ;
end ;

normkind : case abs(kind(y)) of
zerokind : z := y ;
unnormkind, normkind : domult ;
infkind : z := y ;
end ;

infkind : case abs(kind(y)) of
zerokind : makenan( nanmul, z ) ;
unnormkind : if unzero(y) then makenan(nanmul, z) else z := x ;
normkind , infkind : z := x ;
end ;
end ;

z.sign := (x.sign <> y.sign) ; (* exclusive OR signs *)

end ;

procedure divide ( x, y : internal ; var z : internal  ) ;

        (* Does extended internal divide z := x/y *)


procedure divbyzero ; (* for invalid divide exception *)
begin
makezero(z) ;
z.exponent := maxexp ; (* Make Infinity result *)
setex( div0 ) ;
end ;

procedure dodiv ;

        (* carries out divide of x by normalized y *)
        (* x should be nonzero and finite *)
        (* signs are ignored *)

var
i,j : integer ;
carry, sbit, tbit : boolean ;
r : internal ;
rx, ry : roundtype ;
rlast, ylast  : integer ;
xrout, yrout : boolean ;

begin
        (* Preround. *)
xrout := false ; yrout := false ;
case fpstatus.mode.round of
rnear : begin
rx := rnear ;
ry := rnear ;
end ;
rzero : begin
rx := rzero ;
yrout := true ;
end ;
rneg : if x.sign = y.sign then begin
rx := rzero ;
yrout := true ;
end
else begin
xrout := true ;
ry := rzero ;
end ;
rpos: if x.sign = y.sign then begin
xrout := true ;
ry := rzero ;
end
else begin
rx := rzero ;
yrout := true ;
end ;
end ;
if xrout then
        begin
        if x.sign then roundkcs( x, rneg, xprec ) else roundkcs ( x, rpos, xprec )
        end
        else roundkcs( x, rx, xprec) ;
if yrout then begin
        if y.sign then roundkcs( y, rneg, xprec ) else roundkcs ( y, rpos, xprec )
        end
else roundkcs( y, ry, xprec) ;

ylast := lastbit ( y, 0, leastsigbit ) ;
rlast := lastbit ( x, 0, leastsigbit ) + 1 ;
if rlast < (ylast+1) then rlast := ylast + 1  ;

for i := 0 to (rlast-1) do
r.significand[i+1] := x.significand[i] ; (* remainder R := X/2 *)
r.significand[0] := false ;
z.exponent := x.exponent - y.exponent + 1 ;
sbit := false ; (* Sbit contains the sign of the remainder between steps *)

for i := 0 to (stickybit-1) do begin (* for each bit of result *)
tbit := sbit ; (* T bit holds test to decide + or -.  *)
sbit := r.significand[ylast+1] ;
for j := (ylast+1) to (rlast-1) do
r.significand[j] := r.significand[j+1] ;
carry := false ;

if tbit then begin (* If R < 0 then R := 2 * (R+Y) *)
for j := ylast downto 0 do begin
adder(sbit, y.significand[j], tbit, carry) ;
sbit := r.significand[j] ;
r.significand[j] := tbit ;
end ;
adder( sbit, false, sbit, carry) ; (* Sbit now has the sign of R *)
end
else begin (* If R >= 0 then R := 2 * (R-Y) *)
for j := ylast  downto 0 do begin
suber(sbit, y.significand[j], tbit, carry) ;
sbit := r.significand[j] ;
r.significand[j] := tbit ;
end ;
suber(sbit, false, sbit, carry) ; (* Sbit now has the sign of R *)
end ;
r.significand[rlast] := false ; (* result of left shift *)
z.significand[i] := not sbit ; (* Result bit for this step *)
end ;

(* Next check for sticky bit.  Result Z is exact iff
        R = 0 or R <= -2Y  *)

z.significand[stickybit] := true ; (* Result is almost always inexact *)
if sbit then begin (* R < 0 so compute R + 2Y *)
carry := false ;
for j := rlast downto 0 do
adder(r.significand[j], y.significand[j], r.significand[j], carry) ;
z.significand[stickybit] :=  carry ; (* If no carry then R+2Y < 0 *)
end ;
if z.significand[stickybit] then begin (* check if R >=0 or R+2Y >= 0 *)
(* R >= 0 ; Is R = 0 ? *)
z.significand[stickybit] := not zerofield ( r, 0, rlast )  ;
end ;


if not z.significand[0] then begin (* normalize once *)
z.exponent := z.exponent - 1 ;
for i := 1 to stickybit do z.significand[i-1] := z.significand[i] ;
end ;

if (x.exponent < 0) and (y.exponent < 0) and (z.exponent > 0) then begin
        (* Gross underfl has caused integer overfl leading to large
        positive exponent.  *)
right(z,stickybit) ;
z.exponent := minexp ;
setex(underfl) ;
end ;
end ;



begin (* divide *)
case abs(kind(x)) of

zerokind: case abs(kind(y)) of
zerokind: makenan( nandiv, z) ;
unnormkind: if unzero(y) then makenan(nandiv, z) else z := x ;
normkind, infkind : z := x ;
nankind : z := y ;
end ;

unnormkind : case abs(kind(y)) of
zerokind : if unzero(x) then makenan(nandiv, z) else divbyzero ;
unnormkind : makenan( nandiv, z) ;
normkind :  dodiv ;
infkind : makezero(z) ;
nankind : z := y ;
end ;

normkind : case abs(kind(y)) of
zerokind : divbyzero ;
unnormkind : makenan(nandiv,z) ;
normkind : dodiv ;
infkind : makezero(z) ;
nankind : z := y ;
end ;

infkind : case abs(kind(y)) of
zerokind, unnormkind, normkind : z := x ;
infkind : makenan(nandiv,z) ;
nankind : z := y ;
end ;

nankind : picknan( x, y, z) ;
end ;

z.sign := x.sign <> y.sign ; (* Do Exclusive Or *)
end ;


End-Of-File
echo Extracting global.i
cat >global.i <<'End-Of-File'
(* File Calc:Global, Version 24 May 1984.     *)

(* global constant, type, and variable declarations for  calc *)

const
leastsigbit = 63 ;
(* Position of least significant bit in external extended.  *)

maxexp  = 32767 ; (* 2**15-1, maximum internal exponent *)
(* used for INF and NAN *)
minexp  = -maxexp ;
(* -2**15+1, minimum internal exponent, used for zero *)

biasex = 16383 ; (* Extended exponent bias.  *)
maxex = 16384 ; (* Extended maximum unbiased exponent.  *)
minex = -16383 ; (* Extended minimum unbiased exponent.  *)

biased = 1022 ; (* Double exponent bias.  *)
maxed = 1025 ; (* Double maximum unbiased exponent.  *)
mined = -1022 ; (* Double minimum unbiased exponent.  *)

biases = 126 ; (* Single exponent bias.  *)
maxes = 129 ; (* Single maximum unbiased exponent.  *)
mines = -126 ; (* Single minimum unbiased exponent.  *)

zerokind = 0 ; (* KIND of zero *)
unnormkind = 1 ; (* KIND of denormalized or unnormalized *)
normkind = 2 ; (* KIND of normalized number *)
infkind = 3 ; (* KIND of infinity *)
nankind = 4 ; (* KIND of NAN *)
negunnorm = -1 ;
negnorm = -2 ;
neginf = -3 ;
negnan = -4 ;

ordbell = 7 ; (* ASCII code for Bell.  *)

nanfields = 4 ; (* Number of fields in NAN significand.  *)
{
notord = unord ;
}
nansqrt = 1 ;
nanadd = 2 ;
nanint = 3 ;
nandiv = 4 ;
nantrap = 5 ;
nanunord = 6 ;
nanproj = 7 ;
nanmul = 8 ;
nanrem = 9 ;
nanresult = 12 ;
nanascbin = 17 ;
nanascnan = 18 ;
naninteger = 20 ;
nanzero = 21 ;

type

strng = fpstring ; (* Standard string type.  *)

inttype = i16 .. i64 ; (* Types of integer operands.  *)

sxinternal = record (* Special signless single extended internal format.  *)
        exponent : integer ;
        significand : array[0..1] of integer ;
        end ;

pstack = ^ stacktype ;

stacktype = record (* item on stack *)
x : internal ;
next : pstack ;
end ;

nibarray = array[0..3] of boolean ;

var
fpstatus : fpstype ; (* current status *)
storagemode : arithtype ; (* current storage mode *)
        (* Each operand is rounded to this mode after operation.  *)
testflag : boolean ; (* True if DOTEST is to be called.  *)
stack : pstack  ;
digitset, hexset : set of char ;

tensmall : array [ 0 .. 31 ] of internal ; (* Table of small powers of ten *)
tenbig : array [ 0 .. 31 ] of internal ; (* Table of big powers of ten *)
pi, e : internal ; (* Familiar constants.  *)

leftnan, rightnan : array [ 1 .. nanfields ] of 0 .. stickybit ;
        (* Indexes of bitfield boundaries for NAN significands.  *)

xcpnname : array [ exception                               ] of
        packed array [1..2] of char ;
        (* Two character codes for exceptions.  *)


End-Of-File
echo Extracting addsubpas.i
cat >addsubpas.i <<'End-Of-File'
procedure adder (* x, y : boolean ; var z : boolean ; var carry : boolean *);

        (* computes boolean z := x + y with carry in and out *)


begin
z := carry ;
carry := z and x ;
z := ( z <> x ) ;
if carry then z := y else begin
carry := (z and y) ;
z := (z <> y) ;
end ;
end ;

procedure suber (* x, y : boolean ; var z : boolean ; var carry : boolean *) ;


        (* computes boolean z := x - y with carry in and out *)

begin
z := y <> carry ;
carry := carry and y ;
if carry then z := x else begin
carry := (z and (not x)) ;
z := (z <> x) ;
end ;
end ;
End-Of-File
echo Extracting divrem.i
cat >divrem.i <<'End-Of-File'
(* File Calc:Divrem, Version 19 February 1982.  *)

procedure divrem ( x, y : internal ; var q, r : internal  ) ;

        (* Computes q := x DIV y and r := x REM y.  *)

procedure dodivrem ;

        (* Computes DIV and REM for normalized y and normalized or
        unnormalized x.  *)

var
i,j : integer ;
carry, sbit, tbit, zbit : boolean ;
rlast, ylast  : integer ;
nc : integer ;
roundup : boolean ;

begin
        (* Preround. *)
roundkcs( x, fpstatus.mode.round, xprec) ;
roundkcs( y, fpstatus.mode.round, xprec) ;
makezero(q) ; (* Starting assumption.  *)
r := x ; (* Starting assumption.  *)
nc := x.exponent - y.exponent + 1 ; (* Number of cycles to produce result. *)
if nc >= 0 then begin
ylast := lastbit ( y, 0, leastsigbit ) ;
rlast := lastbit ( x, 0, leastsigbit ) + 1 ;
if rlast < (ylast+1) then rlast := ylast + 1  ;

for i := 0 to (rlast-1) do
r.significand[i+1] := x.significand[i] ; (* remainder R := X/2 *)
r.significand[0] := false ;
sbit := false ; (* Sbit contains the sign of the remainder between steps *)

for i := 1 to nc do begin (* for each bit of result *)
tbit := sbit ; (* T bit holds test to decide + or -.  *)
sbit := r.significand[ylast+1] ;
for j := (ylast+1) to (rlast-1) do
r.significand[j] := r.significand[j+1] ;
carry := false ;

if tbit then begin (* If R < 0 then R := 2R+Y *)
for j := ylast downto 0 do begin
adder(sbit, y.significand[j], tbit, carry) ;
sbit := r.significand[j] ;
r.significand[j] := tbit ;
end ;
adder( sbit, false, sbit, carry) ; (* Sbit now has the sign of R *)
end
else begin (* If R >= 0 then R := 2R-Y *)
for j := ylast  downto 0 do begin
suber(sbit, y.significand[j], tbit, carry) ;
sbit := r.significand[j] ;
r.significand[j] := tbit ;
end ;
suber(sbit, false, sbit, carry) ; (* Sbit now has the sign of R *)
end ;
r.significand[rlast] := false ; (* result of left shift *)
if (leastsigbit-nc+i) >= 0 then
q.significand[leastsigbit-nc+i] := not sbit ;
end ;

r.exponent := r.exponent - nc + 1  ;
for i := 0 to (leastsigbit-nc) do
q.significand[i] := false ;  (* Fill in bits.  *)
carry := false ;
if sbit then (* R < 0 so set R := R + Y.  *)
for j := rlast downto 0 do
adder(r.significand[j], y.significand[j], r.significand[j], carry ) ;
if  zerofield( r, 0, rlast) then
        (* R=0 so q is exact and r is zero.  *)
        makezero(r)
        else begin (* At this point
        2R < Y implies q is ok
        2R = Y implies round q to even and set r to +- .5 Y
        Y < 2R < 2Y implies round q upward, set r to r-y.  *)
roundup := false ; (* Roundup controls rounding of q, and thus r.  *)
carry := false ;
zbit := false ;
for j := rlast downto 1 do begin (* Compute sign of 2R - Y .*)
suber(r.significand[j], y.significand[j-1], tbit, carry ) ;
zbit := zbit or tbit ;
end ;
suber(r.significand[0], false, tbit, carry ) ;
zbit := zbit or tbit ; (* zbit is false if 2R = Y.  *)
if not carry then begin (* 2R >= Y.  *)
roundup := zbit (* 2R>Y *) or q.significand[leastsigbit] (* 2R=Y *) ;
end ;
if roundup then begin (* Increment q; reverse r.  *)
carry := true ;
for j := leastsigbit downto 0 do
adder(q.significand[j], false, q.significand[j], carry) ;
carry := false ;
for j := (leastsigbit+1) downto 0 do
suber( y.significand[j], r.significand[j], r.significand[j], carry ) ;
        (* R := Y - R.  *)
r.sign := not r.sign ;
end end ;
donormalize(r) ;
q.exponent := 64 ;
donormalize(q) ;
end ;
end ;

begin (* divrem *)

if (abs(kind(x))=nankind) or (abs(kind(y))=nankind) then begin
picknan( x, y, q ) ;
r := q ;
end
else
        begin (* no nan *)
        donormalize(x) ; (* Rem always normalizes x. *)
        case abs(kind(y)) of

zerokind, unnormkind : begin
makenan( nanrem, q ) ;
r := q ;
end ;

normkind : case abs(kind(x)) of
zerokind : begin
makezero(q) ;
r := x ;
end ;
unnormkind, normkind : dodivrem ;
infkind : begin
q := x ;
makenan( nanrem, r) ;
end ;
end ;

infkind : case abs(kind(x)) of
;zerokind, normkind : begin
makezero(q) ;
r := x ;
end ;
infkind : begin
q := x ;
makenan( nanrem, r) ;
end end end ;
        end   (* not nan *) ;
q.sign := x.sign <> y.sign ; (* EXOR. *)
end ;


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