sources-request@panda.UUCP (09/04/85)
Mod.sources: Volume 3, Issue 4
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 calc.p
cat >calc.p <<'End-Of-File'
program calculator (input, output) ;
(* File calc.p, Version 9 October 1984. *)
(* calc is a calculator style program to demonstrate the proposed
IEEE floating point arithmetic *)
#include 'sane.h'
#include 'oldfplib.h'
#include 'calctest.h'
#include 'calcsingle.h'
#include 'calcdouble.h'
#include 'global.i'
#include 'forward.i'
#include 'init.i'
#include 'divrem.i'
#include 'extra.i'
#include 'storage.i'
#include 'addsubpas.i'
#include 'utility.i'
#include 'function.i'
#include 'hex.i'
#include 'base.i'
procedure store (* var x : internal *) ;
(* Rounds x to current storage mode, setting exceptions accordingly,
then puts result back in internal format. *)
var
yx : cextended ;
yd : cdouble ;
ys : csingle ;
yi : cint64 ;
begin
case storagemode of
i16, i32, i64 : tointeger( storagemode, x, yi ) ;
flt32 : tosingle ( x, ys ) ;
f64 : todouble ( x, yd ) ;
ext80 : toextended ( x, yx ) ;
otherwise
end ;
end ;
procedure commandloop ;
var
c : char ;
s : strng ;
i,j : integer ;
found, exit : boolean ;
ps : pstack ;
badnan, x, y, z, r : internal ;
(* Rule is: x gets the top of stack, y the next,
for use in DOTEST *)
error : boolean ;
cc : conditioncode ;
oldtop : internal ; (* Saves previous top of stack, so we can tell when it
changes. New tops of stack are displayed. *)
heap : ^ integer ; (* Heap marker. *)
yx : cextended ;
yd : cdouble ;
yi : cint64 ;
xs, ys, zs : csingle ;
tx : real ;
es : integer ;
fpe : xcpn ;
buffer : strng ; (* Used to buffer multiple commands. *)
semipos : integer ; (* Used to record end of command. *)
fulldisplay : boolean ; (* Flag set at the end of a calculator operation;
if true, the top of stack will be displayed;
if false, only traps, if any, will be displayed. *)
procedure clear ;
(* Clears stack and heap. *)
begin
stack := nil ;
end ;
procedure docommand ( var found : boolean ) ;
var fpe : xcpn ;
procedure subc ;
var i : integer ;
begin
if sequal(s , 'COMPARE') then begin
found := true ;
popstack(x) ;
popstack(y) ;
compare( 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) ;
pushstack(z) ;
end else
if sequal(s , 'CLEAR') then begin (* CLEAR *) found := true ;
clear end
else if sequal(s , 'CRUNCH') then begin (* Clear stack except for top two entries. *)
found := true ;
popstack(x) ;
popstack(y) ;
clear ;
pushstack(y) ;
pushstack(x) ;
end
;
end ;
procedure subd ;
begin
if sequal(s , 'DUP') then begin (* Duplicate top of stack *)
popstack(x) ;
pushstack(x) ;
pushstack(x) ;
found := true ;
end
else if sequal(s , 'DIV') then begin
found := true ;
popstack(x) ;
popstack(y) ;
divrem( y, x, z, r ) ;
writeln(' REM: ') ;
display(r) ;
pushstack(z) ;
end else if sequal(s, 'DUMP') then begin (* DUMP STACK *)
found := true ;
ps := stack ;
while ps <> nil do begin
display(ps^.x ) ;
ps := ps^.next ;
end ;
end ;
end ;
procedure subRR ;
begin
if sequal(s, 'REV') then begin (* reverse top two entries on stack *)
found := true ;
popstack(x) ;
popstack(y) ;
pushstack(x) ;
pushstack(y) ;
end
else if sequal(s, 'REM') then begin
found := true ;
popstack(x) ;
popstack(y) ;
divrem( y, x, z, r ) ;
writeln(' DIV: ') ;
display(z) ;
pushstack(r) ;
end else if sequal(s, 'RN') then begin
found := true ;
fpstatus.mode.round := rnear ;
end
else if sequal(s, 'RM') then begin
found := true ;
fpstatus.mode.round := rneg ;
end
else if sequal(s, 'RP') then begin
found := true ;
fpstatus.mode.round := rpos ;
end
else if sequal(s, 'RZ') then begin
found := true ;
fpstatus.mode.round := rzero ;
end else if sequal(s, 'R24') then begin
found := true ;
fpstatus.mode.precision := sprec ; (* Round cextended to 24 significant bits. *)
end else if sequal(s, 'R53') then begin
found := true ;
fpstatus.mode.precision := dprec ;
(* Round cextended to 53 significant bits. *)
end end ;
procedure subS ;
begin
{
if sequal(s, 'SOFT') then begin
found := true ;
ffloat_ ; ffunc_ ;
end else if sequal(s, 'SKY') then begin
found := true ;
sfloat_ ; sfunc_ ;
end else }
if sequal(s, 'SCALE') then begin
found := true ;
popstack(x) ;
popstack(y) ;
cscale( y, x, z ) ;
pushstack( z ) ;
end else if sequal(s, 'SQRT') then begin
found := true ;
popstack(x) ;
csqrt( x, z) ;
pushstack(z) ;
end else if sequal(s, 'STOF32') then begin (* Set storage mode. *)
found := true ;
storagemode := flt32 ;
end
else if sequal(s, 'STOF64') then begin
found := true ;
storagemode := f64 ;
end
else if sequal(s, 'STOX80') then begin
found := true ;
storagemode := ext80 ;
end
else if sequal(s, 'STOI16') then begin
found := true ;
storagemode := i16 ;
end else if sequal(s, 'STOI32') then begin
found := true ;
storagemode := i32 ;
end else if sequal(s, 'STOI64') then begin
found := true ;
storagemode := i64
end
end ;
procedure subT ;
var yi : cint64 ;
begin
if sequal(s, 'TOF32') then begin (* Convert to csingle. *)
found := true ;
popstack(x);z:=x ;
tosingle( z, ys) ;
pushstack(z) ;
end else if sequal(s, 'TOF32I') then begin (* Convert to csingle integral. *)
found := true ;
popstack(x);z:=x ;
roundint( z, fpstatus.mode.round, sprec ) ;
tosingle( z, ys) ;
pushstack(z) ;
end else if sequal(s, 'TOF64') then begin (* Convert to cdouble. *)
found := true ;
popstack(x);z:=x ;
todouble( z, yd) ;
pushstack(z) ;
end else if sequal(s, 'TOF64I') then begin (* Convert to cdouble integral. *)
found := true ;
popstack(x);z:=x ;
roundint( z, fpstatus.mode.round, dprec ) ;
todouble( z, yd) ;
pushstack(z) ;
end else if sequal(s, 'TOX80' )then begin (* Convert to cextended. *)
found := true ;
popstack(x);z:=x ;
toextended( z, yx) ;
pushstack(z) ;
end else if sequal(s, 'TOX80I') then begin (* Convert to cextended integral. *)
found := true ;
popstack(x);z:=x ;
roundint( z, fpstatus.mode.round, xprec ) ;
toextended( z, yx) ;
pushstack(z) ;
end else if sequal(s, 'TOI16') then begin (* Convert to 16 bit integer. *)
found := true ;
popstack(x);z:=x ;
tointeger( i16, z, yi) ;
pushstack(z) ;
end else if sequal(s, 'TOI32') then begin (* Convert to 32 bit integer. *)
found := true ;
popstack(x);z:=x ;
tointeger( i32, z, yi) ;
pushstack(z) ;
end else if sequal(s, 'TOI64' )then begin (* Convert to 64 bit integer. *)
found := true ;
popstack(x);z:=x ;
tointeger(i64, z, yi) ;
pushstack(z) ;
end
else if sequal(s, 'TEST') then begin (* Toggle test flag *)
found := true ;
testflag := not testflag ;
end ;
end ;
begin
found := false ;
if length(s) > 0 then case s[1] of
'+' : if length(s)=1 then begin
found := true ;
popstack(x) ;popstack(y) ;
add( y, x, z ) ;
pushstack(z) ;
end ;
'-' : if length(s)=1 then begin
found := true ;
popstack(x) ;popstack(y) ;
r := x ; r.sign := not x.sign ;
add( y, r, z ) ;
pushstack(z) ;
end ;
'*' : if length(s)=1 then begin
found := true ;
popstack(x) ; popstack(y) ;
multiply ( y, x, z) ;
pushstack(z) ;
end ;
'/' : if length(s)=1 then begin
found := true ;
popstack(x) ;popstack(y) ;
divide ( y, x, z) ;
pushstack(z) ;
end ;
'A' : if sequal(s, 'ABS') then begin
found := true ;
popstack(x) ;
z := x ;
z.sign := false ;
pushstack(z) ;
end
else if sequal(s, 'AFF') then begin
found := true ;
fpstatus.mode.clos := affine ;
end ;
'C' : subc ;
'D' : subd ;
'E' : if length(s)=1 then begin
found := true ;
pushstack(e) ;
end else if sequal(s, 'EXIT') then begin (* EXIT *) found := true ; exit := true end ;
'L' : if sequal(s, 'LOGB') then begin
found := true ;
popstack(x) ;
clogb( x, z ) ;
pushstack( z ) ;
end ;
'N' : if sequal(s, 'NEG') then begin (* NEGATE top of stack *)
found := true ;
popstack(x) ;
z := x ;
z.sign := not z.sign ;
pushstack(z) ;
end
else if sequal(s, 'NORM') then begin
found := true ;
fpstatus.mode.norm := normalizing ;
end else if sequal(s, 'NEXT') then begin (* Compute NEXTAFTER function. *)
found := true ;
popstack(x) ;
popstack(y) ;
cnextafter( y, x, z ) ;
pushstack ( z ) ;
end ;
'P' : if sequal(s, 'POP') then begin
found := true ;
if stack <> nil then stack := stack^.next ;
end
else if sequal(s, 'PI') then begin
found := true ;
pushstack(pi) ;
end else if sequal(s, 'PROJ') then begin
found := true ;
fpstatus.mode.clos := proj ;
end ;
'R' : subRr ;
'S' : subS ;
'T' : subT ;
'U' : if sequal(s, 'UNROUNDED') then begin
found := true ;
storagemode := unrounded ;
fpstatus.mode.precision := xprec ;
end ;
'W' : if sequal(s, 'WARN') then begin
found := true ;
fpstatus.mode.norm := warning ;
end ;
otherwise
end ;
if found then writeln( ' Did ',s) ;
if (length(s)=2) and not found then begin (* Is is a trap enable toggle? *)
for fpe := invop to inexact do
if (s[1]=xcpnname[fpe,1]) and (s[2]=xcpnname[fpe,2]) then begin
found := true ; (* Command was name of exception so toggle that trap enable. *)
if fpe in fpstatus.trap then
fpstatus.trap := fpstatus.trap - [fpe] (* If on, turn off. *)
else
fpstatus.trap := fpstatus.trap + [fpe] ; (* If off, turn on. *)
end ;
end ;
if not found then begin (* check for decimal input *)
decbin(s, x, error ) ;
fpstatus.curexcep := fpstatus.curexcep - [invop] ;
badnan.sign := x.sign ; (* Set up BadNaN to compare correctly with x. *)
if (not error) and (not equalinternal(x,badnan)) then begin
found := true ;
pushstack(x) ;
end
end ;
if not found then begin (* check for hex input *)
hexbin(s, x, error ) ;
fpstatus.curexcep := fpstatus.curexcep - [invop] ;
if not error then begin
found := true ;
pushstack(x) ;
end
end ;
if found then begin
if stack <> nil then store(stack^.x) ;
fpstatus.excep := fpstatus.excep + fpstatus.curexcep ; (* Add in current
exceptions. *)
fulldisplay := stack <> nil ;
if fulldisplay then
fulldisplay := (fpstatus.curexcep <> []) or
(not equalinternal( stack^.x, oldtop )) ;
if fulldisplay then begin
displaystatus ;
trapmessage ;
fpstatus.curexcep := [] ;
display(stack^.x) ;
end
else trapmessage
end
else writeln(' Command not recognized: ',s) ;
end ;
begin
clear ;
makenan(nanascnan,badnan) ; (* Create a Bad-NaN Nan for use later. *)
repeat
exit := false ;
fpstatus.excep := [] ;
(* Reset exception flag register for new command strng. *)
writeln(' Command: ') ;
i := 1 ;
while not eoln do
begin
read(c) ;
buffer[i] := c ;
i := i + 1 ;
end ;
buffer[0] := chr(i-1) ;
readln ;
concatchar( buffer, ';' ) ;
while (not exit) and (length(buffer) > 1) do begin (* Get next command. *)
semipos := pos(';', buffer) ; (* Find boundary of next command. *)
copy( buffer, 1, semipos - 1,s ) ; (* Extract next command. *)
delete ( buffer, 1, semipos ) ; (* Remove next command from buffer. *)
if stack <> nil then oldtop := stack^.x ; (* Save old top of stack. *)
fpstatus.curexcep := [] ; (* Reset exception flags for new operation. *)
j := 0 ;
for i := 1 to length(s) do if s[i] <> ' ' then begin (* suppress blanks
and lower case *)
j := j + 1 ;
s[j] := upcase(s[i]) ;
end ;
copy(s,1,j,s) ;
for i := ord(s[0])+1 to maxfpstring do s[i] := ' ' ;
docommand ( found ) ;
if found and testflag then dotest ( s, found, x, y ) ;
end ;
until exit ;
end ;
procedure execute ;
begin
writeln(' Begin IEEE Calculator version 2 September 1985 ') ;
initialize ;
commandloop ;
end ;
#include 'dotest.i'
begin (* Outer block. *)
execute ;
end .
End-Of-File
echo Extracting calcf32.p
cat >calcf32.p <<'End-Of-File'
(* File calcf32.p, Version 9 October 1984. *)
(* This version of the calculator test unit tests 32 bit single precision
IEEE arithmetic accessed by the "shortreal" type. *)
#include 'sane.h'
#include 'oldfplib.h'
#include 'calctest.h'
#include 'calcsingle.h'
type bite = -128..+127 ;
function getbite ( b : bite ) : byt ;
begin
if b >= 0 then getbite := b else getbite := b + 256 ;
end ;
function setbite ( b : byt ) : bite ;
begin
if b < 128 then setbite := b else setbite := b - 256 ;
end ;
procedure swapexcep (* var e : excepset *) ;
var t : excepset ;
begin
e := [] ;
end ;
procedure swaptrap (* var e : excepset *) ;
var t : excepset ;
begin
e := [] ;
end ;
procedure swapmode (* var e : fpmodetype *) ;
var t : fpmodetype ;
begin
t.round := rnear ;
t.clos := affine ;
t.norm := normalizing ;
t.precision := xprec ;
e := t ;
end ;
procedure toreal ( s : csingle ; var r : shortreal ) ;
(* kluge to call a csingle a shortreal *)
type
klugetype = record
case boolean of
false : ( sk : packed array[0..3] of -128..+127 ) ;
true : ( rk : shortreal ) ;
end ;
var
kluge : klugetype ;
i : 0..3 ;
begin
for i := 0 to 3 do kluge.sk[i] := setbite(s[i]) ;
r := kluge.rk ;
end ;
procedure fromreal ( r : shortreal ; var s : csingle ) ;
(* kluge to call a shortreal a csingle *)
type
klugetype = record
case boolean of
false : ( sk : packed array[0..3] of -128..+127 ) ;
true : ( rk : shortreal ) ;
end ;
var
kluge : klugetype ;
i : 0..3 ;
begin
kluge.rk := r ;
for i := 0 to 3 do s[i] := getbite(kluge.sk[i]) ;
end ;
procedure pretest (* var storemode : arithtype *) ;
begin
storemode := flt32 ;
end ;
procedure tstore (* var z : internal *) ;
begin end ;
procedure tabs(* x : internal ; var z : internal *) ;
var
xs : csingle ; xr : shortreal ;
begin
tosingle(x,xs) ; toreal (xs,xr) ;
xr := abs(xr) ;
fromreal(xr,xs) ; unpacksingle(xs,z) ;
end ;
procedure tsqrt(* x : internal ; var z : internal *) ;
var
xs : csingle ; xr : shortreal ;
begin
tosingle(x,xs) ; toreal (xs,xr) ;
fromreal(sqrt(xr),xs) ; unpacksingle(xs,z) ;
end ;
procedure tneg(* x : internal ; var z : internal *) ;
var
xs : csingle ; xr : shortreal ;
begin
tosingle(x,xs) ; toreal (xs,xr) ;
xr := -xr ;
fromreal(xr,xs) ; unpacksingle(xs,z) ;
end ;
procedure tadd (* x, y : internal ; var z : internal *) ;
var
xs,ys,zs : csingle ;
xr,yr,zr : shortreal ;
begin
tosingle(x,xs) ; toreal(xs,xr) ; tosingle(y,ys) ; toreal(ys,yr) ;
zr := xr + yr ;
fromreal(zr,zs) ;
unpacksingle(zs,z) ;
end ;
procedure tsub (* x, y : internal ; var z : internal *) ;
var
xs,ys,zs : csingle ;
xr,yr,zr : shortreal ;
begin
tosingle(x,xs) ; toreal(xs,xr) ; tosingle(y,ys) ; toreal(ys,yr) ;
zr := xr - yr ;
fromreal(zr,zs) ;
unpacksingle(zs,z) ;
end ;
procedure tmul (* x, y : internal ; var z : internal *) ;
var
xs,ys,zs : csingle ;
xr,yr,zr : shortreal ;
begin
tosingle(x,xs) ; toreal(xs,xr) ; tosingle(y,ys) ; toreal(ys,yr) ;
zr := xr * yr ;
fromreal(zr,zs) ;
unpacksingle(zs,z) ;
end ;
procedure tdiv (* x, y : internal ; var z : internal *) ;
var
xs,ys,zs : csingle ;
xr,yr,zr : shortreal ;
begin
tosingle(x,xs) ; toreal(xs,xr) ; tosingle(y,ys) ; toreal(ys,yr) ;
zr := xr / yr ;
fromreal(zr,zs) ;
unpacksingle(zs,z) ;
end ;
procedure trem (* x, y : internal ; var z : internal *) ;
var
xs,ys,zs : csingle ;
xr,yr,zr : shortreal ;
begin
tosingle(x,xs) ; toreal(xs,xr) ; tosingle(y,ys) ; toreal(ys,yr) ;
fromreal(xr - yr * round(xr/yr),zs) ;
unpacksingle(zs,z) ;
end ;
procedure tcompare (* x, y : internal ; var cc : conditioncode *) ;
var
xs,ys,zs : csingle ;
xr,yr,zr : shortreal ;
begin
tosingle(x,xs) ; toreal(xs,xr) ; tosingle(y,ys) ; toreal(ys,yr) ;
write ( ' Tests affirm these predicates: ') ;
if xr=yr then write(' EQ ') ;
IF XR<>YR THEN write(' NE ') ;
IF XR<YR THEN write(' LT ') ;
IF XR<=YR THEN write(' LE ') ;
IF XR>YR THEN write(' GT ') ;
IF XR>=YR THEN write(' GE ') ;
writeln ;
IF xr=yr then cc := equal else
if xr>yr then cc := greater else
if xr<yr then cc := lesser else
cc := notord ;
end ;
procedure tconvert(* x : internal ; var z : internal ; a : arithtype *) ;
var yx : cextended ; yd : cdouble ; ys : csingle ;
yi64 : cint64 ; yi16 : integer ;
xs : csingle ; xr : shortreal ; xl : longint ;
begin
If a=i32 then begin
tosingle(x,xs) ; toreal(xs,xr) ;
xl := round(xr) ;
yi16 := xl ;
writeln(' Intermediate i16 ',yi16) ;
xr := xl ;
fromreal(xr,xs) ; unpacksingle(xs,z) ;
end
else begin
z := x ;
end
end ;
procedure tintconvert(* x : internal ; var z : internal ; a : arithtype *) ;
var yx : cextended ; yd : cdouble ; ys : csingle ;
yi64 : cint64 ; yi16 : integer ;
xs : csingle ; xr : shortreal ; xl : longint ;
begin
If a=i32 then begin
tosingle(x,xs) ; toreal(xs,xr) ;
xl := trunc(xr) ;
yi16 := xl ;
writeln(' Intermediate i16 ',yi16) ;
xr := xl ;
fromreal(xr,xs) ; unpacksingle(xs,z) ;
end
else begin
z := x ;
end
end ;
procedure tdisplay(* x : internal *) ;
var
xs : csingle ; xr : shortreal ;
s : fpstring ; i,j : integer ; error : boolean ;
begin
tosingle(x,xs) ; toreal(xs,xr) ;
{write (' Free ') ;
for i := 1 to 4 do begin
f32_ascii(xr,5*i- 1,0,0,fp_free,s,error ) ;
for j := length(s)+1 to 5*i-1 do write(' ') ;
write(' ',s) ;
end ;
writeln ;
write (' Lisa ') ;
for i := 1 to 4 do begin
f32_ascii(xr,5*i-1,0,0,fp_lisa,s,error ) ;
for j := length(s)+1 to 5*i-1 do write(' ') ;
write(' ',s) ;
end ;
writeln ;
}
writeln(' Efmt ',xr:5, xr:10, xr : 15, xr : 20 ) ;
writeln(' Ffmt ', xr : 5 : 0, xr : 10 : 5, xr : 15 : 7, xr : 20 : 10 ) ;
end ;
procedure tdecbin
(* s : fpstring ; var xout : internal ; var error : boolean *) ;
(* converts decimal fpstring s to internal format *)
(* error is set true if bad format *)
var
r : shortreal ;
xs : csingle ;
next : integer ;
f : text ;
i : integer ;
begin
rewrite(f) ;
for i := 1 to ord(s[0]) do write(f,s[i]) ;
writeln(f) ;
reset(f) ;
readln(f,r) ;
fromreal(r,xs) ; unpacksingle(xs,x) ;
end ;
End-Of-File
echo Extracting calcf64.p
cat >calcf64.p <<'End-Of-File'
(* File calcf64.p, Version 8 October 1984. *)
(* This version of the calculator test unit tests 64 bit double precision
IEEE arithmetic accessed by the "real" type. *)
#include 'sane.h'
#include 'oldfplib.h'
#include 'calctest.h'
#include 'calcdouble.h'
type bite = -128..+127 ;
function getbite ( b : bite ) : byt ;
begin
if b >= 0 then getbite := b else getbite := b + 256 ;
end ;
function setbite ( b : byt ) : bite ;
begin
if b < 128 then setbite := b else setbite := b - 256 ;
end ;
procedure swapexcep (* var e : excepset *) ;
var t : excepset ;
begin
e := [] ;
end ;
procedure swaptrap (* var e : excepset *) ;
var t : excepset ;
begin
e := [] ;
end ;
procedure swapmode (* var e : fpmodetype *) ;
var t : fpmodetype ;
begin
t.round := rnear ;
t.clos := affine ;
t.norm := normalizing ;
t.precision := xprec ;
e := t ;
end ;
procedure toreal ( s : cdouble ; var r : real ) ;
(* kluge to call a cdouble a real *)
type
klugetype = record
case boolean of
false : ( sk : packed array[0..7] of -128..+127 ) ;
true : ( rk : real ) ;
end ;
var
kluge : klugetype ;
i : 0..7 ;
begin
for i := 0 to 7 do kluge.sk[i] := setbite(s[i]) ;
r := kluge.rk ;
end ;
procedure fromreal ( r : real ; var s : cdouble ) ;
(* kluge to call a real a cdouble *)
type
klugetype = record
case boolean of
false : ( sk : packed array[0..7] of -128..+127 ) ;
true : ( rk : real ) ;
end ;
var
kluge : klugetype ;
i : 0..7 ;
begin
kluge.rk := r ;
for i := 0 to 7 do s[i] := getbite(kluge.sk[i]) ;
end ;
procedure pretest (* var storemode : arithtype *) ;
begin
storemode := f64 ;
end ;
procedure tstore (* var z : internal *) ;
begin end ;
procedure tabs(* x : internal ; var z : internal *) ;
var
xs : cdouble ; xr : real ;
begin
todouble(x,xs) ; toreal (xs,xr) ;
xr := abs(xr) ;
fromreal(xr,xs) ; unpackdouble(xs,z) ;
end ;
procedure tsqrt(* x : internal ; var z : internal *) ;
var
xs : cdouble ; xr : real ;
begin
todouble(x,xs) ; toreal (xs,xr) ;
fromreal(sqrt(xr),xs) ; unpackdouble(xs,z) ;
end ;
procedure tneg(* x : internal ; var z : internal *) ;
var
xs : cdouble ; xr : real ;
begin
todouble(x,xs) ; toreal (xs,xr) ;
xr := -xr ;
fromreal(xr,xs) ; unpackdouble(xs,z) ;
end ;
procedure tadd (* x, y : internal ; var z : internal *) ;
var
xs,ys,zs : cdouble ;
xr,yr,zr : real ;
begin
todouble(x,xs) ; toreal(xs,xr) ; todouble(y,ys) ; toreal(ys,yr) ;
zr := xr + yr ;
fromreal(zr,zs) ;
unpackdouble(zs,z) ;
end ;
procedure tsub (* x, y : internal ; var z : internal *) ;
var
xs,ys,zs : cdouble ;
xr,yr,zr : real ;
begin
todouble(x,xs) ; toreal(xs,xr) ; todouble(y,ys) ; toreal(ys,yr) ;
zr := xr - yr ;
fromreal(zr,zs) ;
unpackdouble(zs,z) ;
end ;
procedure tmul (* x, y : internal ; var z : internal *) ;
var
xs,ys,zs : cdouble ;
xr,yr,zr : real ;
begin
todouble(x,xs) ; toreal(xs,xr) ; todouble(y,ys) ; toreal(ys,yr) ;
zr := xr * yr ;
fromreal(zr,zs) ;
unpackdouble(zs,z) ;
end ;
procedure tdiv (* x, y : internal ; var z : internal *) ;
var
xs,ys,zs : cdouble ;
xr,yr,zr : real ;
begin
todouble(x,xs) ; toreal(xs,xr) ; todouble(y,ys) ; toreal(ys,yr) ;
zr := xr / yr ;
fromreal(zr,zs) ;
unpackdouble(zs,z) ;
end ;
procedure trem (* x, y : internal ; var z : internal *) ;
var
xs,ys,zs : cdouble ;
xr,yr,zr : real ;
begin
todouble(x,xs) ; toreal(xs,xr) ; todouble(y,ys) ; toreal(ys,yr) ;
fromreal(xr - yr * round(xr/yr),zs) ;
unpackdouble(zs,z) ;
end ;
procedure tcompare (* x, y : internal ; var cc : conditioncode *) ;
var
xs,ys,zs : cdouble ;
xr,yr,zr : real ;
begin
todouble(x,xs) ; toreal(xs,xr) ; todouble(y,ys) ; toreal(ys,yr) ;
write ( ' Tests affirm these predicates: ') ;
if xr=yr then write(' EQ ') ;
IF XR<>YR THEN write(' NE ') ;
IF XR<YR THEN write(' LT ') ;
IF XR<=YR THEN write(' LE ') ;
IF XR>YR THEN write(' GT ') ;
IF XR>=YR THEN write(' GE ') ;
writeln ;
IF xr=yr then cc := equal else
if xr>yr then cc := greater else
if xr<yr then cc := lesser else
cc := notord ;
end ;
procedure tconvert(* x : internal ; var z : internal ; a : arithtype *) ;
var yx : cextended ; yd : cdouble ; ys : csingle ;
yi64 : cint64 ; yi16 : integer ;
xs : cdouble ; xr : real ; xl : longint ;
begin
If a=i32 then begin
todouble(x,xs) ; toreal(xs,xr) ;
xl := round(xr) ;
yi16 := xl ;
writeln(' Intermediate i16 ',yi16) ;
xr := xl ;
fromreal(xr,xs) ; unpackdouble(xs,z) ;
end
else begin
z := x ;
end
end ;
procedure tintconvert(* x : internal ; var z : internal ; a : arithtype *) ;
var yx : cextended ; yd : cdouble ; ys : csingle ;
yi64 : cint64 ; yi16 : integer ;
xs : cdouble ; xr : real ; xl : longint ;
begin
If a=i32 then begin
todouble(x,xs) ; toreal(xs,xr) ;
xl := trunc(xr) ;
yi16 := xl ;
writeln(' Intermediate i16 ',yi16) ;
xr := xl ;
fromreal(xr,xs) ; unpackdouble(xs,z) ;
end
else begin
z := x ;
end
end ;
procedure tdisplay(* x : internal *) ;
var
xs : cdouble ; xr : real ;
s : fpstring ; i,j : integer ; error : boolean ;
begin
todouble(x,xs) ; toreal(xs,xr) ;
{write (' Free ') ;
for i := 1 to 4 do begin
f32_ascii(xr,5*i- 1,0,0,fp_free,s,error ) ;
for j := length(s)+1 to 5*i-1 do write(' ') ;
write(' ',s) ;
end ;
writeln ;
write (' Lisa ') ;
for i := 1 to 4 do begin
f32_ascii(xr,5*i-1,0,0,fp_lisa,s,error ) ;
for j := length(s)+1 to 5*i-1 do write(' ') ;
write(' ',s) ;
end ;
writeln ;
}
writeln(' Efmt ',xr:5, xr:10, xr : 15, xr : 20 ) ;
writeln(' Ffmt ', xr : 5 : 0, xr : 10 : 5, xr : 15 : 7, xr : 20 : 10 ) ;
end ;
procedure tdecbin
(* s : fpstring ; var xout : internal ; var error : boolean *) ;
(* converts decimal fpstring s to internal format *)
(* error is set true if bad format *)
var
r : real ;
xs : cdouble ;
next : integer ;
f : text ;
i : integer ;
begin
rewrite(f) ;
for i := 1 to ord(s[0]) do write(f,s[i]) ;
writeln(f) ;
reset(f) ;
readln(f,r) ;
fromreal(r,xs) ; unpackdouble(xs,x) ;
end ;
End-Of-File
echo Extracting calctest.p
cat >calctest.p <<'End-Of-File'
(* File calctest.p, Version 5 October 1984. *)
#include 'sane.h'
#include 'oldfplib.h'
#include 'calctest.h'
procedure pretest (* var storemode : arithtype *) ;
begin
end ;
procedure tstore (* var z : internal *) ;
begin
end ;
procedure tadd (* x, y : internal ; var z : internal *) ;
begin
end ;
procedure tsub (* x, y : internal ; var z : internal *) ;
begin
end ;
procedure tmul (* x, y : internal ; var z : internal *) ;
begin
end ;
procedure tdiv (* x, y : internal ; var z : internal *) ;
begin
end ;
procedure trem (* x, y : internal ; var z : internal *) ;
begin
end ;
procedure tcompare (* x, y : internal ; var cc : conditioncode *) ;
begin
end ;
procedure tconvert (* x : internal ; var z : internal ; a : arithtype *) ;
begin
end ;
procedure tintconvert
(* x : internal ; var z : internal ; a : arithtype *) ;
begin
end ;
procedure tabs (* x : internal ; var z : internal *) ;
begin
end ;
procedure tsqrt (* x : internal ; var z : internal *) ;
begin
end ;
procedure tneg (* x : internal ; var z : internal *) ;
begin
end ;
procedure tdisplay (* x : internal *) ;
begin
end ;
procedure tdecbin
(* s : fpstring ; var x : internal ; var error : boolean *) ;
begin
end ;
procedure swapexcep (* var e : excepset *) ;
begin
end ;
procedure swaptrap (* var e : excepset *) ;
begin
end ;
procedure swapmode (* var e : fpmodetype *) ;
begin
end ;
End-Of-File
echo ""
echo "End of Kit"
exitsources-request@panda.UUCP (09/04/85)
Mod.sources: Volume 3, Issue 5
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 calcdouble.h
cat >calcdouble.h <<'End-Of-File'
procedure todouble ( var x : internal ; var y : cdouble ) ; external ;
procedure unpackdouble ( y : cdouble ; var x : internal ) ; external ;
End-Of-File
echo Extracting calcsingle.h
cat >calcsingle.h <<'End-Of-File'
procedure tosingle ( var x : internal ; var y : csingle ) ; external ;
procedure unpacksingle ( y : csingle ; var x : internal ) ; external ;
End-Of-File
echo Extracting calctest.h
cat >calctest.h <<'End-Of-File'
(* File calctest.h, Version 5 October 1984. *)
(* This version of the calculator test unit is a dummy to provide
constant and type declarations only. *)
(* Global constant, type, and variable declarations for Calc. *)
const
stickybit = 66 ; (* position of sticky bit in internal representation *)
type
arithtype = ( i16, i32, i64, flt32, f64, ext80, unrounded ) ;
(* types of arithmetic operands *)
fpmodetype = record (* floating point mode record *)
round : rmode (* roundmodetype *) ;
precision : extprec ;
clos : closure(*type*) ;
norm : denorm ;
end ;
fpstype = record (* complete status of floating point unit *)
mode : fpmodetype ;
curexcep : excepset ; (* Set of exceptions generated by current op. *)
excep : excepset ;
trap : excepset end ;
internal = record (* internal extended format *)
(* unlike external extended, most significant bit represents 0.5,
not 1.0 *)
sign : boolean ;
exponent : integer ; (* range is -2**15 to 2**15-1 *)
significand : array [0..stickybit] of boolean ;
(* bit stickybit-2 is guard ;
bit (stickybit-1) is round ; bit stickybit is sticky *)
end ;
(* Following are temporary calculator internal types which use
logical bytes, which may not be the same as the physical bytes
specified in x80modes. *)
csingle = array [0..3] of byt ;
cdouble = array [0..7] of byt ;
cextended = array [0..9] of byt ;
cint64 = array [0..7] of byt ;
procedure pretest ( var storemode : arithtype ) ; external ;
procedure swapmode ( var e : fpmodetype ) ; external ;
procedure swaptrap ( var e : excepset ) ; external ;
procedure swapexcep ( var e : excepset ) ; external ;
procedure tneg ( x : internal ; var z : internal ) ; external ;
procedure tabs ( x : internal ; var z : internal ) ; external ;
procedure tsqrt ( x : internal ; var z : internal ) ; external ;
procedure tadd ( x, y : internal ; var z : internal ) ; external ;
procedure tsub ( x, y : internal ; var z : internal ) ; external ;
procedure tmul ( x, y : internal ; var z : internal ) ; external ;
procedure tdiv ( x, y : internal ; var z : internal ) ; external ;
procedure trem ( x, y : internal ; var z : internal ) ; external ;
procedure tcompare ( x, y : internal ; var cc : conditioncode ) ; external ;
procedure tstore ( storagemode : arithtype ; var z : internal ) ; external ;
procedure tconvert ( x : internal ; var z : internal ; a : arithtype ) ; external ;
procedure tintconvert ( x : internal ; var z : internal ; a : arithtype ) ; external ;
procedure tdisplay ( x : internal ) ; external ;
procedure tdecbin ( s : fpstring ; var x : internal ; var error : boolean ) ; external ;
procedure ffloat_ ; external ;
procedure ffunc_ ; external ;
procedure sfloat_ ; external ;
procedure sfunc_ ; external ;
End-Of-File
echo Extracting oldfplib.h
cat >oldfplib.h <<'End-Of-File'
const
maxfpstring = 80 ;
invop = invalid;
overfl = overflow;
underfl = underflow;
div0 = divbyzero;
inxact = inexact;
cvtovfl = invalid ;
type
byt = 0..255 ;
fpstring = packed array [0..maxfpstring] of char ;
roundtype = ( rnear, rzero, rpos, rneg, rout ) ;
rmode = rnear .. rneg ;
closure = (proj, affine) ;
denorm = ( warning, normalizing ) ;
extprec = ( xprec, sprec, dprec ) ;
xcpn = exception ;
excepset = set of exception ;
fp_cc = ( equal, lesser, greater, notord ) ;
conditioncode = fp_cc ;
End-Of-File
echo Extracting sane.h
cat >sane.h <<'End-Of-File'
type
longint = integer ;
integer = -32768..32767 ;
single = array [0..1] of integer ;
double = array [0..3] of integer ;
comp = array [0..3] of integer ;
extended = array [0..4] of integer ;
environ = integer ;
rounddir = ( tonearest, upward, downward, towardzero ) ;
relop = ( gt, lt, gl, eq, ge, le, gel, unord ) ;
exception = ( invalid, underflow, overflow, divbyzero, inexact ) ;
numclass = ( snan, qnan, infinite, zero, normal, denormal ) ;
roundprecision = ( extprecision, dblprecision, realprecision ) ;
procedure SetRnd ( r : rounddir ) ; external ;
function GetRnd : rounddir ; external ;
procedure SetXcp ( x : exception; onoff : boolean ); external ;
function TestXcp ( x : exception ) : boolean ; external ;
End-Of-File
echo ""
echo "End of Kit"
exitsources-request@panda.UUCP (09/04/85)
Mod.sources: Volume 3, Issue 6
Submitted by: decvax!decwrl!sun!dgh!dgh (David Hough)
#! /bin/sh
: make a directory, cd to it, and run this through sh
echo If this kit is complete, "End of Kit" will echo at the end
echo Extracting base.i
cat >base.i <<'End-Of-File'
(* File base.i, Version 8 October 1984. *)
procedure mpyten ( var x : internal ; n : integer ) ;
(* Multiplies x by 10**n, using table of powers of ten. *)
var
n1, n2 : integer ;
begin
n1 := abs(n) div 32 ;
n2 := abs(n) mod 32 ;
if n1 < 32 then begin
if n > 0 then begin
if n2 > 0 then multiply( x, tensmall[n2], x ) ;
if n1 > 0 then multiply( x, tenbig[n1], x) ;
end
else if n < 0 then begin
if n2 > 0 then divide ( x, tensmall[n2], x ) ;
if n1 > 0 then divide ( x, tenbig[n1], x ) ;
end ;
end
else begin (* n is too big. *)
if n > 0 then begin
makeinf(x) ;
setex ( overfl ) ;
end
else begin
makezero(x) ;
setex ( underfl ) ;
end ;
end ;
end ;
procedure xtodec ( x : internal ; r : roundtype ;
var s : strng ) ;
(* Converts abs(x) to an integral value, then that
value is converted to a strng of ASCII digits. *)
var
tf : excepset ;
acc : -20..+20 ; (* divide accumulator *)
i, j : integer ;
carry : boolean ;
nib : nibarray ;
last : integer ;
begin
roundint ( x, r, xprec ) ;
fpstatus.curexcep := fpstatus.curexcep - [inxact] ;
(* Don't care about inxact. *)
if kind(x) = zerokind then makeucsdstring('0 ',s) else begin
s[0] := chr(0) ;
last := x.exponent - 1 ; (* Last is the last active bit of the accumulator. *)
repeat (* Get one digit per cycle until there's nothing left. *)
(* For each digit, do a NON RESTORING divide by TEN. *)
acc := - 10 ;
for j := 0 to last do begin
(* Do one divide minicycle for each bit of dividend, plus one extra. *)
acc := acc + acc ; (* Double remainder. *)
if x.significand[j] then acc := acc + 1 ; (* Shift in next bit of dividend. *)
if acc < 0 then begin (* Do add-ten cycle. *)
acc := acc + 10 ;
end
else begin (* Do subtract-ten cycle. *)
acc := acc - 10 ;
end ;
x.significand[j] := acc >= 0 ; (* Sign of this step determines quotient bit
and whether add or subtract next time. *)
(* End around complement rotate for quotient bit. *)
end (* Divide mini-cycle. *) ;
if acc < 0 then begin
(* Remainder is negative so add 10. *)
acc := acc + 10 ;
end ;
precatchar(' ',s) ;
s[1] := chr(ord('0') + acc ) ;
j := firstbit ( x, 0, last ) ;
if j <= last then
begin
left(x, j) ;
end ;
last := last - j ;
until last < 0 ;
(* Keep doing divide cycles until the quotient is zero. *)
end ;
end ;
procedure subdec (* x : internal ; p1, p2 : integer ; var s: strng *) ;
(* s receives a strng of decimal digits representing the integer in
x.significand[p1]..x.significand[p2], right justified. *)
var
j, i : integer ;
nib : nibarray ;
begin
if p2 = stickybit then begin (* Avoid trying to donormalize sticky bit. *)
for i := p1 to p2 do
x.significand[i-1] := x.significand[i] ;
p1 := p1 - 1 ;
p2 := p2 - 1 ;
end ;
for i := 0 to (p1-1) do
x.significand[i] := false ; (* Clear bits outside field. *)
for i := (p2+1) to stickybit do
x.significand[i] := false ;
x.exponent := p2 + 1 ;
donormalize(x) ;
xtodec( x, rnear, s) ;
end ;
procedure findbinten ( x : internal ; n : integer ;
var s : strng ; var p : integer ) ;
(* Converts x into s and p, s a strng of exactly n significant
digits, so
x .=. s * 10^p *)
var
e1, e2, dp : integer ;
tf : excepset ;
t : internal ;
cc : conditioncode ;
i : integer ;
norm : boolean ;
et, dt : integer ;
fraction : integer ;
spill : boolean ;
begin
(* First ESTIMATE p.
We want a p such that
10^(n-1) <= int(abs(x)*10^p) < 10^n
For a first guess, use
n - log10( 2**e * 1+f )
which we approximate by
n - ((77+(1/16))/256) * (e + f )
e is broken into two pieces to get benefit of
22 bit product. *)
norm := x.significand[0] ; (* If x is not normalized then we don't force
n significant digits in E format. *)
e2 := (x.exponent-1) div 256 ; (* e = e1 + 256*e2 *)
e1 := (x.exponent-1) - 256 * e2 ;
fraction := xbyte(x,1,8) ;
e1 := 77 * e1 + 16 * e2 ; (* First order contribution from e1 and second
order from e2. *)
if norm then e1 := e1 + ( 77 * fraction ) div 256 ;
(* If normalized add in a contribution for the fraction.
If not normalized, assume significand of 1.00000... *)
et := e1 div 256 ; dt := e1 - et * 256 ;
(* We are never sure how Pascal div and mod work on negative
numbers. *)
if dt > 0 then et := et + 1 ;
(* but we try to get et as close as possible anyway to
the ceiling of e1/256. *)
e2 := 77 * e2 ;
p := n - (et + e2) ;
(* Now remedy flaws in approximation by altering p as necessary. *)
tf := fpstatus.curexcep ; (* Save exceptions. *)
dp := 0 ; (* Assume no correction required. *)
repeat
fpstatus.curexcep := tf ;
(* Restore exception status. Don't want inexact flag set
for inappropriate p. *)
p := p + dp ; (* Correct p. *)
dp := 0 ; (* Assume no correction required. *)
t := x ;
mpyten(t, p) ; (* Multiply t by appropriate power. *)
roundint( t, fpstatus.mode.round, xprec ) ;
t.sign := false ; (* Use absolute values for comparison. *)
compare ( t, tensmall[n], cc ) ; (* t must not exceed n sig digits. *)
case cc of
otherwise ;
greater : dp := -1 ; (* If too big, correct and repeat. *)
equal : begin
if norm or (n=1) then begin (* We want only n digits. *)
p := p - 1 ; (* If almost, avoid full repeat of process. *)
t := tensmall[n-1] ;
end
else begin (* If not normalized we want n-1 digits. *)
p := p - 2 ;
t := tensmall[n-2] ;
end ;
end ;
lesser : begin
compare(t,tensmall[n-1],cc) ;
if norm then begin (* If normalized, insure enough sig digits. *)
if cc = lesser then dp := + 1 ; (* If not enough digits, correct *)
end (* Need exactly n digits. *)
else begin (* If unnormalized, want no more than n-1 digits. *)
if cc <> lesser then dp := -1 ; (* Try again for less than n. *)
end ;
end ;
end ;
t.sign := x.sign ;
{
if dp <> 0 then writeln(' Power of ten correction: ',dp) ;
}
spill := ([underfl, overfl] * fpstatus.curexcep ) <> [] ;
until (dp = 0) or (kind(t)=zerokind) or spill ;
(* Repeat until no correction necessary or over/underfl occurs. *)
(* or the number rounds to normalized zero. *)
if spill then
begin
(* String of asterisks if overfl. *)
s[0] := chr(0) ;
while length(s) < n do concatchar(s,'*') ;
end
else
begin
xtodec(*2*) ( t, fpstatus.mode.round, s(*, n*) ) ; (* Get strng. *)
end ;
while length(s) < n do precatchar( '0', s ) ;
(* Add unnormalizing digits. *)
p := - p ;
end ;
procedure decint ( s : strng ; var x : internal ; var e : integer ) ;
(* DECINT converts a strng s of decimal digits into x and e
such that
s = x * 10^e
If s >= 2^(leastsigbit+1) then the sticky bit may be set. *)
const
last = 69 ; (* last bit of decint accumulator *)
var
i,j : integer ;
acc : array[0..last] of boolean ; (* Accumulator long enough to hold 20
significant digits and a sticky bit *)
carry, zero : boolean ;
n : nibarray ;
procedure tenmult ; (* multiplies acc by 10 *)
var
i : integer ;
carry : boolean ;
begin
for i := 0 to (last-3) do acc[i] := acc[i+3] ; (* multiply acc by 8 *)
for i := (last-2) to last do acc[i] := false ;
carry := false ;
for i := (last-1) downto 2 do
adder( acc[i], acc[i-2], acc[i], carry) ;
for i := 1 downto 0 do
adder( acc[i], false, acc[i], carry) ;
end ;
begin (* decint *)
for i := 0 to last do acc[i] := false ;
e := 0 ;
zero := true ;
for j := 1 to length(s) do begin
if s[j] = '0' then begin
if not zero then e := e + 1 ;
end
else begin
if not zero then begin
e := e + 1 ;
while (e>0) and (not acc[0]) and (not acc[1]) and (not acc[2]) and
(not acc[3]) do begin (* multiply by ten *)
tenmult ;
e := e - 1 ;
end ;
end ;
zero := false ;
if e > 0 then acc[last] := true (* set sticky bit on *)
else begin (* add digit *)
hexnibble( s[j], n) ;
carry := false ;
for i := (last-1) downto (last-4) do
adder( acc[i], n[i+4-last], acc[i], carry) ;
for i := (last-5) downto 0 do
adder( acc[i], false, acc[i], carry) ;
end ;
end ;
end ;
if zero then makezero(x) else begin
(* Now acc[0..(last-1)] represents an integer value,
acc[last] a sticky bit, to be multiplied by 10^e *)
i := 0 ;
while ( i < last ) and not acc[i] do i := i + 1 ;
(* search for first nonzero bit *)
x.exponent := last - i ; (* Set exponent of result *)
for j := 0 to (last-i-1) do acc[j] := acc[j+i] ;
(* normalize *)
for j := (last-i) to (last-1) do acc[j] := false ;
for i := 0 to (stickybit-1) do x.significand[i] := acc[i] ;
x.significand[stickybit] := acc[last] ;
for i := stickybit to (last-1) do
x.significand[stickybit] := x.significand[stickybit] or acc[i] ;
end ;
end ;
procedure putdec (* s : strng ; p1, p2 : integer ;
var x : internal ; var error : boolean *) ;
(* Interprets s as a decimal integer, puts value in bits
p1..p2 of x.significand.
Sets Error if any significant bits don't fit in field. *)
var
i, j : integer ;
y : internal ;
e : integer ;
f : excepset ;
begin
decint( s, y, e ) ;
f := fpstatus.curexcep ;
mpyten ( y, e ) ;
fpstatus.curexcep := f ; (* Ignore exceptions that may arise. *)
e := y.exponent ;
error := (e > (leastsigbit+1)) ;
(* Bad news if s too big to fit in 64 bits. *)
if not error then begin
if kind(y) = zerokind then begin
for i := p1 to p2 do
x.significand[i] := false ; (* Set up zero field. *)
end
else begin
if (p2-p1+1) >= e then begin (* y fits in field. *)
for i := p2 downto (p2+1-e) do
x.significand[i] := y.significand[i+e-1-p2] ;
for i := (p2-e) downto p1 do
x.significand[i] := false ;
end
else begin
for i := p2 downto p1 do
x.significand[i] := y.significand[i+e-1-p2] ;
for i := (p1-p2+e-2) downto 0 do
error := error or y.significand[i] ; (* Check for lost significant bits. *)
end end end
end ;
procedure todecint ( x : internal ; var s : strng ) ;
(* if x is an integer less than 2**15,
then s receives the decimal digits representing x.
Otherwise s is set to empty. *)
var
i, k : integer ;
begin
s[0] := chr(0) ;
if kind(x) = zerokind then makeucsdstring('0 ',s) else
if (abs(kind(x)) = normkind) and (x.exponent <= 15) and (x.exponent >= 1)
then begin
if zerofield ( x, x.exponent, stickybit ) then begin (* it's all integer *)
k := 0 ;
for i := 0 to (x.exponent-1) do begin (* Accumulate k. *)
k := k + k ;
if x.significand[i] then k := k + 1 ;
end ;
while k > 0 do begin
precatchar( chr(ord('0') + (k mod 10)), s) ;
k := k div 10 ;
end ;
if x.sign then precatchar( '-', s ) ;
end end
end ;
procedure bindec (* x : internal ; var s : strng *) ;
(* converts x to decimal format *)
var
e, i, j, k : integer ;
nib : nibarray ;
t : strng ;
tf : excepset ;
ns : integer ;
begin
case abs(kind(x)) of
zerokind : if x.sign
then makeucsdstring('-0',s) else makeucsdstring('0 ',s) ;
unnormkind, normkind : begin
todecint ( x, s ) ;
if length(s) < 1 then begin (* Can't represent as integer; too bad. *)
ns := 19 ;
case storagemode of (* Set number of significant digits output. *)
flt32 : ns := 9 ;
f64 : ns := 17 ;
otherwise ;
end ;
findbinten( x, ns, s, e ) ;
tf := fpstatus.curexcep ;
fpstatus.curexcep := [] ;
if not (overfl in tf) then begin
if e <> 0 then begin (* x not an integer so write it in E format. *)
e := e + ns-1 ;
for i := length(s) downto 2 do s[i+1] := s[i] ;
s[2] := '.' ;
s[0] := chr(length(s)+1) ;
if e <> 0 then begin
concatchar(s, 'E') ;
if e > 0 then concatchar(s, '+') ;
intdec(e, t) ;
s := concat(s,t) ;
end ;
end ;
end ;
if x.sign then precatchar('-',s) ;
end ;
end ;
infkind, nankind : nanascii ( x, false, s ) ;
otherwise
end ;
end ;
procedure decbin (* s : strng ; var x : internal ; var error : boolean *) ;
(* converts decimal strng s to internal format *)
(* error is set true if bad format *)
type
stringclass = (nonnumeric, truezero, nonzero) ; (* types of strng *)
var
class : stringclass ;
i, k, min : integer ;
e1, e2 : integer ;
sub : strng ;
t : strng ;
esign : boolean ;
nib : nibarray ;
ee, ie : integer ;
procedure checkadd ( x, y : integer ; var z : integer ; var error : boolean ) ;
(* Computes z := x + y except if z overflows error is set to true
and z is set to maxexp-1 or minexp+1 *)
begin
error := false ;
z := x + y ;
if (x>0) and (y>0) and (z<=0) then begin
z := maxexp - 1 ;
error := true ;
end
else if (x<0) and (y<0) and (z>=0) then begin
z := minexp + 1 ;
error := true ;
end ;
end ;
procedure bump ; (* removes first character from strng t *)
begin
delete (t,1,1)
end ;
begin
class := nonnumeric ;
error := false ;
esign := false ;
x.sign := false ;
x.exponent := 0 ;
ee := 0 ; ie := 0 ;
for i := 0 to stickybit do x.significand[i] := false ;
sub[0] := chr(0) ; (* substring for accumulating significant digits *)
t[0] := chr(0) ;
for i := 1 to length(s) do if s[i] <> ' ' then concatchar(t,upcase(s[i])) ;
concatchar(t,'!') ; (* this marks the end of the input strng *)
if t[1] = '+' then bump else if t[1] = '-' then begin (* handle negative *)
x.sign := true ;
bump
end ;
while t[1] = '0' do begin
class := truezero ;
bump ; (* delete leading zeros *)
end ;
while t[1] in digitset do begin (* digits before point *)
class := nonzero ;
concatchar(sub, t[1]) ;
bump
end ;
if t[1] = '.' then begin (* check for point *)
bump ;
while t[1] in digitset do begin (* process digits after point *)
if (t[1] <> '0') or (class = nonzero) then class := nonzero
else class := truezero ;
concatchar( sub, t[1]) ;
ie := ie - 1 ;
bump ;
end ;
end ;
ee := 0 ;
if t[1] = 'E' then bump ; (* handle E for exponent *)
if t[1] = '+' then bump else if t[1]='-' then begin (* exponent sign *)
esign := true ;
bump
end ;
while t[1] in digitset do begin (* exponent digits *)
if ee > ((maxexp - (ord(t[1])-ord('0'))) div 10 ) then begin
error := true ;
ee := maxexp - 1 ;
end else
begin
ee := 10 * ee + ord(t[1]) - ord('0') ;
end ; bump end ;
if class = truezero then x.exponent := minexp else begin
if esign then ee := -ee ;
checkadd(ee,ie,ee,error) ; (* ee := ee + ie *)
if not error then begin
(* Minimize ee if possible by adding zeros to sub *)
ie := 19 - length(sub) ; (* Maximum number of zeros to add. *)
if (ee>0) and (ie>0) then begin (* Go ahead and add. *)
if ee < ie then ie := ee ; (* Only add enough to reduce ee to 0. *)
ee := ee - ie ;
for i := 1 to ie do concatchar( sub, '0') ;
end ;
decint ( sub, x, ie ) ; (* Convert substring to x and ie *)
checkadd( ee, ie, ee, error ) ; (* Add in ie to exponent. *)
end ;
if not error then
mpyten ( x, ee ) ; (* Adjust x by appropriate power of ten. *)
end ;
if class = nonnumeric then
(* the following code checks for INFs and NANs *)
begin
NANer ( s, false, x, error )
end
else
if length(t) > 1 then error := true ;
if error then
begin
makenan(nanascbin,x) ;
end ;
end ;
procedure display (* x : internal *) ;
(* displays x in decimal and binary *)
begin
write(' Hex: ') ; displayhex(x) ;
write(' Dec: ') ; displaydec(x) ;
end ;
End-Of-File
echo Extracting utility.i
cat >utility.i <<'End-Of-File'
(* File utility.i, Version 8 October 1984. *)
function length ( var x : strng ) : integer ;
begin (* concat *)
length := ord(x[0]) ;
end (* concat *) ;
procedure displayhex ( x : internal ) ;
var s : strng ;
i : integer ;
begin
binhex(x,s) ;
for i := 1 to length(s) do write(s[i]) ;
writeln ;
end ;
procedure displaydec ( x : internal ) ;
var s : strng ;
i : integer ;
begin
bindec(x,s) ;
for i := 1 to length(s) do write(s[i]) ;
writeln ;
end ;
procedure concatchar ( var s : strng ; c : char ) ;
(* concatenates the character c onto the end of s *)
var
ls : integer ;
begin
ls := ord(s[0]) + 1 ;
s[ls] := c ;
s[0] := chr(ls) ;
end ;
function upcase ( c : char ) : char ;
begin
if ('a' <= c) and (c <= 'z') then upcase := chr(ord(c)-32) else upcase := c
end ;
function stackspace : integer ;
(* Computes number of stack entries left.
As this number approaches zero, operation becomes dangerous. *)
var space : integer ;
begin
stackspace := 10000 ;
end ;
procedure hexnibble ( h : char ; var n : nibarray ) ;
(* Converts ASCII hexit h into a nibarray *)
var
i, w : integer ;
begin
if h in digitset then w := ord(h)-ord('0') else w := ord(h) - ord('A') + 10 ;
for i := 3 downto 0 do begin
n[i] := odd(w) ;
w := w div 2 ;
end ;
end ;
function nibblehex (* n : nibarray ) : char *) ;
(* converts a nibarray into a hexit ASCII character *)
var
i, w : integer ;
c : char ;
begin
w := 0 ;
for i := 0 to 3 do begin
w := w + w ;
if n[i] then w := w + 1 ;
end ;
if w < 10 then c := chr(ord('0') + w) else c := chr(ord('A') + w - 10 ) ;
nibblehex := c ;
end ;
procedure displayexcep ( es : excepset ) ;
(* Displays exception names that are enabled. *)
var i : xcpn ;
begin
for i := invop to inexact do
if i in es then write(' ',xcpnname[i],' ') ;
end ;
procedure displaystatus ;
(* Displays current mode, trap, exception flags. *)
begin
write(' Modes: ') ;
case fpstatus.mode.round of
rneg : write(' RM ') ;
rpos : write(' RP ') ;
rzero : write(' RZ ') ;
otherwise
end ;
case fpstatus.mode.precision of
sprec: write(' R24 ') ;
dprec: write(' R53 ') ;
otherwise
end ;
if fpstatus.mode.clos = proj then write(' PROJ ') ;
if fpstatus.mode.norm = warning then write(' WARN ') ;
case storagemode of
i16 : write(' I16 ') ;
i32 : write(' I32 ') ;
i64 : write(' I64 ') ;
flt32 : write(' F32 ') ;
f64 : write(' F64 ') ;
ext80 : write(' X80 ') ;
otherwise
end ;
writeln ;
if fpstatus.trap <> [] then begin (* Write out trap line. *)
write(' Traps: ') ;
displayexcep( fpstatus.trap ) ;
writeln ;
end ;
if fpstatus.excep <> [] then begin (* Write out exceptions. *)
write(' Exceptions: ') ;
displayexcep( fpstatus.excep ) ;
writeln ;
end ;
end ;
procedure trapmessage ;
(* Announces name of trap that would occur now. *)
var
tset : excepset ;
f : xcpn ;
begin
tset := fpstatus.trap * fpstatus.curexcep ;
if tset <> [] then begin
f := invop ; (* Start with highest priority exception. *)
while (f <> cvtovfl) and not (f in tset) do f := succ(f) ;
if f in tset then writeln( xcpnname[f],' TRAP occurs. ') ;
end ;
end ;
procedure setex (* e : xcpn *) ;
(* Turns on exception flag in curexcep. *)
begin
fpstatus.curexcep := fpstatus.curexcep + [e] ;
end ;
function zerofield (* x : internal ; p1, p2 : integer ) : boolean *) ;
(* Returns true if x.significand[p1..p2] is all zeros. *)
var i : integer ;
begin
i := p1 ;
while ( i < p2 ) and not x.significand[i] do i := i + 1 ;
zerofield := ( i >= p2 ) and not x.significand[p2] ;
(* Can't test bit p2 in main loop ; would cause range error if
p2 were stickybit, on subsequent test. *)
end ;
function firstbit (* x : internal ; p1, p2 : integer ) : integer *) ;
(* Returns index of leftmost onebit in field
x.significand[p1..p2].
If field is zero, returns p2+1. *)
var i : integer ;
begin
i := p1 ;
while ( i < p2 ) and not x.significand[i] do i := i + 1 ;
if ( i >= p2 ) and not x.significand[p2] then i := p2+1 ;
(* Can't test bit p2 in main loop ; would cause range error if
p2 were stickybit, on subsequent test. *)
firstbit := i ;
end ;
function lastbit (* x : internal ; p1, p2 : integer ) : integer *) ;
(* Returns index of rightmost nonzero bit in field
x.significand[p1..p2].
If field is zero, returns p1-1. *)
var i : integer ;
begin
i := p2 ;
while ( i > p1 ) and not x.significand[i] do i := i - 1 ;
if ( i <= p1 ) and not x.significand[p1] then i := p1 - 1 ;
(* Can't test bit p1 in main loop ; would cause range error if
p1 were zero, on subsequent test. *)
lastbit := i ;
end ;
function kind (* x : internal ) : integer *) ;
(* returns kind(x) but all NANs have kind=4 in order to fit int16 *)
var
i, tkind : integer ;
begin
if x.exponent = maxexp then begin (* inf or nan *)
if zerofield ( x, 1, stickybit ) then tkind := ord(infkind)
else tkind := ord(nankind) ;
end
else
if (x.exponent <= minexp) and zerofield(x, 0, stickybit)
then tkind := ord(zerokind)
else if x.significand[0] = true then tkind := ord(normkind)
else tkind := ord(unnormkind) ;
if x.sign then tkind := -tkind ;
kind := tkind ;
end ;
procedure makezero (* var x : internal *) ;
(* makes x into a zero. Does not change sign of x *)
var i : integer ;
begin
x.exponent := minexp ;
for i := 0 to stickybit do x.significand[i] := false ;
end ;
procedure makeinf (* var x : internal *) ;
(* makes x into a infinity. Does not change sign of x *)
var i : integer ;
begin
x.exponent := maxexp ;
for i := 0 to stickybit do x.significand[i] := false ;
end ;
procedure makenan (* n : integer ; var x : internal *) ;
(* makes x a NAN and inserts n in its more significant field.
Sets NV Operand flag in curexcep. *)
var
i : integer ;
begin
x.exponent := maxexp ;
for i := 0 to stickybit do x.significand[i] := false ;
i := 15 ;
while n <> 0 do begin
x.significand[i] := odd(n) ;
n := n div 2 ;
i := i - 1 ;
end ;
setex( invop ) ;
end ;
function unzero (* x : internal ) : boolean *) ;
(* returns TRUE if x is an unnormalized zero,
FALSE otherwise *)
begin
unzero := zerofield( x, 0, stickybit ) and (x.exponent > minexp) ;
end ;
procedure pushstack (* x : internal *) ;
(* pushes x on stack *)
(* In case of NV exception and trapping NAN, makes the NAN
non-trapping. *)
var
p : pstack ;
begin
if stackspace >= 3 then begin
new(p) ;
if (invop in fpstatus.curexcep) and (abs(kind(x))=nankind) then
x.significand[1] := false ; (* By convention bit 1 determines trapping/
non-trapping. *)
p^.x := x ;
p^.next := stack ;
stack := p ;
end else
writeln(' ERROR: not enough space for push! ') ;
end ;
procedure popstack (* var x : internal *) ;
(* pops stack to x, or sets x to 0 if stack is empty *)
begin
if stack=nil then begin
x.sign := false ;
makezero(x) ;
end
else begin
x := stack^.x ;
stack := stack^.next ;
if (abs(kind(x))=nankind) and x.significand[1] then (* It's a Trapping NAN. *)
setex( invop ) ;
end
end ;
procedure donormalize (* var x : internal *) ;
(* normalizes x *)
(* Unnormalized zeros are set to normalized zeros.
a INFs and NANs are not changed *)
var
i, j : integer ;
begin
if x.exponent < maxexp then begin
i := firstbit( x, 0, stickybit ) ;
if i > stickybit then x.exponent := minexp (* zero *) else
if i > 0 then begin
x.exponent := x.exponent - i ;
for j := i to stickybit do
x.significand[j-i] := x.significand[j] ;
for j := ((stickybit+1)-i) to (stickybit-1) do
x.significand[j] := x.significand[stickybit] ;
{
if (x.significand[stickybit]) and (i>1) then
writeln(i,' ERROR: Normalizing Sticky Bit ') ;
(* It's OK to shift sticky bit into next to last position because
during rounding the last two positions are stuck together.
It is definitely not OK to shift a sticky bit any further left. *)
}
end ;
end ;
end ;
procedure right (* var x : internal ; n : integer *) ;
(* does sticky right shift of internal x *)
var
i : integer ;
begin
{
if (0 > n) or (n > stickybit) then
writeln(' Funny Right ',n) ;
}
if n > stickybit then n := stickybit ; (* It's all the same for large n. *)
x.significand[stickybit] := not zerofield( x, (stickybit-n), stickybit) ;
for i := (stickybit-1) downto n do
x.significand[i] := x.significand[i-n] ;
for i := (n-1) downto 0 do
x.significand[i] := false ;
end ;
procedure left (* var x : internal ; n : integer *) ;
(* Lefts shifts significand of x, n times *)
var i : integer ;
begin
{
if (0 > n) or ( n > stickybit) then
writeln(' Funny left ',n) ;
}
if n > stickybit then n := stickybit ; (* All the same for large n. *)
for i := 0 to (stickybit-1-n) do
x.significand[i] := x.significand[i+n] ;
{
if x.significand[stickybit] and (n>1) then
writeln(n,' Error: LEFT shift of STICKY bit ') ;
}
for i := (stickybit-n) to (stickybit-1) do
x.significand[i] := x.significand[stickybit] ;
end ;
procedure roundkcs (* var x : internal ; r : roundtype ;
p : precisetype *) ;
(* Rounds x according to rounding mode and rounding precision.
Sets Inexact flag in curexcep if appropriate. *)
var i : integer ;
akx : integer ;
procedure dorn ; (* round to nearest *)
var
i : integer ;
carry : boolean ;
begin
carry := true ;
i := (leastsigbit+1) ;
while (i>=0) and carry do begin
adder(x.significand[i], false, x.significand[i], carry) ;
i := i-1 ;
end ;
if carry then begin (* carry out of most significant bit occurred. *)
x.significand[0] := true ;
x.exponent := x.exponent + 1 ;
end
else begin (* Check for ambiguous case *)
if zerofield( x, leastsigbit+1, stickybit )
then (* It is the ambiguous case. *)
x.significand[leastsigbit] := false ; (* force round to even. *)
end ;
end ;
procedure doro ; (* Round away from zero (Outward Round) *)
var
i : integer ;
carry : boolean ;
begin
carry := false ;
for i := (leastsigbit+1) to stickybit do carry := carry or x.significand[i] ;
if carry then begin (* propagate a carry *)
i := leastsigbit ;
while (i >= 0) and carry do begin
adder(x.significand[i], false, x.significand[i], carry) ;
i := i - 1 ;
end ;
if carry then begin (* Carry out occurred, so renormalize. *)
x.significand[0] := true ;
x.exponent := x.exponent + 1 ;
end ;
end ;
end ;
begin (* round *)
akx := abs(kind(x)) ;
if akx in [unnormkind, normkind] then begin
case p of
sprec: begin
right(x, 40) ;
x.exponent := x.exponent + 40 ;
end ;
dprec : begin
right(x, 11) ;
x.exponent := x.exponent + 11 ;
end ;
otherwise
end ;
for i := (leastsigbit+1) to stickybit do if x.significand[i] then
begin
setex( inxact ) ;
end ;
case r of
rnear : begin
dorn ;
end ;
rneg : if x.sign then doro ;
rpos : if not x.sign then doro ;
otherwise
end ;
for i := (leastsigbit+1) to stickybit do x.significand[i] := false ;
(* Eliminate G, R, and S bits. *)
case p of
sprec: begin
left(x,39) ;
x.exponent := x.exponent - 39 ;
if not x.significand[0] then
begin
left( x, 1) ;
x.exponent := x.exponent - 1 ;
end ;
end ;
dprec: begin
left(x,10) ;
x.exponent := x.exponent - 10 ;
if not x.significand[0] then
begin
left(x,1) ;
x.exponent := x.exponent - 1 ;
end ;
end ;
otherwise
end ;
end ;
end ;
procedure roundint (* var x : internal ; r : roundtype ; p : precisetype *) ;
(* Rounds x to an integral value in accordance with modes. *)
var
akx, i, count : integer ;
begin
akx := abs(kind(x)) ;
if akx in [unnormkind, normkind] then begin
if (x.exponent >= (leastsigbit+1)) then count := 0
else if x.exponent <= (leastsigbit+1-stickybit) then count := stickybit
else count := (leastsigbit+1) - x.exponent ;
(* Compute shift count of bits to get rid of. *)
case p of (* But allow for rounding to shorter precisions, too. *)
sprec: if count < 40 then count := 40 ;
dprec: if count < 11 then count := 11 ;
otherwise
end ;
if count > 0 then right ( x, count) ;
roundkcs( x, r, xprec) ; (* Do rounding. *)
if count > leastsigbit then begin (* Limit left shifts
for 0 < x < 1 which must be rounded either to 0 or 1. *)
count := leastsigbit ;
x.exponent := 1 ;
end ;
if count > 0 then begin
left(x, count-1 ) ;
if x.significand[0] then x.exponent := x.exponent + 1 (* Rounding carry out. *)
else left(x, 1) ;
end ;
if zerofield ( x, 0, stickybit ) then x.exponent := minexp ;
(* No significant bits left so make it a true zero. *)
end end ;
procedure picknan (* x, y : internal ; var z : internal *) ;
(* Sets z to whichever of x or y is a NAN.
If both are NANs, sets z to the one with the
greatest significand. *)
var i : integer ;
begin
if abs(kind(x)) = nankind then
if abs(kind(y)) = nankind then begin
i := 0 ;
while (i <= leastsigbit) and (x.significand[i] = y.significand[i]) do
i := i + 1 ;
if x.significand[i] then z := x else z := y ;
end
else z := x else z := y
end ;
function equalinternal ( x1, x2 : internal ) : boolean ;
(* Returns true if x1 = x2 *)
var
t : boolean ;
i : integer ;
begin
t := (x1.sign=x2.sign) and (x1.exponent=x2.exponent) ;
if t then for i := 0 to stickybit do
t := t and (x1.significand[i]=x2.significand[i]) ;
equalinternal := t ;
end ;
function concat ( x, y : strng ) : strng ;
var t : strng ;
i, lx, ly : integer ;
begin (* concat *)
t := x ;
lx := ord(x[0]) ; ly := ord(y[0]) ;
for i := 1 to ly do t[lx+i] := y[i] ;
t[0] := chr(lx+ly) ;
concat := t ;
end (* concat *) ;
procedure precatchar ( c : char ; var x : strng ) ;
var i, ls : integer ;
begin (* precatchar *)
ls := ord(x[0]) ;
for i := ls downto 1 do x[i+1] := x[i] ;
x[1] := c ;
x[0] := chr(ls+1) ;
end (* precatchar *) ;
procedure delete ( var x : strng ; index, count : integer ) ;
var i : integer ;
begin (* delete *)
for i := index+count to length(x) do
x[i-count] := x[i] ;
x[0] := chr(length(x)-count) ;
end (* delete *) ;
procedure makeucsdstring( x : strng; var t : strng );
(* Converts a constant string to UCSD form. *)
var i, l : integer ;
begin
l := 0 ;
while (33 <= ord(x[l])) and (ord(x[l]) <= 126) do l := l + 1 ;
for i := l downto 1 do t[i] := x[i-1] ;
t[0] := chr(l) ;
end ;
procedure copy ( s : strng; index, count : integer ; var t : strng ) ;
var i : integer ; l : integer ;
begin (* copy *)
t[0] := chr(count) ;
for i := 1 to count do t[i] := s[index+i-1] ;
end (* copy *) ;
procedure insert( s : strng ; var d : strng ; index : integer ) ;
var i,ld,ls : integer ;
begin (* insert *)
ls := ord(s[0]) ; ld := ord(d[0]) ;
for i := ld downto index do d[i+ls] := d[i] ;
for i := ls downto 1 do d[index+i-1] := s[i] ;
d[0] := chr(ls+ld) ;
end (* insert *) ;
function pos ( c : char ; s : strng ) : integer ;
var i, l : integer ;
begin (* pos *)
l := ord(s[0]) ;
i := 1 ;
while (i <= l) and (s[i] <> c) do i := i + 1 ;
if i <= l then pos := i else pos := 0 ;
end (* pos *) ;
function sequal ( s, c : strng ) : boolean ;
(* Compares UCSD string s to C string c and returns true if equal. *)
var
i, ls, lc : integer ;
begin (* sequal *)
lc := 0;
ls := ord(s[0]) ;
while (33 <= ord(c[lc])) and (ord(c[lc]) <= 126) do lc := lc + 1 ;
if lc <> ls then
begin
sequal := false ;
end
else
begin
i := 1 ;
while (i <= ls) and (s[i] = c[i-1]) do i := i + 1 ;
sequal := i > ls ;
end ;
end (* sequal *) ;
End-Of-File
echo Extracting init.i
cat >init.i <<'End-Of-File'
(* File init.i, Version 4 February 1985. *)
PROCEDURE initialize ;
(* does all the initializing of tables and variables. *)
var error : boolean ;
procedure inittensmall ;
(* procedure to initialize the table of small powers of ten *)
var
i : integer ;
j : integer ;
x : internal ;
carry : boolean ;
last : integer ;
error : boolean ;
begin
(* Make 10^0=1 *)
for i := 1 to stickybit do x.significand[i] := false ;
x.sign := false ;
x.exponent := 1 ;
x.significand[0] := true ;
tensmall[0] := x ;
(* Make other exact powers of ten. *)
last := 0 ; (* Last non-zero bit. *)
for j := 1 to 28 do begin
x.exponent := x.exponent + 3 ; (* Multiply by 8 first. *)
last := last + 2 ; (* At least 2 more significant bits. *)
carry := false ;
for i := last downto 2 do
adder(x.significand[i-2], x.significand[i],x.significand[i],carry) ;
for i := 1 downto 0 do
adder(false, x.significand[i],x.significand[i],carry) ;
if carry then begin (* Overflowed slightly. *)
x.exponent := x.exponent + 1 ;
last := last + 1 ;
for i := last downto 1 do
x.significand[i] := x.significand[i-1] ;
x.significand[0] := true ;
end ;
tensmall[j] := x ;
end ;
hexbin(' .a18f 07d7 36b9 0be5 4 h 97', tensmall[29], error ) ;
hexbin(' .c9f2 c9cd 0467 4ede c h 100', tensmall[30], error ) ;
hexbin(' .fc6f 7c40 4581 2296 4 h 103', tensmall[31], error ) ;
end ;
procedure inittenbig ;
(* procedure to initalize the table of large powers of ten *)
var
error : boolean ;
begin
tenbig[0] := tensmall[0] ;
hexbin(' .9dc5 ada8 2b70 b59d c h 107',tenbig[1], error ) ;
hexbin(' .c278 1f49 ffcf a6d5 4 h 213',tenbig[2], error ) ;
hexbin(' .efb3 ab16 c59b 14a2 c h 319',tenbig[3], error ) ;
hexbin(' .93ba 47c9 80e9 8cdf c h 426',tenbig[4], error ) ;
hexbin(' .b616 a12b 7fe6 17aa 4 h 532',tenbig[5], error ) ;
hexbin(' .e070 f78d 3927 556a c h 638',tenbig[6], error ) ;
hexbin(' .8a52 96ff e33c c92f c h 745',tenbig[7], error ) ;
hexbin(' .aa7e ebfb 9df9 de8d c h 851',tenbig[8], error ) ;
hexbin(' .d226 fc19 5c6a 2f8c 4 h 957',tenbig[9], error ) ;
hexbin(' .8184 2f29 f2cc e375 c h 1064',tenbig[10], error ) ;
hexbin(' .9fa4 2700 db90 0ad2 4 h 1170',tenbig[11], error ) ;
hexbin(' .c4c5 e310 aef8 aa17 4 h 1276',tenbig[12], error ) ;
hexbin(' .f28a 9c07 e9b0 9c58 c h 1382',tenbig[13], error ) ;
hexbin(' .957a 4ae1 ebf7 f3d3 c h 1489',tenbig[14],error) ;
hexbin(' .b83e d8dc 0795 a262 4 h 1595',tenbig[15], error) ;
end ;
procedure inittenhuge ;
var
error : boolean ;
begin
hexbin(' .e319 a0ae a60e 91c6 c h 1701',tenbig[16], error) ;
hexbin(' .8bf6 1451 432d 7bc2 c h 1808', tenbig[17], error) ;
hexbin(' .ac83 fb89 6b67 95fc c h 1914', tenbig[18], error) ;
hexbin(' .d4a4 4fb4 b8fa 79af c h 2020', tenbig[19], error) ;
hexbin(' .830c f791 e54a 9d1c c h 2127', tenbig[20], error) ;
hexbin(' .a188 4b69 ade2 4964 4 h 2233', tenbig[21], error) ;
hexbin(' .c71a a36a 1f8f 01cb c h 2339', tenbig[22], error) ;
hexbin(' .f56a 298f 4370 28f3 4 h 2445', tenbig[23], error) ;
hexbin(' .973f 9ca8 cd00 a68c 4 h 2552', tenbig[24], error) ;
hexbin(' .ba6d 9b40 d7cc 9ecc c h 2658', tenbig[25], error) ;
hexbin(' .e5ca 5a0b 8d73 7f0e 4 h 2764', tenbig[26], error) ;
hexbin(' .8d9e 89d1 1346 bda5 4 h 2871', tenbig[27], error) ;
hexbin(' .ae8f 2b2c e3d5 dbe9 c h 2977', tenbig[28], error) ;
hexbin(' .d729 3020 5a0c 1b2f c h 3083', tenbig[29], error) ;
hexbin(' .849a 672a 0d2e cfd1 c h 3190', tenbig[30], error) ;
hexbin(' .a372 2c13 41fa 93de 4 h 3296', tenbig[31], error) ;
end ;
begin
digitset := [ '0' .. '9' ] ;
hexset := digitset + [ 'A' .. 'F' ] ;
stack := nil ;
storagemode := unrounded ;
testflag := false ;
fpstatus.mode.round := rnear ;
fpstatus.mode.precision := xprec ;
fpstatus.mode.clos := affine ;
fpstatus.mode.norm := normalizing ;
fpstatus.curexcep := [] ;
fpstatus.excep := [] ;
fpstatus.trap := [] ;
leftnan[1] := 0 ;
leftnan[2] := 24 ;
leftnan[3] := 53 ;
leftnan[4] := leastsigbit + 1 ;
rightnan[1] := leftnan[2] - 1 ;
rightnan[2] := leftnan[3] - 1 ;
rightnan[3] := leftnan[4] - 1 ;
rightnan[4] := stickybit ;
xcpnname[cvtovfl] := 'IV' ;
xcpnname[overfl] := 'OV' ;
xcpnname[underfl] := 'UN' ;
xcpnname[div0] := 'D0' ;
xcpnname[invop] := 'NO' ;
xcpnname[inxact] := 'NX' ;
inittensmall ;
inittenbig ;
inittenhuge ;
(*
decbin ( ' 3.1415926535 89793 23846 26433', pi, error ) ;
decbin ( ' 2.7182818284 59045 23536 02874', e, error ) ;
*)
hexbin ( ' .c90fdaa22168c234c h 2', pi, error ) ;
decbin ( ' 2.7182818284 59045 23536 02874', e, error ) ;
end ;
End-Of-File
echo ""
echo "End of Kit"
exitsources-request@panda.UUCP (09/04/85)
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"
exitsources-request@panda.UUCP (09/04/85)
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