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

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

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

Mod.sources:  Volume 3, Issue 7
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 extra.i
cat >extra.i <<'End-Of-File'

(* File extra.i, version 9 October 1984 *)

procedure csqrt ( x : internal ; var z : internal  ) ;

        (* Computes z := sqrt(x).  *)

procedure dosqrt ;

        (* Does SQRT for normalized positive x.  *)

var
i, j : integer ;
r : internal ;
carry : boolean ;
sbit, vbit, orbit : boolean ;

begin
roundkcs ( x, fpstatus.mode.round, xprec ) ; (* Pre-round.  *)
r := x ; (* R will be the remainder for the nonrestoring binary square root  *)
z.sign := false ; (* Result is never negative since x is positive donormalize  *)
if odd(r.exponent) then begin
r.exponent := r.exponent + 1 ; (* Make exponent even.  *)
right( r, 1 ) ; (* And make fraction 0.25 <= r <= 0.5  *)
end ;
z.exponent := r.exponent div 2 ;
sbit := false ; (* Sign bit of remainder, initially positive.  *)
carry := false  ;
                (* Subtract 0.25 to start the fun.  *)
suber(r.significand[1], true, r.significand[1], carry) ;
suber(r.significand[0], false, r.significand[0], carry) ;

                (* Now do main loop.
                Ri fits in i+1 bits.
                Zi fits in i-1 bits.  *)

for i := 1 to (leastsigbit+2) do
if sbit then begin (* R is negative so add:
                        Zi+1 := 2 Zi
                        Ri+1 := 4 Ri + 4 Zi+1 + 3    *)
z.significand[i-1] := false ; (* Set result bit.  *)
vbit := r.significand[0] ;  (* Catch overfl.  *)
left(r,1) ; (* Multiply R by 2.  *)
carry := false ;
adder( r.significand[i+1], true, r.significand[i+1], carry) ;
(* Add 3*2**-i-2 *)
adder(r.significand[i], true, r.significand[i], carry) ;
for j := (i-1) downto 0 do (* Add Zi+1.  *)
adder(r.significand[j], z.significand[j], r.significand[j], carry ) ;
adder( vbit, false, vbit, carry ) ;
adder ( sbit, false, sbit, carry ) ; (* Sets sign of r.  *)
end

else begin (* R is >= 0 so subtract:
                        Zi+1 := 2 Zi + 1
                        Ri+1 := 4 Ri - 4 Zi+1 - 1   *)
z.significand[i-1] := true ; (* Set result bit.  *)
vbit := r.significand[0] ;
left(r,1) ;
carry := false ;
suber( r.significand[i+1], true, r.significand[i+1], carry ) ;
        (* Subtract 1 *)
suber(r.significand[i], false, r.significand[i], carry ) ;
for j := (i-1) downto 0 do (* Subtract Zi+1 *)
suber( r.significand[j], z.significand[j], r.significand[j], carry ) ;
suber( vbit, false, vbit, carry ) ;
suber( sbit, false, sbit, carry ) ;
end ;

