IEEE Calculator (part 5 of 6)
sources-request at panda.UUCP
sources-request at panda.UUCP
Wed Sep 4 12:12:50 AEST 1985
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
More information about the Mod.sources
mailing list