IEEE Calculator (part 4 of 6)
sources-request at panda.UUCP
sources-request at panda.UUCP
Wed Sep 4 12:11:44 AEST 1985
Mod.sources: Volume 3, Issue 6
Submitted by: decvax!decwrl!sun!dgh!dgh (David Hough)
#! /bin/sh
: make a directory, cd to it, and run this through sh
echo If this kit is complete, "End of Kit" will echo at the end
echo Extracting base.i
cat >base.i <<'End-Of-File'
(* File base.i, Version 8 October 1984. *)
procedure mpyten ( var x : internal ; n : integer ) ;
(* Multiplies x by 10**n, using table of powers of ten. *)
var
n1, n2 : integer ;
begin
n1 := abs(n) div 32 ;
n2 := abs(n) mod 32 ;
if n1 < 32 then begin
if n > 0 then begin
if n2 > 0 then multiply( x, tensmall[n2], x ) ;
if n1 > 0 then multiply( x, tenbig[n1], x) ;
end
else if n < 0 then begin
if n2 > 0 then divide ( x, tensmall[n2], x ) ;
if n1 > 0 then divide ( x, tenbig[n1], x ) ;
end ;
end
else begin (* n is too big. *)
if n > 0 then begin
makeinf(x) ;
setex ( overfl ) ;
end
else begin
makezero(x) ;
setex ( underfl ) ;
end ;
end ;
end ;
procedure xtodec ( x : internal ; r : roundtype ;
var s : strng ) ;
(* Converts abs(x) to an integral value, then that
value is converted to a strng of ASCII digits. *)
var
tf : excepset ;
acc : -20..+20 ; (* divide accumulator *)
i, j : integer ;
carry : boolean ;
nib : nibarray ;
last : integer ;
begin
roundint ( x, r, xprec ) ;
fpstatus.curexcep := fpstatus.curexcep - [inxact] ;
(* Don't care about inxact. *)
if kind(x) = zerokind then makeucsdstring('0 ',s) else begin
s[0] := chr(0) ;
last := x.exponent - 1 ; (* Last is the last active bit of the accumulator. *)
repeat (* Get one digit per cycle until there's nothing left. *)
(* For each digit, do a NON RESTORING divide by TEN. *)
acc := - 10 ;
for j := 0 to last do begin
(* Do one divide minicycle for each bit of dividend, plus one extra. *)
acc := acc + acc ; (* Double remainder. *)
if x.significand[j] then acc := acc + 1 ; (* Shift in next bit of dividend. *)
if acc < 0 then begin (* Do add-ten cycle. *)
acc := acc + 10 ;
end
else begin (* Do subtract-ten cycle. *)
acc := acc - 10 ;
end ;
x.significand[j] := acc >= 0 ; (* Sign of this step determines quotient bit
and whether add or subtract next time. *)
(* End around complement rotate for quotient bit. *)
end (* Divide mini-cycle. *) ;
if acc < 0 then begin
(* Remainder is negative so add 10. *)
acc := acc + 10 ;
end ;
precatchar(' ',s) ;
s[1] := chr(ord('0') + acc ) ;
j := firstbit ( x, 0, last ) ;
if j <= last then
begin
left(x, j) ;
end ;
last := last - j ;
until last < 0 ;
(* Keep doing divide cycles until the quotient is zero. *)
end ;
end ;
procedure subdec (* x : internal ; p1, p2 : integer ; var s: strng *) ;
(* s receives a strng of decimal digits representing the integer in
x.significand[p1]..x.significand[p2], right justified. *)
var
j, i : integer ;
nib : nibarray ;
begin
if p2 = stickybit then begin (* Avoid trying to donormalize sticky bit. *)
for i := p1 to p2 do
x.significand[i-1] := x.significand[i] ;
p1 := p1 - 1 ;
p2 := p2 - 1 ;
end ;
for i := 0 to (p1-1) do
x.significand[i] := false ; (* Clear bits outside field. *)
for i := (p2+1) to stickybit do
x.significand[i] := false ;
x.exponent := p2 + 1 ;
donormalize(x) ;
xtodec( x, rnear, s) ;
end ;
procedure findbinten ( x : internal ; n : integer ;
var s : strng ; var p : integer ) ;
(* Converts x into s and p, s a strng of exactly n significant
digits, so
x .=. s * 10^p *)
var
e1, e2, dp : integer ;
tf : excepset ;
t : internal ;
cc : conditioncode ;
i : integer ;
norm : boolean ;
et, dt : integer ;
fraction : integer ;
spill : boolean ;
begin
(* First ESTIMATE p.
We want a p such that
10^(n-1) <= int(abs(x)*10^p) < 10^n
For a first guess, use
n - log10( 2**e * 1+f )
which we approximate by
n - ((77+(1/16))/256) * (e + f )
e is broken into two pieces to get benefit of
22 bit product. *)
norm := x.significand[0] ; (* If x is not normalized then we don't force
n significant digits in E format. *)
e2 := (x.exponent-1) div 256 ; (* e = e1 + 256*e2 *)
e1 := (x.exponent-1) - 256 * e2 ;
fraction := xbyte(x,1,8) ;
e1 := 77 * e1 + 16 * e2 ; (* First order contribution from e1 and second
order from e2. *)
if norm then e1 := e1 + ( 77 * fraction ) div 256 ;
(* If normalized add in a contribution for the fraction.
If not normalized, assume significand of 1.00000... *)
et := e1 div 256 ; dt := e1 - et * 256 ;
(* We are never sure how Pascal div and mod work on negative
numbers. *)
if dt > 0 then et := et + 1 ;
(* but we try to get et as close as possible anyway to
the ceiling of e1/256. *)
e2 := 77 * e2 ;
p := n - (et + e2) ;
(* Now remedy flaws in approximation by altering p as necessary. *)
tf := fpstatus.curexcep ; (* Save exceptions. *)
dp := 0 ; (* Assume no correction required. *)
repeat
fpstatus.curexcep := tf ;
(* Restore exception status. Don't want inexact flag set
for inappropriate p. *)
p := p + dp ; (* Correct p. *)
dp := 0 ; (* Assume no correction required. *)
t := x ;
mpyten(t, p) ; (* Multiply t by appropriate power. *)
roundint( t, fpstatus.mode.round, xprec ) ;
t.sign := false ; (* Use absolute values for comparison. *)
compare ( t, tensmall[n], cc ) ; (* t must not exceed n sig digits. *)
case cc of
otherwise ;
greater : dp := -1 ; (* If too big, correct and repeat. *)
equal : begin
if norm or (n=1) then begin (* We want only n digits. *)
p := p - 1 ; (* If almost, avoid full repeat of process. *)
t := tensmall[n-1] ;
end
else begin (* If not normalized we want n-1 digits. *)
p := p - 2 ;
t := tensmall[n-2] ;
end ;
end ;
lesser : begin
compare(t,tensmall[n-1],cc) ;
if norm then begin (* If normalized, insure enough sig digits. *)
if cc = lesser then dp := + 1 ; (* If not enough digits, correct *)
end (* Need exactly n digits. *)
else begin (* If unnormalized, want no more than n-1 digits. *)
if cc <> lesser then dp := -1 ; (* Try again for less than n. *)
end ;
end ;
end ;
t.sign := x.sign ;
{
if dp <> 0 then writeln(' Power of ten correction: ',dp) ;
}
spill := ([underfl, overfl] * fpstatus.curexcep ) <> [] ;
until (dp = 0) or (kind(t)=zerokind) or spill ;
(* Repeat until no correction necessary or over/underfl occurs. *)
(* or the number rounds to normalized zero. *)
if spill then
begin
(* String of asterisks if overfl. *)
s[0] := chr(0) ;
while length(s) < n do concatchar(s,'*') ;
end
else
begin
xtodec(*2*) ( t, fpstatus.mode.round, s(*, n*) ) ; (* Get strng. *)
end ;
while length(s) < n do precatchar( '0', s ) ;
(* Add unnormalizing digits. *)
p := - p ;
end ;
procedure decint ( s : strng ; var x : internal ; var e : integer ) ;
(* DECINT converts a strng s of decimal digits into x and e
such that
s = x * 10^e
If s >= 2^(leastsigbit+1) then the sticky bit may be set. *)
const
last = 69 ; (* last bit of decint accumulator *)
var
i,j : integer ;
acc : array[0..last] of boolean ; (* Accumulator long enough to hold 20
significant digits and a sticky bit *)
carry, zero : boolean ;
n : nibarray ;
procedure tenmult ; (* multiplies acc by 10 *)
var
i : integer ;
carry : boolean ;
begin
for i := 0 to (last-3) do acc[i] := acc[i+3] ; (* multiply acc by 8 *)
for i := (last-2) to last do acc[i] := false ;
carry := false ;
for i := (last-1) downto 2 do
adder( acc[i], acc[i-2], acc[i], carry) ;
for i := 1 downto 0 do
adder( acc[i], false, acc[i], carry) ;
end ;
begin (* decint *)
for i := 0 to last do acc[i] := false ;
e := 0 ;
zero := true ;
for j := 1 to length(s) do begin
if s[j] = '0' then begin
if not zero then e := e + 1 ;
end
else begin
if not zero then begin
e := e + 1 ;
while (e>0) and (not acc[0]) and (not acc[1]) and (not acc[2]) and
(not acc[3]) do begin (* multiply by ten *)
tenmult ;
e := e - 1 ;
end ;
end ;
zero := false ;
if e > 0 then acc[last] := true (* set sticky bit on *)
else begin (* add digit *)
hexnibble( s[j], n) ;
carry := false ;
for i := (last-1) downto (last-4) do
adder( acc[i], n[i+4-last], acc[i], carry) ;
for i := (last-5) downto 0 do
adder( acc[i], false, acc[i], carry) ;
end ;
end ;
end ;
if zero then makezero(x) else begin
(* Now acc[0..(last-1)] represents an integer value,
acc[last] a sticky bit, to be multiplied by 10^e *)
i := 0 ;
while ( i < last ) and not acc[i] do i := i + 1 ;
(* search for first nonzero bit *)
x.exponent := last - i ; (* Set exponent of result *)
for j := 0 to (last-i-1) do acc[j] := acc[j+i] ;
(* normalize *)
for j := (last-i) to (last-1) do acc[j] := false ;
for i := 0 to (stickybit-1) do x.significand[i] := acc[i] ;
x.significand[stickybit] := acc[last] ;
for i := stickybit to (last-1) do
x.significand[stickybit] := x.significand[stickybit] or acc[i] ;
end ;
end ;
procedure putdec (* s : strng ; p1, p2 : integer ;
var x : internal ; var error : boolean *) ;
(* Interprets s as a decimal integer, puts value in bits
p1..p2 of x.significand.
Sets Error if any significant bits don't fit in field. *)
var
i, j : integer ;
y : internal ;
e : integer ;
f : excepset ;
begin
decint( s, y, e ) ;
f := fpstatus.curexcep ;
mpyten ( y, e ) ;
fpstatus.curexcep := f ; (* Ignore exceptions that may arise. *)
e := y.exponent ;
error := (e > (leastsigbit+1)) ;
(* Bad news if s too big to fit in 64 bits. *)
if not error then begin
if kind(y) = zerokind then begin
for i := p1 to p2 do
x.significand[i] := false ; (* Set up zero field. *)
end
else begin
if (p2-p1+1) >= e then begin (* y fits in field. *)
for i := p2 downto (p2+1-e) do
x.significand[i] := y.significand[i+e-1-p2] ;
for i := (p2-e) downto p1 do
x.significand[i] := false ;
end
else begin
for i := p2 downto p1 do
x.significand[i] := y.significand[i+e-1-p2] ;
for i := (p1-p2+e-2) downto 0 do
error := error or y.significand[i] ; (* Check for lost significant bits. *)
end end end
end ;
procedure todecint ( x : internal ; var s : strng ) ;
(* if x is an integer less than 2**15,
then s receives the decimal digits representing x.
Otherwise s is set to empty. *)
var
i, k : integer ;
begin
s[0] := chr(0) ;
if kind(x) = zerokind then makeucsdstring('0 ',s) else
if (abs(kind(x)) = normkind) and (x.exponent <= 15) and (x.exponent >= 1)
then begin
if zerofield ( x, x.exponent, stickybit ) then begin (* it's all integer *)
k := 0 ;
for i := 0 to (x.exponent-1) do begin (* Accumulate k. *)
k := k + k ;
if x.significand[i] then k := k + 1 ;
end ;
while k > 0 do begin
precatchar( chr(ord('0') + (k mod 10)), s) ;
k := k div 10 ;
end ;
if x.sign then precatchar( '-', s ) ;
end end
end ;
procedure bindec (* x : internal ; var s : strng *) ;
(* converts x to decimal format *)
var
e, i, j, k : integer ;
nib : nibarray ;
t : strng ;
tf : excepset ;
ns : integer ;
begin
case abs(kind(x)) of
zerokind : if x.sign
then makeucsdstring('-0',s) else makeucsdstring('0 ',s) ;
unnormkind, normkind : begin
todecint ( x, s ) ;
if length(s) < 1 then begin (* Can't represent as integer; too bad. *)
ns := 19 ;
case storagemode of (* Set number of significant digits output. *)
flt32 : ns := 9 ;
f64 : ns := 17 ;
otherwise ;
end ;
findbinten( x, ns, s, e ) ;
tf := fpstatus.curexcep ;
fpstatus.curexcep := [] ;
if not (overfl in tf) then begin
if e <> 0 then begin (* x not an integer so write it in E format. *)
e := e + ns-1 ;
for i := length(s) downto 2 do s[i+1] := s[i] ;
s[2] := '.' ;
s[0] := chr(length(s)+1) ;
if e <> 0 then begin
concatchar(s, 'E') ;
if e > 0 then concatchar(s, '+') ;
intdec(e, t) ;
s := concat(s,t) ;
end ;
end ;
end ;
if x.sign then precatchar('-',s) ;
end ;
end ;
infkind, nankind : nanascii ( x, false, s ) ;
otherwise
end ;
end ;
procedure decbin (* s : strng ; var x : internal ; var error : boolean *) ;
(* converts decimal strng s to internal format *)
(* error is set true if bad format *)
type
stringclass = (nonnumeric, truezero, nonzero) ; (* types of strng *)
var
class : stringclass ;
i, k, min : integer ;
e1, e2 : integer ;
sub : strng ;
t : strng ;
esign : boolean ;
nib : nibarray ;
ee, ie : integer ;
procedure checkadd ( x, y : integer ; var z : integer ; var error : boolean ) ;
(* Computes z := x + y except if z overflows error is set to true
and z is set to maxexp-1 or minexp+1 *)
begin
error := false ;
z := x + y ;
if (x>0) and (y>0) and (z<=0) then begin
z := maxexp - 1 ;
error := true ;
end
else if (x<0) and (y<0) and (z>=0) then begin
z := minexp + 1 ;
error := true ;
end ;
end ;
procedure bump ; (* removes first character from strng t *)
begin
delete (t,1,1)
end ;
begin
class := nonnumeric ;
error := false ;
esign := false ;
x.sign := false ;
x.exponent := 0 ;
ee := 0 ; ie := 0 ;
for i := 0 to stickybit do x.significand[i] := false ;
sub[0] := chr(0) ; (* substring for accumulating significant digits *)
t[0] := chr(0) ;
for i := 1 to length(s) do if s[i] <> ' ' then concatchar(t,upcase(s[i])) ;
concatchar(t,'!') ; (* this marks the end of the input strng *)
if t[1] = '+' then bump else if t[1] = '-' then begin (* handle negative *)
x.sign := true ;
bump
end ;
while t[1] = '0' do begin
class := truezero ;
bump ; (* delete leading zeros *)
end ;
while t[1] in digitset do begin (* digits before point *)
class := nonzero ;
concatchar(sub, t[1]) ;
bump
end ;
if t[1] = '.' then begin (* check for point *)
bump ;
while t[1] in digitset do begin (* process digits after point *)
if (t[1] <> '0') or (class = nonzero) then class := nonzero
else class := truezero ;
concatchar( sub, t[1]) ;
ie := ie - 1 ;
bump ;
end ;
end ;
ee := 0 ;
if t[1] = 'E' then bump ; (* handle E for exponent *)
if t[1] = '+' then bump else if t[1]='-' then begin (* exponent sign *)
esign := true ;
bump
end ;
while t[1] in digitset do begin (* exponent digits *)
if ee > ((maxexp - (ord(t[1])-ord('0'))) div 10 ) then begin
error := true ;
ee := maxexp - 1 ;
end else
begin
ee := 10 * ee + ord(t[1]) - ord('0') ;
end ; bump end ;
if class = truezero then x.exponent := minexp else begin
if esign then ee := -ee ;
checkadd(ee,ie,ee,error) ; (* ee := ee + ie *)
if not error then begin
(* Minimize ee if possible by adding zeros to sub *)
ie := 19 - length(sub) ; (* Maximum number of zeros to add. *)
if (ee>0) and (ie>0) then begin (* Go ahead and add. *)
if ee < ie then ie := ee ; (* Only add enough to reduce ee to 0. *)
ee := ee - ie ;
for i := 1 to ie do concatchar( sub, '0') ;
end ;
decint ( sub, x, ie ) ; (* Convert substring to x and ie *)
checkadd( ee, ie, ee, error ) ; (* Add in ie to exponent. *)
end ;
if not error then
mpyten ( x, ee ) ; (* Adjust x by appropriate power of ten. *)
end ;
if class = nonnumeric then
(* the following code checks for INFs and NANs *)
begin
NANer ( s, false, x, error )
end
else
if length(t) > 1 then error := true ;
if error then
begin
makenan(nanascbin,x) ;
end ;
end ;
procedure display (* x : internal *) ;
(* displays x in decimal and binary *)
begin
write(' Hex: ') ; displayhex(x) ;
write(' Dec: ') ; displaydec(x) ;
end ;
End-Of-File
echo Extracting utility.i
cat >utility.i <<'End-Of-File'
(* File utility.i, Version 8 October 1984. *)
function length ( var x : strng ) : integer ;
begin (* concat *)
length := ord(x[0]) ;
end (* concat *) ;
procedure displayhex ( x : internal ) ;
var s : strng ;
i : integer ;
begin
binhex(x,s) ;
for i := 1 to length(s) do write(s[i]) ;
writeln ;
end ;
procedure displaydec ( x : internal ) ;
var s : strng ;
i : integer ;
begin
bindec(x,s) ;
for i := 1 to length(s) do write(s[i]) ;
writeln ;
end ;
procedure concatchar ( var s : strng ; c : char ) ;
(* concatenates the character c onto the end of s *)
var
ls : integer ;
begin
ls := ord(s[0]) + 1 ;
s[ls] := c ;
s[0] := chr(ls) ;
end ;
function upcase ( c : char ) : char ;
begin
if ('a' <= c) and (c <= 'z') then upcase := chr(ord(c)-32) else upcase := c
end ;
function stackspace : integer ;
(* Computes number of stack entries left.
As this number approaches zero, operation becomes dangerous. *)
var space : integer ;
begin
stackspace := 10000 ;
end ;
procedure hexnibble ( h : char ; var n : nibarray ) ;
(* Converts ASCII hexit h into a nibarray *)
var
i, w : integer ;
begin
if h in digitset then w := ord(h)-ord('0') else w := ord(h) - ord('A') + 10 ;
for i := 3 downto 0 do begin
n[i] := odd(w) ;
w := w div 2 ;
end ;
end ;
function nibblehex (* n : nibarray ) : char *) ;
(* converts a nibarray into a hexit ASCII character *)
var
i, w : integer ;
c : char ;
begin
w := 0 ;
for i := 0 to 3 do begin
w := w + w ;
if n[i] then w := w + 1 ;
end ;
if w < 10 then c := chr(ord('0') + w) else c := chr(ord('A') + w - 10 ) ;
nibblehex := c ;
end ;
procedure displayexcep ( es : excepset ) ;
(* Displays exception names that are enabled. *)
var i : xcpn ;
begin
for i := invop to inexact do
if i in es then write(' ',xcpnname[i],' ') ;
end ;
procedure displaystatus ;
(* Displays current mode, trap, exception flags. *)
begin
write(' Modes: ') ;
case fpstatus.mode.round of
rneg : write(' RM ') ;
rpos : write(' RP ') ;
rzero : write(' RZ ') ;
otherwise
end ;
case fpstatus.mode.precision of
sprec: write(' R24 ') ;
dprec: write(' R53 ') ;
otherwise
end ;
if fpstatus.mode.clos = proj then write(' PROJ ') ;
if fpstatus.mode.norm = warning then write(' WARN ') ;
case storagemode of
i16 : write(' I16 ') ;
i32 : write(' I32 ') ;
i64 : write(' I64 ') ;
flt32 : write(' F32 ') ;
f64 : write(' F64 ') ;
ext80 : write(' X80 ') ;
otherwise
end ;
writeln ;
if fpstatus.trap <> [] then begin (* Write out trap line. *)
write(' Traps: ') ;
displayexcep( fpstatus.trap ) ;
writeln ;
end ;
if fpstatus.excep <> [] then begin (* Write out exceptions. *)
write(' Exceptions: ') ;
displayexcep( fpstatus.excep ) ;
writeln ;
end ;
end ;
procedure trapmessage ;
(* Announces name of trap that would occur now. *)
var
tset : excepset ;
f : xcpn ;
begin
tset := fpstatus.trap * fpstatus.curexcep ;
if tset <> [] then begin
f := invop ; (* Start with highest priority exception. *)
while (f <> cvtovfl) and not (f in tset) do f := succ(f) ;
if f in tset then writeln( xcpnname[f],' TRAP occurs. ') ;
end ;
end ;
procedure setex (* e : xcpn *) ;
(* Turns on exception flag in curexcep. *)
begin
fpstatus.curexcep := fpstatus.curexcep + [e] ;
end ;
function zerofield (* x : internal ; p1, p2 : integer ) : boolean *) ;
(* Returns true if x.significand[p1..p2] is all zeros. *)
var i : integer ;
begin
i := p1 ;
while ( i < p2 ) and not x.significand[i] do i := i + 1 ;
zerofield := ( i >= p2 ) and not x.significand[p2] ;
(* Can't test bit p2 in main loop ; would cause range error if
p2 were stickybit, on subsequent test. *)
end ;
function firstbit (* x : internal ; p1, p2 : integer ) : integer *) ;
(* Returns index of leftmost onebit in field
x.significand[p1..p2].
If field is zero, returns p2+1. *)
var i : integer ;
begin
i := p1 ;
while ( i < p2 ) and not x.significand[i] do i := i + 1 ;
if ( i >= p2 ) and not x.significand[p2] then i := p2+1 ;
(* Can't test bit p2 in main loop ; would cause range error if
p2 were stickybit, on subsequent test. *)
firstbit := i ;
end ;
function lastbit (* x : internal ; p1, p2 : integer ) : integer *) ;
(* Returns index of rightmost nonzero bit in field
x.significand[p1..p2].
If field is zero, returns p1-1. *)
var i : integer ;
begin
i := p2 ;
while ( i > p1 ) and not x.significand[i] do i := i - 1 ;
if ( i <= p1 ) and not x.significand[p1] then i := p1 - 1 ;
(* Can't test bit p1 in main loop ; would cause range error if
p1 were zero, on subsequent test. *)
lastbit := i ;
end ;
function kind (* x : internal ) : integer *) ;
(* returns kind(x) but all NANs have kind=4 in order to fit int16 *)
var
i, tkind : integer ;
begin
if x.exponent = maxexp then begin (* inf or nan *)
if zerofield ( x, 1, stickybit ) then tkind := ord(infkind)
else tkind := ord(nankind) ;
end
else
if (x.exponent <= minexp) and zerofield(x, 0, stickybit)
then tkind := ord(zerokind)
else if x.significand[0] = true then tkind := ord(normkind)
else tkind := ord(unnormkind) ;
if x.sign then tkind := -tkind ;
kind := tkind ;
end ;
procedure makezero (* var x : internal *) ;
(* makes x into a zero. Does not change sign of x *)
var i : integer ;
begin
x.exponent := minexp ;
for i := 0 to stickybit do x.significand[i] := false ;
end ;
procedure makeinf (* var x : internal *) ;
(* makes x into a infinity. Does not change sign of x *)
var i : integer ;
begin
x.exponent := maxexp ;
for i := 0 to stickybit do x.significand[i] := false ;
end ;
procedure makenan (* n : integer ; var x : internal *) ;
(* makes x a NAN and inserts n in its more significant field.
Sets NV Operand flag in curexcep. *)
var
i : integer ;
begin
x.exponent := maxexp ;
for i := 0 to stickybit do x.significand[i] := false ;
i := 15 ;
while n <> 0 do begin
x.significand[i] := odd(n) ;
n := n div 2 ;
i := i - 1 ;
end ;
setex( invop ) ;
end ;
function unzero (* x : internal ) : boolean *) ;
(* returns TRUE if x is an unnormalized zero,
FALSE otherwise *)
begin
unzero := zerofield( x, 0, stickybit ) and (x.exponent > minexp) ;
end ;
procedure pushstack (* x : internal *) ;
(* pushes x on stack *)
(* In case of NV exception and trapping NAN, makes the NAN
non-trapping. *)
var
p : pstack ;
begin
if stackspace >= 3 then begin
new(p) ;
if (invop in fpstatus.curexcep) and (abs(kind(x))=nankind) then
x.significand[1] := false ; (* By convention bit 1 determines trapping/
non-trapping. *)
p^.x := x ;
p^.next := stack ;
stack := p ;
end else
writeln(' ERROR: not enough space for push! ') ;
end ;
procedure popstack (* var x : internal *) ;
(* pops stack to x, or sets x to 0 if stack is empty *)
begin
if stack=nil then begin
x.sign := false ;
makezero(x) ;
end
else begin
x := stack^.x ;
stack := stack^.next ;
if (abs(kind(x))=nankind) and x.significand[1] then (* It's a Trapping NAN. *)
setex( invop ) ;
end
end ;
procedure donormalize (* var x : internal *) ;
(* normalizes x *)
(* Unnormalized zeros are set to normalized zeros.
a INFs and NANs are not changed *)
var
i, j : integer ;
begin
if x.exponent < maxexp then begin
i := firstbit( x, 0, stickybit ) ;
if i > stickybit then x.exponent := minexp (* zero *) else
if i > 0 then begin
x.exponent := x.exponent - i ;
for j := i to stickybit do
x.significand[j-i] := x.significand[j] ;
for j := ((stickybit+1)-i) to (stickybit-1) do
x.significand[j] := x.significand[stickybit] ;
{
if (x.significand[stickybit]) and (i>1) then
writeln(i,' ERROR: Normalizing Sticky Bit ') ;
(* It's OK to shift sticky bit into next to last position because
during rounding the last two positions are stuck together.
It is definitely not OK to shift a sticky bit any further left. *)
}
end ;
end ;
end ;
procedure right (* var x : internal ; n : integer *) ;
(* does sticky right shift of internal x *)
var
i : integer ;
begin
{
if (0 > n) or (n > stickybit) then
writeln(' Funny Right ',n) ;
}
if n > stickybit then n := stickybit ; (* It's all the same for large n. *)
x.significand[stickybit] := not zerofield( x, (stickybit-n), stickybit) ;
for i := (stickybit-1) downto n do
x.significand[i] := x.significand[i-n] ;
for i := (n-1) downto 0 do
x.significand[i] := false ;
end ;
procedure left (* var x : internal ; n : integer *) ;
(* Lefts shifts significand of x, n times *)
var i : integer ;
begin
{
if (0 > n) or ( n > stickybit) then
writeln(' Funny left ',n) ;
}
if n > stickybit then n := stickybit ; (* All the same for large n. *)
for i := 0 to (stickybit-1-n) do
x.significand[i] := x.significand[i+n] ;
{
if x.significand[stickybit] and (n>1) then
writeln(n,' Error: LEFT shift of STICKY bit ') ;
}
for i := (stickybit-n) to (stickybit-1) do
x.significand[i] := x.significand[stickybit] ;
end ;
procedure roundkcs (* var x : internal ; r : roundtype ;
p : precisetype *) ;
(* Rounds x according to rounding mode and rounding precision.
Sets Inexact flag in curexcep if appropriate. *)
var i : integer ;
akx : integer ;
procedure dorn ; (* round to nearest *)
var
i : integer ;
carry : boolean ;
begin
carry := true ;
i := (leastsigbit+1) ;
while (i>=0) and carry do begin
adder(x.significand[i], false, x.significand[i], carry) ;
i := i-1 ;
end ;
if carry then begin (* carry out of most significant bit occurred. *)
x.significand[0] := true ;
x.exponent := x.exponent + 1 ;
end
else begin (* Check for ambiguous case *)
if zerofield( x, leastsigbit+1, stickybit )
then (* It is the ambiguous case. *)
x.significand[leastsigbit] := false ; (* force round to even. *)
end ;
end ;
procedure doro ; (* Round away from zero (Outward Round) *)
var
i : integer ;
carry : boolean ;
begin
carry := false ;
for i := (leastsigbit+1) to stickybit do carry := carry or x.significand[i] ;
if carry then begin (* propagate a carry *)
i := leastsigbit ;
while (i >= 0) and carry do begin
adder(x.significand[i], false, x.significand[i], carry) ;
i := i - 1 ;
end ;
if carry then begin (* Carry out occurred, so renormalize. *)
x.significand[0] := true ;
x.exponent := x.exponent + 1 ;
end ;
end ;
end ;
begin (* round *)
akx := abs(kind(x)) ;
if akx in [unnormkind, normkind] then begin
case p of
sprec: begin
right(x, 40) ;
x.exponent := x.exponent + 40 ;
end ;
dprec : begin
right(x, 11) ;
x.exponent := x.exponent + 11 ;
end ;
otherwise
end ;
for i := (leastsigbit+1) to stickybit do if x.significand[i] then
begin
setex( inxact ) ;
end ;
case r of
rnear : begin
dorn ;
end ;
rneg : if x.sign then doro ;
rpos : if not x.sign then doro ;
otherwise
end ;
for i := (leastsigbit+1) to stickybit do x.significand[i] := false ;
(* Eliminate G, R, and S bits. *)
case p of
sprec: begin
left(x,39) ;
x.exponent := x.exponent - 39 ;
if not x.significand[0] then
begin
left( x, 1) ;
x.exponent := x.exponent - 1 ;
end ;
end ;
dprec: begin
left(x,10) ;
x.exponent := x.exponent - 10 ;
if not x.significand[0] then
begin
left(x,1) ;
x.exponent := x.exponent - 1 ;
end ;
end ;
otherwise
end ;
end ;
end ;
procedure roundint (* var x : internal ; r : roundtype ; p : precisetype *) ;
(* Rounds x to an integral value in accordance with modes. *)
var
akx, i, count : integer ;
begin
akx := abs(kind(x)) ;
if akx in [unnormkind, normkind] then begin
if (x.exponent >= (leastsigbit+1)) then count := 0
else if x.exponent <= (leastsigbit+1-stickybit) then count := stickybit
else count := (leastsigbit+1) - x.exponent ;
(* Compute shift count of bits to get rid of. *)
case p of (* But allow for rounding to shorter precisions, too. *)
sprec: if count < 40 then count := 40 ;
dprec: if count < 11 then count := 11 ;
otherwise
end ;
if count > 0 then right ( x, count) ;
roundkcs( x, r, xprec) ; (* Do rounding. *)
if count > leastsigbit then begin (* Limit left shifts
for 0 < x < 1 which must be rounded either to 0 or 1. *)
count := leastsigbit ;
x.exponent := 1 ;
end ;
if count > 0 then begin
left(x, count-1 ) ;
if x.significand[0] then x.exponent := x.exponent + 1 (* Rounding carry out. *)
else left(x, 1) ;
end ;
if zerofield ( x, 0, stickybit ) then x.exponent := minexp ;
(* No significant bits left so make it a true zero. *)
end end ;
procedure picknan (* x, y : internal ; var z : internal *) ;
(* Sets z to whichever of x or y is a NAN.
If both are NANs, sets z to the one with the
greatest significand. *)
var i : integer ;
begin
if abs(kind(x)) = nankind then
if abs(kind(y)) = nankind then begin
i := 0 ;
while (i <= leastsigbit) and (x.significand[i] = y.significand[i]) do
i := i + 1 ;
if x.significand[i] then z := x else z := y ;
end
else z := x else z := y
end ;
function equalinternal ( x1, x2 : internal ) : boolean ;
(* Returns true if x1 = x2 *)
var
t : boolean ;
i : integer ;
begin
t := (x1.sign=x2.sign) and (x1.exponent=x2.exponent) ;
if t then for i := 0 to stickybit do
t := t and (x1.significand[i]=x2.significand[i]) ;
equalinternal := t ;
end ;
function concat ( x, y : strng ) : strng ;
var t : strng ;
i, lx, ly : integer ;
begin (* concat *)
t := x ;
lx := ord(x[0]) ; ly := ord(y[0]) ;
for i := 1 to ly do t[lx+i] := y[i] ;
t[0] := chr(lx+ly) ;
concat := t ;
end (* concat *) ;
procedure precatchar ( c : char ; var x : strng ) ;
var i, ls : integer ;
begin (* precatchar *)
ls := ord(x[0]) ;
for i := ls downto 1 do x[i+1] := x[i] ;
x[1] := c ;
x[0] := chr(ls+1) ;
end (* precatchar *) ;
procedure delete ( var x : strng ; index, count : integer ) ;
var i : integer ;
begin (* delete *)
for i := index+count to length(x) do
x[i-count] := x[i] ;
x[0] := chr(length(x)-count) ;
end (* delete *) ;
procedure makeucsdstring( x : strng; var t : strng );
(* Converts a constant string to UCSD form. *)
var i, l : integer ;
begin
l := 0 ;
while (33 <= ord(x[l])) and (ord(x[l]) <= 126) do l := l + 1 ;
for i := l downto 1 do t[i] := x[i-1] ;
t[0] := chr(l) ;
end ;
procedure copy ( s : strng; index, count : integer ; var t : strng ) ;
var i : integer ; l : integer ;
begin (* copy *)
t[0] := chr(count) ;
for i := 1 to count do t[i] := s[index+i-1] ;
end (* copy *) ;
procedure insert( s : strng ; var d : strng ; index : integer ) ;
var i,ld,ls : integer ;
begin (* insert *)
ls := ord(s[0]) ; ld := ord(d[0]) ;
for i := ld downto index do d[i+ls] := d[i] ;
for i := ls downto 1 do d[index+i-1] := s[i] ;
d[0] := chr(ls+ld) ;
end (* insert *) ;
function pos ( c : char ; s : strng ) : integer ;
var i, l : integer ;
begin (* pos *)
l := ord(s[0]) ;
i := 1 ;
while (i <= l) and (s[i] <> c) do i := i + 1 ;
if i <= l then pos := i else pos := 0 ;
end (* pos *) ;
function sequal ( s, c : strng ) : boolean ;
(* Compares UCSD string s to C string c and returns true if equal. *)
var
i, ls, lc : integer ;
begin (* sequal *)
lc := 0;
ls := ord(s[0]) ;
while (33 <= ord(c[lc])) and (ord(c[lc]) <= 126) do lc := lc + 1 ;
if lc <> ls then
begin
sequal := false ;
end
else
begin
i := 1 ;
while (i <= ls) and (s[i] = c[i-1]) do i := i + 1 ;
sequal := i > ls ;
end ;
end (* sequal *) ;
End-Of-File
echo Extracting init.i
cat >init.i <<'End-Of-File'
(* File init.i, Version 4 February 1985. *)
PROCEDURE initialize ;
(* does all the initializing of tables and variables. *)
var error : boolean ;
procedure inittensmall ;
(* procedure to initialize the table of small powers of ten *)
var
i : integer ;
j : integer ;
x : internal ;
carry : boolean ;
last : integer ;
error : boolean ;
begin
(* Make 10^0=1 *)
for i := 1 to stickybit do x.significand[i] := false ;
x.sign := false ;
x.exponent := 1 ;
x.significand[0] := true ;
tensmall[0] := x ;
(* Make other exact powers of ten. *)
last := 0 ; (* Last non-zero bit. *)
for j := 1 to 28 do begin
x.exponent := x.exponent + 3 ; (* Multiply by 8 first. *)
last := last + 2 ; (* At least 2 more significant bits. *)
carry := false ;
for i := last downto 2 do
adder(x.significand[i-2], x.significand[i],x.significand[i],carry) ;
for i := 1 downto 0 do
adder(false, x.significand[i],x.significand[i],carry) ;
if carry then begin (* Overflowed slightly. *)
x.exponent := x.exponent + 1 ;
last := last + 1 ;
for i := last downto 1 do
x.significand[i] := x.significand[i-1] ;
x.significand[0] := true ;
end ;
tensmall[j] := x ;
end ;
hexbin(' .a18f 07d7 36b9 0be5 4 h 97', tensmall[29], error ) ;
hexbin(' .c9f2 c9cd 0467 4ede c h 100', tensmall[30], error ) ;
hexbin(' .fc6f 7c40 4581 2296 4 h 103', tensmall[31], error ) ;
end ;
procedure inittenbig ;
(* procedure to initalize the table of large powers of ten *)
var
error : boolean ;
begin
tenbig[0] := tensmall[0] ;
hexbin(' .9dc5 ada8 2b70 b59d c h 107',tenbig[1], error ) ;
hexbin(' .c278 1f49 ffcf a6d5 4 h 213',tenbig[2], error ) ;
hexbin(' .efb3 ab16 c59b 14a2 c h 319',tenbig[3], error ) ;
hexbin(' .93ba 47c9 80e9 8cdf c h 426',tenbig[4], error ) ;
hexbin(' .b616 a12b 7fe6 17aa 4 h 532',tenbig[5], error ) ;
hexbin(' .e070 f78d 3927 556a c h 638',tenbig[6], error ) ;
hexbin(' .8a52 96ff e33c c92f c h 745',tenbig[7], error ) ;
hexbin(' .aa7e ebfb 9df9 de8d c h 851',tenbig[8], error ) ;
hexbin(' .d226 fc19 5c6a 2f8c 4 h 957',tenbig[9], error ) ;
hexbin(' .8184 2f29 f2cc e375 c h 1064',tenbig[10], error ) ;
hexbin(' .9fa4 2700 db90 0ad2 4 h 1170',tenbig[11], error ) ;
hexbin(' .c4c5 e310 aef8 aa17 4 h 1276',tenbig[12], error ) ;
hexbin(' .f28a 9c07 e9b0 9c58 c h 1382',tenbig[13], error ) ;
hexbin(' .957a 4ae1 ebf7 f3d3 c h 1489',tenbig[14],error) ;
hexbin(' .b83e d8dc 0795 a262 4 h 1595',tenbig[15], error) ;
end ;
procedure inittenhuge ;
var
error : boolean ;
begin
hexbin(' .e319 a0ae a60e 91c6 c h 1701',tenbig[16], error) ;
hexbin(' .8bf6 1451 432d 7bc2 c h 1808', tenbig[17], error) ;
hexbin(' .ac83 fb89 6b67 95fc c h 1914', tenbig[18], error) ;
hexbin(' .d4a4 4fb4 b8fa 79af c h 2020', tenbig[19], error) ;
hexbin(' .830c f791 e54a 9d1c c h 2127', tenbig[20], error) ;
hexbin(' .a188 4b69 ade2 4964 4 h 2233', tenbig[21], error) ;
hexbin(' .c71a a36a 1f8f 01cb c h 2339', tenbig[22], error) ;
hexbin(' .f56a 298f 4370 28f3 4 h 2445', tenbig[23], error) ;
hexbin(' .973f 9ca8 cd00 a68c 4 h 2552', tenbig[24], error) ;
hexbin(' .ba6d 9b40 d7cc 9ecc c h 2658', tenbig[25], error) ;
hexbin(' .e5ca 5a0b 8d73 7f0e 4 h 2764', tenbig[26], error) ;
hexbin(' .8d9e 89d1 1346 bda5 4 h 2871', tenbig[27], error) ;
hexbin(' .ae8f 2b2c e3d5 dbe9 c h 2977', tenbig[28], error) ;
hexbin(' .d729 3020 5a0c 1b2f c h 3083', tenbig[29], error) ;
hexbin(' .849a 672a 0d2e cfd1 c h 3190', tenbig[30], error) ;
hexbin(' .a372 2c13 41fa 93de 4 h 3296', tenbig[31], error) ;
end ;
begin
digitset := [ '0' .. '9' ] ;
hexset := digitset + [ 'A' .. 'F' ] ;
stack := nil ;
storagemode := unrounded ;
testflag := false ;
fpstatus.mode.round := rnear ;
fpstatus.mode.precision := xprec ;
fpstatus.mode.clos := affine ;
fpstatus.mode.norm := normalizing ;
fpstatus.curexcep := [] ;
fpstatus.excep := [] ;
fpstatus.trap := [] ;
leftnan[1] := 0 ;
leftnan[2] := 24 ;
leftnan[3] := 53 ;
leftnan[4] := leastsigbit + 1 ;
rightnan[1] := leftnan[2] - 1 ;
rightnan[2] := leftnan[3] - 1 ;
rightnan[3] := leftnan[4] - 1 ;
rightnan[4] := stickybit ;
xcpnname[cvtovfl] := 'IV' ;
xcpnname[overfl] := 'OV' ;
xcpnname[underfl] := 'UN' ;
xcpnname[div0] := 'D0' ;
xcpnname[invop] := 'NO' ;
xcpnname[inxact] := 'NX' ;
inittensmall ;
inittenbig ;
inittenhuge ;
(*
decbin ( ' 3.1415926535 89793 23846 26433', pi, error ) ;
decbin ( ' 2.7182818284 59045 23536 02874', e, error ) ;
*)
hexbin ( ' .c90fdaa22168c234c h 2', pi, error ) ;
decbin ( ' 2.7182818284 59045 23536 02874', e, error ) ;
end ;
End-Of-File
echo ""
echo "End of Kit"
exit
More information about the Mod.sources
mailing list