z.significand[stickybit-1] := false ; (* This bit isn't used.  *)

        (* Determine sticky bit.  Z is exact iff
        Rn + 4 Zn + 1 <= 0   *)

carry := false ; orbit := false ;
adder( r.significand[leastsigbit+3], true, vbit, carry ) ; (* Add 1.  *)
orbit := orbit or vbit ;
adder( r.significand[leastsigbit+2], false, vbit, carry ) ;
orbit := orbit or vbit ;
for j := (leastsigbit+1) downto 0 do begin
adder( r.significand[j], z.significand[j], vbit, carry ) ;
orbit := orbit or vbit ;
end ;
adder(sbit, false, vbit, carry ) ;
orbit := orbit or vbit ;
adder( sbit, false, sbit, carry ) ;
z.significand[stickybit] := orbit and (not sbit) ;
                (* Inexact if result of test is positive.  *)
end ;


begin (* csqrt*)

case kind(x) of
negnan, nankind : z := x ;
neginf, negnorm, negunnorm, unnormkind : makenan(nansqrt, z) ;
zerokind : z := x ;
normkind : dosqrt ;
infkind : if fpstatus.mode.clos = affine then z := x else
makenan(nansqrt, z ) ;
otherwise
end ;

end ;

procedure clogb ( x : internal ; var z : internal ) ;

        (* Sets y to the unbiased exponent of x.  *)

var
yi : cint64 ;
i, k : integer ;

begin
case abs(kind(x)) of

zerokind : begin
makeinf(z) ;
z.sign := true ;
end ;

unnormkind, normkind : begin
for i := 0 to 5 do yi[i] := 0 ;
k := x.exponent - 1 ; (* -1 because binary point is to left of bit 0.  *)
yi[6] := abs(k) div 256 ;
yi[7] := abs(k) mod 256 ;
unpackinteger ( yi, z, i16 ) ;
z.sign := k < 0 ;
end ;

infkind : begin
makeinf(z) ;
z.sign := false ;
end ;

nankind : z := x ;
otherwise
end ;
end ;

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

        (* Sets z to the next machine number after x in the direction of
        y.  *)

var
cc : conditioncode ;
i : integer ;
rnd : roundtype ;
moveright : boolean ;
t : internal ;

begin
roundkcs(x, fpstatus.mode.round, xprec ) ; (* Preround.  *)
roundkcs(y, fpstatus.mode.round, xprec ) ;
z := x ; (* Default result.  *)
compare( x, y, cc ) ;
if cc in [lesser,greater] then
        begin (* x <> y *)
        moveright := cc = lesser ; (* If x < y then move x to right (+INF) *)
        rnd := fpstatus.mode.round ;
        if moveright then fpstatus.mode.round := rpos else
                fpstatus.mode.round := rneg ;
        case abs(kind(x)) of
        zerokind : begin (* zero *)
                   z.significand[leastsigbit] := true ;
                   z.sign := not moveright ;
                   end   (* zero *) ;
        infkind :  begin (* inf *)
                   z.exponent := maxexp - 1 ;
                   for i := 0 to leastsigbit do z.significand[i] := true ;
                   z.sign := moveright ;
                   end   (* inf *) ;
        unnormkind, normkind :
                if unzero(x) then z.exponent := x.exponent - 1
        else
                begin (* Do add *)
                makezero(t) ;
                t.significand[leastsigbit] := true ;
                t.sign := not moveright ;
                add(x, t, z) ;
                end   (* Do add *) ;
        otherwise
    end (* case *) ;
        roundkcs( z, fpstatus.mode.round, fpstatus.mode.precision ) ;
        store(z) ;
     fpstatus.mode.round := rnd ; (* Force special rounding mode on store.  *)
        end   (* x <> y *) ;
fpstatus.curexcep := fpstatus.curexcep - [inxact] ; (* Don't want inxact
        on a NEXT operation.  *)
end ;


procedure complement ( var x : internal ; var v : boolean ) ;

        (* Complements x.significand, treating it as a 64 bit integer.
        v is a carry out bit.  *)

var
carry : boolean ;
i : integer ;

begin
carry := false ;
for i := leastsigbit downto 0 do
suber( false, x.significand[i], x.significand[i], carry ) ;
v := carry ;
end ;

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

        (* Sets z to x * 2 **int(y).  *)

var
rx, ry : roundtype ;

procedure doscale ;

        (* Carries out scaling for proper x and y.  *)

var
xe : internal ;
i, k : integer ;
v, v2, carry : boolean ;
s : strng ;
irs : integer ;

begin
z := x ; (* Now all we have to do is set the exponent.  *)
xe.sign := x.exponent < 0 ; (* xe will contain exponent of x expanded.  *)
k := abs(x.exponent) ;
for i := leastsigbit downto 0 do begin
xe.significand[i] := odd(k) ;
k := k div 2   ;
end ;

if xe.sign then complement( xe, v2 ) ;

if y.exponent > 64 then begin (* Substitute for huge y.  *)
y.exponent := 64 ;
y.significand[0] := true ;
end ;
if y.exponent < (64-stickybit) then irs := stickybit
        (* Look out for 16 bit integer overfl.  *)
else irs := 64 - y.exponent ; (* Set up count for right shift.  *)
right( y, irs  ) ; (* Align significand of y as an integer.  *)
if y.sign then complement(y, v) ;
carry := false ;
for i := stickybit downto 0 do
adder( xe.significand[i], y.significand[i], xe.significand[i], carry ) ;
adder( v, v2, xe.sign, carry ) ;

if xe.sign then complement( xe, v ) ;

v := not zerofield( xe, 0, 48 ) ; (* v is now an overfl flag.  *)
k := 0 ;
for i := 49 to leastsigbit do begin
k := k + k ;
if xe.significand[i] then k := k + 1 ;
end ;
if xe.sign then k := -k ; (* Set up correct negative exponent.  *)
v := v or (k=maxexp) or (k=minexp) ;
if v then begin (* Exponent overfl.  *)
if xe.sign then begin (* Floating underfl.  *)
makezero(z) ;
setex ( underfl ) ;
end
else begin (* Floating overfl.  *)
makeinf(z) ;
setex ( overfl ) ;
end
end
else z.exponent := k ;
end ;

begin (* Scale.  *)
if (abs(kind(x))=nankind) or (abs(kind(x))=nankind) then
picknan(x, y, z ) else begin
rx := fpstatus.mode.round ; (* Default.  *)
ry := rx ;
case rx of
rneg : if x.sign then ry := rpos ;
rpos : if x.sign then ry := rneg ;
rzero : ry := rneg ;
otherwise
end ;

roundkcs(x, rx, xprec) ;
roundint(y, ry, xprec) ;
donormalize(y) ;

case abs(kind(x)) of

zerokind : case abs(kind(y)) of

zerokind, normkind : z := x ;
infkind : if (fpstatus.mode.clos = affine) and
(kind(y) = neginf) then z := x else
makenan( nanmul, z) ; (* 2 **INF = NAN,  2**+INF = +INF, 2**-INF = 0 *)
end ;

unnormkind, normkind : case abs(kind(y)) of
zerokind, normkind : doscale ;
infkind : if fpstatus.mode.clos = proj then makenan(nanmul, z)
else if x.sign then makezero(z)
else makeinf(z) ;
end ;

infkind : case abs(kind(y)) of
zerokind, normkind : z := x ;
infkind : if (fpstatus.mode.clos=proj) or (kind(x)=neginf) then
makenan(nanmul, z)
else z := x ;
end ;

otherwise
end ;
z.sign := x.sign ;
end ;
end ;


End-Of-File
echo Extracting storage.i
cat >storage.i <<'End-Of-File'
(* File storage.i, Version 9 October 1984.  *)

function xbyte ( x : internal ; p1, p2 : integer ) : BYT ;

        (* Converts bits
        x.significand[p1..p2]
        into a BYT value.  *)

var
b : BYT ;
i : integer ;

begin
b := 0 ;
for i := p1 to p2 do
if x.significand[i] then b := b + b + 1 else b := b + b ;
xbyte := b ;
end ;

procedure ibytes ( k : integer ; var b1, b2 : BYT ) ;

        (* Converts 16 bit integer into two BYT values.  *)

var
neg : boolean ;

begin
neg := k < 0 ;
if neg then k := ( k + 16384 ) + 16384  ; (* Remove most significant bit.  *)
b1 := k div 256 ;
b2 := k mod 256 ;
if neg then b1 := b1 + 128 ; (* Restore most significant bit.  *)
end ;

procedure bytehex ( b : BYT  ; var s  : strng )  ;

        (* Converts BYT to two hex digits.  *)

var
nib : nibarray ;
i,j : integer ;
w : BYT ;

begin
s[0] := chr(2) ;
w := b ;
for j := 2 downto 1 do begin
for i := 3 downto 0 do begin
nib[i] := odd(w) ;
w := w div 2 ;
end ;
s[j] := nibblehex(nib) ;
end ;
end ;

procedure bytex ( b : BYT ; var x : internal ; p1, p2 : integer ) ;

        (* Inserts BYT b into
        x.significand[p1..p2] *)

var i : integer ;

begin
for i := p2 downto p1 do begin
x.significand[i] := odd(ord(b)) ;
b := b div 2 ;
end ;
end ;

procedure unpackextended ( y : cextended ; var x : internal ) ;

        (* Unpacks cextended into internal.  *)

var
zero : boolean ;
i : integer ;

begin
x.sign := (y[0] >= 128) ;
if x.sign then y[0] := y[0] - 128 ; (* Remove sign bit.  *)
x.exponent := (256*y[0] + y[1]) - biasex ;
for i := 2 to 9 do bytex( y[i], x, (8*i-16), (8*i-9) ) ;
for i := (leastsigbit+1) to stickybit do
x.significand[i] := false ;

if x.exponent >= maxex then x.exponent := maxexp ; (* INF/NAN *)
if x.exponent <= minex then begin
zero := y[2]=0 ;
for i := 3 to 9 do
zero := zero and (y[i]=0) ;
if zero then x.exponent := minexp else begin
x.exponent := minex + 1 ;
        (* Add offset for cextended denormalized.  *)
if (fpstatus.mode.norm = normalizing)  then begin
donormalize(x) ;
end
        (* Normalize denormalized operand in normalizing mode.  *)
end
end ;

end ;

procedure toextended ( var x : internal ; var y : cextended  ) ;

        (* Converts x to cextended y.  *)

var i : integer ;
s : strng ;
special : boolean ;
y0,y1 : BYT ;

begin
case abs(kind(x)) of
otherwise ;
zerokind : x.exponent := minex ;

unnormkind, normkind : begin
if x.exponent <= minex then begin (* Underflow.  *)
if underfl in fpstatus.trap then begin (* Trap enabled.  *)
setex ( underfl ) ;
x.exponent := x.exponent + 24576 ;
if x.exponent <= minex then begin (* Severe underfl - give invalid result. *)
makenan(nanresult,x) ;
end ;
end
else begin (* Trap disabled.  *)
right( x, minex + 1   - x.exponent ) ;
x.exponent := minex ;
roundkcs( x, fpstatus.mode.round, fpstatus.mode.precision ) ;
if inxact in fpstatus.curexcep  then
        setex ( underfl ) ; (* Signal.  *)
end ;
end ;

roundkcs( x, fpstatus.mode.round, fpstatus.mode.precision ) ;
if (x.exponent >= maxex)  then begin (* Overflow.  *)
if overfl in fpstatus.trap then begin (* Trap enabled.  *)
setex ( overfl ) ;
x.exponent := x.exponent - 24576 ;
if x.exponent >= maxex then (* Severe overfl - give invalid result.  *)
begin
makenan(nanresult,x) ;
end ;
end

else begin (* Trap disabled.  *)
setex ( inxact ) ;
setex( overfl ) ;
case fpstatus.mode.round of
rneg : special := not x.sign ;
rpos : special := x.sign ;
rnear : special := false ;
rzero : special := true ;
otherwise
end ;
if special then begin (* Special case roundings.  *)
x.exponent := maxex - 1 ;
        (* Round normalized to largest normalized number.
        Round unnormalized to largest exponent, same significand.  *)
if x.significand[0] then
for i := 0 to leastsigbit do x.significand[i] := true ;
end
else begin (* Normal case - set INF.  *)
x.exponent := maxex ;
for i := 0 to leastsigbit do x.significand[i] := false ;
end ;
end
end ;
if abs(kind(x)) = nankind then begin
setex(invop) ;
fpstatus.curexcep := fpstatus.curexcep - [inxact] ;
x.exponent := maxex ;
end end ;

infkind, nankind : x.exponent := maxex ;
end ;

for i := 2 to 9 do (* Pack significand.  *)
y[i] := xbyte ( x, (8*i-16), (8*i-9) ) ;
ibytes ( x.exponent + biasex, y0, y1 ) ; (* Pack exponent.  *)
y[0] := y0 ; y[1] := y1 ;
if x.sign then y[0] := y[0] + 128 ; (* Pack sign bit.  *)

write(' Extended format: ') ;
for i := 0 to 9 do begin
bytehex( y[i], s ) ;
write(s[1],s[2], ' ') ;
end ;
writeln ;

unpackextended ( y, x) ;

end ;

procedure unpackdouble (* y : cdouble ; var x : internal *) ;

        (* Unpacks cdouble into internal.  *)

var
i : integer ;
zero : boolean ;

begin
x.sign := y[0] >= 128 ;
if x.sign then y[0] := y[0] - 128 ;
x.exponent := (16*y[0] + (y[1] div 16)) - biased ;
bytex ( y[1] mod 16, x, 1, 4 ) ;
for i := 2 to 7 do bytex ( y[i], x, (8*i-11), (8*i-4) ) ;
for i := 53 to stickybit do x.significand[i] := false ;

if x.exponent >= maxed then begin
x.exponent := maxexp ;
x.significand[0] := false ;
end
else if x.exponent <= mined then  begin
x.significand[0] := false ;
if zerofield( x, 1, 52 ) then x.exponent := minexp (* Normal Zero.  *)
else x.exponent := x.exponent + 1 ; (* Offset for denormalized numbers.  *)
if (fpstatus.mode.norm = normalizing)  then donormalize(x)
        (* Normalize denormalized operand in normalizing mode.  *)
end
else x.significand[0] := true ; (* Insert leading bit. *)
end ;

procedure todouble (* var x : internal ; var y : cdouble  *) ;

        (* Converts x to cdouble y.  *)

var
i : integer ;
s : strng ;
special : boolean ;
y0,y1 : BYT ;

begin
case abs(kind(x)) of
otherwise ;
zerokind : x.exponent := mined ;

unnormkind, normkind : begin
if x.exponent <= mined then begin (* Underflow.  *)
if underfl in fpstatus.trap then begin (* Trap enabled.  *)
setex ( underfl ) ;
x.exponent := x.exponent + 1536 ;
if( x.exponent <= mined) or not x.significand[0]  then begin (* Severe underfl.  *)
makenan(nanresult,x)
end ;
end
else begin (* Trap disabled.  *)
right( x, mined + 1 - x.exponent ) ;
x.exponent := mined+1 ;
roundkcs( x, fpstatus.mode.round, dprec ) ;
if inxact in fpstatus.curexcep  then setex ( underfl ) ; (* Signal.  *)
end ;
end ;

roundkcs( x, fpstatus.mode.round, dprec ) ;
if (x.exponent >= maxed) and x.significand[0] then begin (* Overflow.  *)
if overfl in fpstatus.trap then begin (* Trap enabled.  *)
setex ( overfl ) ;
x.exponent := x.exponent - 1536 ;
if x.exponent >= maxed then begin (* Severe overfl.  *)
makenan(nanresult,x)
end ;
end

else begin (* Trap disabled.  *)
setex ( inxact ) ;
setex( overfl ) ;
case fpstatus.mode.round of
rneg : special := not x.sign ;
rpos : special := x.sign ;
rnear : special := false ;
rzero : special := true ;
otherwise
end ;
if special then begin (* Special case roundings.  *)
x.exponent := maxed - 1 ; (* Round to largest normalized number.  *)
for i := 0 to leastsigbit do x.significand[i] := true ;
end
else begin (* Normal case - set INF.  *)
x.exponent := maxed ;
for i := 0 to leastsigbit do x.significand[i] := false ;
end ;
end
end ;

if (x.exponent=(mined+1)) and (not x.significand[0]) then
        x.exponent := mined ; (* Look for  denormalized number,
         which may have resulted from an underfl, but might not have.  *)

if (abs(kind(x))=nankind) or (  (x.exponent > mined) and (x.exponent < maxed)
and not x.significand[0]) then begin
(* Invalid Result.  *)
makenan( nanresult, x ) ;
setex ( invop ) ;
fpstatus.curexcep := fpstatus.curexcep - [  inxact ] ;
x.exponent := maxed ;
end ;
end ;

infkind, nankind :
        begin (* inf/nan *)
        x.exponent := maxed ;
        for i := 53 to leastsigbit do
                if x.significand[i] then x.significand[52] := true ;
                (* OR together least significant bits of NAN *)
        end   (* inf/nan *) ;
end (* case *);

ibytes (( x.exponent + biased) * 16, y0, y1 ) ;
        (* Pack exponent *)
y[0] := y0 ; y[1] := y1 ;
if x.sign then y[0] := y[0] + 128 ; (* Pack sign.  *)
y[1] := y[1] + xbyte( x, 1, 4 ) ;
for i := 2 to 7 do
y[i] := xbyte ( x, 8 * i - 11, 8 * i - 4 ) ; (* Pack significand.  *)

write(' Double format: ') ;
for i := 0 to 7 do begin
bytehex( y[i], s ) ;
write(s[1],s[2], ' ') ;
end ;
writeln ;

unpackdouble( y, x ) ;
end ;

procedure unpacksingle (* y : csingle ; var x : internal *) ;

        (* Unpacks csingle into internal.  *)

var
i : integer ;
zero : boolean ;

begin
x.sign := y[0] >= 128 ;
if x.sign then y[0] := y[0] - 128 ;
x.exponent := (2*y[0] + (y[1] div 128)) - biases ;
bytex ( y[1] mod 128, x, 1, 7 ) ;
for i := 2 to 3 do bytex ( y[i], x, (8*i-8), (8*i-1) ) ;
for i := 24 to stickybit do x.significand[i] := false ;

if x.exponent >= maxes then begin
x.exponent := maxexp ;
x.significand[0] := false ;
end
else if x.exponent <= mines then  begin
x.significand[0] := false ;
if zerofield( x, 1, 23 ) then x.exponent := minexp (* Normal Zero.  *)
else x.exponent := x.exponent + 1 ; (* Offset for denormalized numbers.  *)
if (fpstatus.mode.norm = normalizing)   then donormalize(x)
        (* Normalize denormalized operand in normalizing mode.  *)
end
else x.significand[0] := true ; (* Insert leading bit. *)
end ;

procedure tosingle (* var x : internal ; var y : csingle  *) ;

        (* Converts x to csingle y.  *)

var
i : integer ;
s : strng ;
special : boolean ;
y0,y1 : BYT ;

begin
case abs(kind(x)) of
otherwise ;
zerokind : x.exponent := mines ;

unnormkind, normkind : begin
if x.exponent <= mines then begin (* Underflow.  *)
if underfl in fpstatus.trap then begin (* Trap enabled.  *)
setex ( underfl ) ;
x.exponent := x.exponent + 192 ;
if (  x.exponent <= mines) or (not x.significand[0])
then begin (* Severe underfl.  *)
makenan(nanresult,x) ;
end ;
end
else begin (* Trap disabled.  *)
right( x, mines + 1 - x.exponent ) ;
x.exponent := mines+1 ;
roundkcs( x, fpstatus.mode.round, sprec ) ;
if inxact in fpstatus.curexcep  then setex ( underfl ) ; (* Signal.  *)
end ;
end ;

roundkcs( x, fpstatus.mode.round, sprec ) ;
if (x.exponent >= maxes) and x.significand[0] then begin (* Overflow.  *)
if overfl in fpstatus.trap then begin (* Trap enabled.  *)
setex ( overfl ) ;
x.exponent := x.exponent - 192 ;
if x.exponent >= maxes then begin (* Severe overfl.  *)
makenan(nanresult,x) ;
end ;
end

else begin (* Trap disabled.  *)
setex ( inxact ) ;
setex( overfl ) ;
case fpstatus.mode.round of
rneg : special := not x.sign ;
rpos : special := x.sign ;
rnear : special := false ;
rzero : special := true ;
otherwise
end ;
if special then begin (* Special case roundings.  *)
x.exponent := maxes - 1 ; (* Round to largest normalized number.  *)
for i := 0 to leastsigbit do x.significand[i] := true ;
end
else begin (* Normal case - set INF.  *)
x.exponent := maxes ;
for i := 0 to leastsigbit do x.significand[i] := false ;
end ;
end
end ;
if  ( (x.exponent=(mines+1)) and (not x.significand[0]))
then
        x.exponent := mines ; (* Look for  denormalized number.  *)

if (abs(kind(x))=nankind) or (  (x.exponent > mines) and (x.exponent < maxes)
and not x.significand[0] )  then begin
        (* Invalid Result.  *)
        makenan( nanresult, x ) ;
        setex ( invop ) ;
        fpstatus.curexcep := fpstatus.curexcep - [inxact] ;
        x.exponent := maxes ;
        end ;
        end ;

infkind, nankind :
        begin (* inf/nan *)
        x.exponent := maxes ;
        for i := 24 to leastsigbit do
                if x.significand[i] then x.significand[23] := true ;
                (* OR together least significant bits of NAN *)
        end   (* inf/nan *) ;
end (* case *);

ibytes (( x.exponent + biases) * 128, y0, y1 ) ;
        (* Pack exponent *)
y[0] := y0 ; y[1] := y1 ;
if x.sign then y[0] := y[0] + 128 ; (* Pack sign.  *)
y[1] := y[1] + xbyte( x, 1, 7 ) ;
for i := 2 to 3 do
y[i] := xbyte ( x, 8 * i - 8 , 8 * i - 1 ) ; (* Pack significand.  *)

write(' Single format: ') ;
for i := 0 to 3 do begin
bytehex( y[i], s ) ;
write(s[1],s[2], ' ') ;
end ;
writeln ;

unpacksingle( y, x ) ;
end ;

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

        (* Unpacks integer in y according to itype.
        The significant bytes are presumed to be on the right.  *)

var i, msy : integer ;
carry : boolean ;
es : excepset ;

begin
case itype of
i16 : msy := 6 ;
i32 : msy := 4 ;
i64 : msy := 0 ;
otherwise
end ;
x.sign := y[msy] >= 128 ;
if x.sign then (* Expand negative.  *)
for i := 0 to (msy-1) do y[i] := 255
else
for i := 0 to (msy-1) do y[i] := 0 ;
for i := 0 to 7 do bytex( y[i], x, 8*i, 8*i+7) ;
if x.sign then begin
carry := false ;
for i:= leastsigbit downto 0 do
suber( false, x.significand[i], x.significand[i], carry ) ;
end ;
for i := (leastsigbit+1) to stickybit do x.significand[i] := false ;
x.exponent := 64 ;
donormalize(x) ;
if (itype = i64) and (x.exponent = 64) then
        begin (* It was really a NAN *)
        es := fpstatus.curexcep ;
        makenan(naninteger, x) ;
        x.sign := false ; (* Default is a positive NAN.  *)
        fpstatus.curexcep := es ; (* Don't let makenan set NV.  *)
        end   (* It was really a NAN *) ;
end ;

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

        (* Converts x into integer value of type i-type.  *)

var
i, imax : integer ;
s : strng ;
carry : boolean ;

procedure i64nan ;
        (* Creates an int64 nan *)
var i : integer ;
begin (* i64nan *)
x.significand[0] := true ;
for i := 1 to stickybit do x.significand[i] := false ;
end   (* i64nan *) ;

begin
case itype of
i16 : imax := 16 ;
i32 : imax := 32 ;
i64 : imax := 64 ;
otherwise
end ;

case abs(kind(x)) of
otherwise ;
unnormkind, normkind : begin
roundint( x, fpstatus.mode.round, xprec) ;
donormalize(x) ;
if kind(x) <> zerokind then begin
if x.exponent < 64 then right( x, 64 - x.exponent ) ;
if x.exponent > 64 then
        begin
        left ( x, x.exponent - 64 ) ;
        end ;
if (x.exponent >= imax) and (* Exclude case of max negative integer.  *)
((x.exponent <> imax) or (not x.sign) or
(lastbit(x,leastsigbit-imax+1,leastsigbit) > (leastsigbit-imax+1)))
then begin
x.significand[leastsigbit+1-imax] := false ; (* Turn off bit to allow room
        for sign bit.  *)
setex ( cvtovfl ) ;
end ;
if (itype=i64) and (x.exponent >= imax) then
        begin (* overflowed to nan *)
        i64nan ;
        setex(cvtovfl) ; (* Might not have been set for -2^63.  *)
        end   (* overflowed to nan *) ;
end
end
 ;

infkind : begin
setex ( cvtovfl  ) ;
if itype = i64 then i64nan else
        begin (* not i64 *)
        for i := leastsigbit downto (leastsigbit - imax + 2 )
        do x.significand[i] := true ;
        x.significand[leastsigbit-imax+1] := false ;
        end   (* not i64 *) ;
end ;
nankind : begin
if itype = i64 then i64nan else
        begin (* not i64 *)
        setex ( invop ) ;
        for i := leastsigbit downto (leastsigbit - imax + 2  )
        do x.significand[i] := false ;
        x.significand[leastsigbit-imax+1] := true ;
        end   (* not i64 *) ;
end ;

end ;

if x.sign then begin (* Complement.  *)
carry := false ;
for i := leastsigbit downto (leastsigbit - imax + 1) do
suber( false, x.significand[i], x.significand[i], carry ) ;
end ;

for i := 0 to 7 - (imax div 8) do y[i] := 0 ;
for i := (8 - (imax div 8)) to 7 do
y[i] := xbyte( x, leastsigbit - 63 + 8*i, leastsigbit - 56 + 8*i ) ;

write(' Integer format: ') ;
for i := (8 - (imax div 8)) to 7 do  begin
bytehex(y[i],s) ;
write(s[1],s[2],' ') ;
end ;
writeln ;

unpackinteger( y, x, itype ) ;

end ;


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

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

var
ztrue, z, r : internal ;
cc : conditioncode ;
ps : pstack ;
error : boolean ;
i, k: integer ;
yi : cint64 ;
ms : fpmodetype ; es, ts : excepset ;

procedure subRR ;

begin
if sequal(s , 'REM') then begin
found := true ;
trem( y, x,  z ) ;
end
end ;

procedure subS ;

var
xr,yr,zr :real ;

begin
if sequal(s , 'SCALE') then begin
found := true ;


cscale( y, x,  z ) ;

end else if sequal(s , 'SQRT') then begin
found := true ;

tsqrt( x, z) ;

end
end ;

procedure subT ;

var yi : cint64 ;

begin
if sequal(s , 'TEST') then begin
found := true ;
pretest( storagemode )  ;
end
else if sequal(s , 'TOF32') then begin (* Convert to single.  *)
found := true ;
tconvert(x,z,flt32) ;
end else if sequal(s , 'TOF32I') then begin (* Convert to single integral.  *)
found := true ;
tintconvert(x,z,flt32) ;
end else if sequal(s , 'TOF64') then begin (* Convert to double.  *)
found := true ;
tconvert(x,z,f64) ;
end else if sequal(s , 'TOF64I') then begin (* Convert to double integral.  *)
found := true ;
tintconvert(x,z,f64) ;
end else  if sequal(s , 'TOX80') then begin (* Convert to extended.  *)
found := true ;
tconvert(x,z,ext80) ;
end else if sequal(s , 'TOX80I') then begin (* Convert to extended integral.  *)
found := true ;
tintconvert(x,z,ext80) ;
end else if sequal(s , 'TOI16') then begin (* Convert to 16 bit integer.  *)
found := true ;
tconvert(x,z,i16) ;
end else if sequal(s , 'TOI32') then begin (* Convert to 32 bit integer.  *)
found := true ;
tconvert(x,z,i32) ;
end else if sequal(s , 'TOI64') then begin (* Convert to 64 bit integer.  *)
found := true ;
tconvert(x,z,i64) ;
end  ;
end ;


begin
writeln(' BEGIN TEST ') ;
makezero(z) ; (* Define default "computed result" for those operations
        that don't return any.  *)
if stack = nil then makezero(ztrue) else ztrue := stack^.x ;
if not sequal(s,'TEST') then begin (* Not ready to do these mode switches until
initialization has been accomplished.  *)

ms := fpstatus.mode ;
swapmode(ms) ;
ts := fpstatus.trap ;
swaptrap(ts) ;
es := fpstatus.excep ;
swapexcep(es) ;
end ;
found := false ;
if length(s) > 0 then case s[1] of

'+' : if length(s)=1 then begin
found := true ;

tadd( y, x,  z ) ;

end ;

'-' : if length(s)=1 then begin
found := true ;

tsub( y, x,  z ) ;

end ;

'*' : if length(s)=1 then begin
found := true ;
tmul (y, x, z) ;

end ;

'/' : if length(s) = 1 then begin
found := true ;
tdiv ( y, x,  z) ;

end ;
'A' : if sequal(s , 'ABS') then begin
found := true ;
tabs(x,z) ;
end
 ;

'C' : if sequal(s , 'COMPARE') then begin
found := true ;
tcompare( y, x,  cc) ;
write(' Compare result: ') ;
case cc of
lesser : writeln(' < ') ;
equal : writeln(' = ' ) ;
greater : writeln(' > ') ;
notord : writeln(' Unordered ') ;
end ;
for i := 0 to 6 do yi[i] := 0 ;
yi[7] := ord(cc) ;
unpackinteger(yi, z, i16);
end ;

'L' : if sequal(s , 'LOGB') then begin
found := true ;
clogb( x,  z ) ;
end ;

'N' : if sequal(s , 'NEG') then begin (* NEGATE top of stack *)
found := true ;
tneg(x,z) ;
end
else if sequal(s , 'NEXT') then begin (* Compute NEXTAFTER function.  *)
found := true ;
cnextafter( y, x,  z ) ;

end ;

'R' : subRr ;
'S' : subS ;
'T' : subT ;

otherwise

end ;

if found then writeln( ' Did ',s) ;

if not found then begin (* check for decimal input *)
tdecbin(s, z, error ) ;
if not error then begin
found := true ;

end
end ;
if sequal(s,'TEST') then writeln(' Begin TEST Mode ')
else begin
if  found then begin
tstore(storagemode,z) ;
swapexcep(es) ;
if (es=fpstatus.excep) and (equalinternal(z,ztrue)) then
writeln(' OK! ')
else
begin
if es <> fpstatus.excep then
        begin
        write(chr(ordbell),' DIFFERENT FLAGS: ') ;
        displayexcep(es) ;
        writeln ;
        end ;
if not equalinternal( z, ztrue ) then
        begin
        writeln(chr(ordbell),' DIFFERENT RESULT: ') ;
        display(z) ;
        end ;
end ;
tdisplay(z) ;
writeln(' END TEST  ') ;
end
else  writeln(' Command not tested: ',s) ;
end ;
end ;



End-Of-File
echo Extracting hex.i
cat >hex.i <<'End-Of-File'
(* File hex.i, Version 8 October 1984 *)

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

                (* Interprets s as a hex 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 ;
nib : nibarray ;

begin
error := false ;
for i := p1 to p2 do x.significand[i] := false ; (* Clear field.  *)
i := p2 + 1 - 4 * length(s) ;
while i < p2 do begin
hexnibble( s[1], nib ) ;
delete ( s, 1, 1 ) ;
for j := 0 to 3 do if nib[j] then begin
if (i+j) < p1 then error := true else x.significand[i+j] := true ;
end ;
i := i + 4 ;
end ;
end ;

procedure intdec ( i : integer ; var s : strng ) ;
        (* converts 16 bit integer to decimal strng *)
var
sign : boolean ;
t : strng ;

begin
if i = 0 then
    begin
    s[0] := chr(1) ;
    s[1] := '0' ;
    end
    else begin
t[0] := chr(1) ;
s[0] := chr(0) ;
sign := false ;
if i < 0 then if i < -32767 then begin
makeucsdstring(' -32768',s) ; i := 0 end
else begin
sign := true ; i := -i end ;
while i <> 0 do begin
t[1] := chr( ord('0') + i mod 10 ) ;
s := concat ( t, s ) ;
i := i div 10 ;
end ;
if sign then
    begin
    t[1] := '-' ;
    s := concat( t, s ) ;
    end ;
end
end ;

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

begin
i := p1 ;
while ( i < p2 ) and not x.significand[i] do i := i + 1 ;
        (* Find most significant non-zero bit in field.  *)
if ( i >= p2 ) and not x.significand[p2] then
    begin
    s[0] := chr(1) ;
    s[1] := '0' ;
    end
    else begin
s[0] := chr(0) ;
i := p2 - 3 - 4 * (( p2 - i ) div 4 ) ;
        (* Start at left end of nibarray containing most significant bit.  *)
while i < p2 do begin
for j := 0 to 3 do
if (i+j) < p1 then nib[j] := false else nib[j] := x.significand[i+j] ;
concatchar( s, nibblehex(nib)) ;
i := i + 4 ;
end ;
end ;
end ;

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

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

var
i, npoint : integer ;
nib : nibarray ;
integral : boolean ;
t : strng ;

begin
s[0] := chr(0) ;
if kind(x) = zerokind then
    begin
    s[0] := chr(1) ; s[1] := '0' ;
    end
    else
if (abs(kind(x)) = normkind) and (x.exponent <= 16) and (x.exponent >= 1)
then begin
if zerofield ( x, x.exponent, stickybit ) then begin (* it's all integer *)
subhex ( x, 0, x.exponent - 1, s ) ;
if x.sign then
    begin
    t[0] := chr(1) ;
    t[1] := '-' ;
    s := concat( t, s ) ;
    end ;
end end
end ;

procedure nanascii ( x : internal ; ishex : boolean ; var s : strng ) ;

        (* Converts an INF or NAN into strng s, using hex for numeric
        field values if ishex is true, and decimal if ishex is false.  *)

var t,t1 : strng ;
k : integer ;

begin
case kind(x) of
neginf : makeucsdstring('--',s) ;
infkind : makeucsdstring('++',s) ;
negnan, nankind : begin
makeucsdstring('NaN''',s) ;
if x.sign then
    begin
    t[1] := '-' ;
    s := concat( t, s ) ;
    end ;
if ishex then
        begin (* ishex nan *)
        subhex ( x, 1, 15, t ) ;
        if not zerofield(x,16,leastsigbit) then
                begin (* Extra stuff *)
                concatchar(t,':') ; (* Colon delimits extra stuff.  *)
                for k := 4 to 15 do
                        begin (* Add hexit.  *)
                        subhex(x,4*k,4*k+3,t1) ;
                        t := concat(t,t1) ;
                        end   (* Add hexit.  *) ;
                while t[length(t)] = '0' do
                        delete (t,length(t),1) ; (* Clear trailing zeros. *)
                end   (* Extra stuff *) ;
        end   (* ishex nan *)
else
        if zerofield( x, 1, 15 ) then makeucsdstring('0.',t) else
                begin (* Decimal Nan, non zero *)
                subdec ( x, 1, 15, t ) ;
                concatchar(t,'.') ; (* . Distinguishes decimal NAN from hex *)
                end   (* Decimal Nan, non zero *) ;
s := concat ( s, t) ;
concatchar(s, '''') ;
end ;
otherwise
end ;
end ;

procedure binhex (* x : internal ; var s  : strng *)(* forward *)  ;
(* converts x to hex format *)

var
i, j, k : integer ;
nib : nibarray ;
t : strng ;

begin
case abs(kind(x)) of
zerokind : if x.sign then
    begin
    s[0] := chr(1) ; s[1] := '0' ;
    end
    else
    begin
    s[0] := chr(2) ; s[1] := '-' ; s[2] := '0' ;
    end ;

unnormkind, normkind : begin
tohexint(x, s) ;
if length(s) > 0 then
    begin
    makeucsdstring('H ',t) ; s := concat(s, t) ;
    end
    else
    begin
s[0] := chr(1) ;
s[1] := '.' ;
for i := 0 to 3 do begin
for j := 0 to 3 do begin
for k := 0 to 3 do
nib[k] := x.significand[k+4*j+16*i] ;
concatchar(s, nibblehex(nib)) ;
end ;
concatchar( s, ' ' ) ;
end ;
nib[0] := x.significand[64] ;
nib[1] := x.significand[65] or x.significand[66] ;
nib[2] := false ;
nib[3] := false ;
concatchar(s, nibblehex(nib)) ;

while( (s[length(s)] = ' ') or( s[length(s)] = '0')) and
(length(s) > 2) do delete(s,length(s),1) ; (* delete trailing 0 and blank *)
makeucsdstring('H ',t) ;
s := concat(s,t) ;
if x.exponent <> 0 then begin
if x.exponent > 0 then concatchar(s, '+') ;
intdec(x.exponent, t) ;
s := concat(s,t) ;
end ;
if x.sign then
    begin
    makeucsdstring('- ',t) ;
    s := concat(t,s) ;
    end ;
end end ;

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

otherwise
end ;
end ;

procedure NANer ( s : strng ; ishex : boolean ;
        var x : internal ; var error : boolean ) ;
        (* Checks for strng in proper INF or NAN format.
        If ishex is true, interprets numeric constants in hex;
        If ishex is false, interprets them in decimal.  *)
var
i, k : integer ;
t, snan : strng ;
nminus, ndot, nplus : integer ;
dset : set of char ;
err : boolean ;

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

begin
error := false ;
t[0] := chr(0) ;
for i := 1 to length(s) do if s[i] <> ' ' then concatchar(t,upcase(s[i])) ;
concatchar(t,'z') ;

nminus := 0 ;  nplus := 0 ;
for i := 1 to length(t) do case t[i] of
'-' : nminus := nminus + 1 ;
'+' : nplus := nplus + 1 ;
otherwise
end ;
if (nplus >= 2) and (nplus>=( length(t)-1)) then begin (* plus infinity *)
x.exponent := maxexp ;
makeucsdstring('z ',t) ;
end ;
if (nminus >= 2) and (nminus=( length(t)-1) ) then begin (* minus inf *)
x.exponent := maxexp ;
makeucsdstring('-z',t) ;
end ;
x.sign := t[1]='-' ; (* Check sign *)
if x.sign then bump else if t[1]='+' then bump ;
if (length(t) >= 3)
 then (* check for NAN *)
if (t[1]='N') and (t[2]='A') and (t[3]='N')  then
        begin (* Nan processing *)
        bump ; bump ; bump ;
        x.exponent := maxexp ;
        if t[1]='''' then
                begin (* Process significand string *)
                bump ; (* Remove ' *)
                if ishex then dset := hexset else dset := digitset ;
                snan[0] := chr(0) ;
                while t[1] = '0' do bump ;
                while t[1] in dset do begin (* Accumulate field value. *)
                concatchar( snan, t[1] ) ;
                bump ;
                end ;
                if ishex then
                puthex( snan, 1, 15, x, error )
                else
                putdec( snan, 1, 15, x, error ) ;
                if ishex then
                        begin (* Extra Hex Processing.  *)
                        if t[1] = ':' then
                                begin (* Extra hex stuff *)
                                bump ;
                                k := 16 ;
                                snan[0] := chr(1) ;
                                snan[1] := ' ' ;
                                while (k <= (leastsigbit-3)) and
                                                (t[1] in dset) do
                                        begin
                                        snan[1] := t[1] ;
                                        puthex(snan,k,k+3,x,err) ;
                                        k := k + 4 ;
                                        bump ;
                                        end ;
                                end   (* Extra hex stuff *) ;
                        if t[1]='''' then bump ; (* Absorb final delimiter.  *)
                        end   (* Extra Hex Processing.  *)
                else
                        begin (* Extra Dec Processing *)
                        if t[1]='.' then
                                begin (* Decimal Point Found *)
                                bump ; (* Absorb decimal point.  *)
                                if t[1]='''' then bump ;
                                        (* Absorb final delimiter.  *)
                                end   (* Decimal Point Found *) ;
                        end   (* Extra Dec Processing *) ;
                if length(t) > 1 then
                        begin (* Extra characters *)
                        error := true ;
                        while (length(t)>1) and (t[1]<>'''') do bump ;
                        if t[1]='''' then bump ;
                        end   (* Extra characters *) ;
                end   (* Process significand string *) ;

        if error or zerofield( x, 1, leastsigbit ) then
                begin
                error := false ;
                makenan(nanascnan,x) ;
                (* NAN  format without significand is invalid. *)
                end ;
        end   (* Nan Processing *);
if length(t) > 1 then
        begin
        error := true ;
        end ;
end  (* NANer *) ;

procedure hexbin (* s : strng ; var x : internal ; var error : boolean *) ;
(* converts hex 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 ;
sigpoint : integer ;
t, snan : strng ;
esign : boolean ;
nib : nibarray ;
ee : integer ;

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 ;
for i := 0 to stickybit do x.significand[i] := false ;
sigpoint := 0 ;
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 hexset do begin (* digits before point *)
class := nonzero ;
hexnibble(t[1], nib) ;
if sigpoint <= (stickybit-4) then min := 3 else min := (stickybit-1)-sigpoint ;
for i := 0 to min do x.significand[sigpoint+i] := nib[i] ;
for i := (stickybit-sigpoint) to 3 do x.significand[stickybit] := x.significand[stickybit] or nib[i] ;
x.exponent := x.exponent + 4 ;
if x.significand[0] then begin
if sigpoint <= (stickybit-4) then sigpoint := sigpoint + 4 else sigpoint := stickybit
end else begin (* donormalize *)
donormalize(x) ;
sigpoint := x.exponent ;
end ;
bump
end ;
if t[1] = '.' then begin (* check for point *)
bump ;
while t[1] in hexset do begin (* process digits after point *)
if (t[1] <> '0') or (class = nonzero) then class := nonzero
else class := truezero ;
hexnibble(t[1], nib) ;
if sigpoint <= (stickybit-4) then min := 3 else min := (stickybit-1)-sigpoint ;
for i := 0 to min do x.significand[sigpoint+i] := nib[i] ;
for i := (stickybit-sigpoint) to 3 do
x.significand[stickybit] := x.significand[stickybit] or nib[i] ;
if x.significand[0] then begin
if sigpoint <= (stickybit-4) then sigpoint := sigpoint + 4 else
sigpoint := stickybit
end else if t[1] = '0' then x.exponent := x.exponent - 4 else
begin (* donormalize *)
sigpoint := x.exponent ;
donormalize(x) ;
sigpoint := 4 + x.exponent - sigpoint ;
end ;
bump ;
end ;
end ;
if t[1] = 'H' then bump ; (* handle H for Hex *)
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 ;
if (x.exponent >= 0 ) and (ee > 0 ) then if x.exponent >= (maxexp - ee)
then begin
error := true ;
x.exponent := maxexp - 1 ;
end ;
if (x.exponent < 0) and ( ee < 0 ) then if x.exponent <= (minexp - ee)
then begin
error := true ;
x.exponent := minexp + 1 ;
end ;
if not error then x.exponent := x.exponent + ee ;
end ;
if class = nonnumeric  then
        (* the following code checks for INFs and NANs *)
NANer ( s, true, x, error )
else
if ( length(t) > 1) then error := true  ;
if error then
        begin (* Erroneous input *)
        makenan(nanascbin,x) ;
        end
end ;



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