-- FloatFns.mesa
-- Copywrite Xerox Corporation 1980
DIRECTORY
RealDefs: FROM "RealDefs",
Mopcodes: FROM "Mopcodes" USING [zPOP],
FloatFnDefs: FROM "FloatFnDefs";
FloatFns: PROGRAM IMPORTS RealDefs EXPORTS FloatFnDefs =
BEGIN
FloatNum: TYPE = RECORD [m2: CARDINAL,sign:[0..1],exp:[0..400B),m1:[0..200B)];
Precision: CARDINAL ← 5;
LoopCount: INTEGER ← 3;
magic: CARDINAL ← 200B-22;
LowHalf: PROCEDURE [l:LONG INTEGER] RETURNS [i:INTEGER]=MACHINE CODE
BEGIN Mopcodes.zPOP; END;
--compute e↑n using continued fractions
--See ACM Algorithm 229, Morelock, James C., "Elementary functions
--by continued fractions"
Exp: PUBLIC PROCEDURE [x: REAL] RETURNS [REAL] =
BEGIN
fl: FloatNum;
Ln2,LogBase2ofE,r,f,kludge: REAL;
k: LONG INTEGER;
i,ikludge: INTEGER;
one:REAL←1;
IF x = 0 THEN RETURN[one];
fl ← [sign: 0,exp: 200B,m1: 61B,m2: 71027B];
Ln2←LOOPHOLE[fl,REAL];
fl ← [sign: 0,exp: 201B,m1: 70B,m2: 125074B];
LogBase2ofE←LOOPHOLE[fl,REAL];
x←x*LogBase2ofE;--exponent of 2 (instead of e)
k←RealDefs.Fix[x];--integer part
x←x-k;--fractional part
x←x*Ln2;--now result = 2↑k * exp(x) and 0<=x<Ln2
r←x*x;
--f←4*LoopCount+2;
i←4*LoopCount+2;f←i;
FOR i DECREASING IN [1..LoopCount] DO
ikludge←4*i-2;
kludge←ikludge;
f←kludge + r/f;
ENDLOOP;
--unscaled result is (f+x)/(f-x)
x←(f+x)/(f-x);
fl ← LOOPHOLE[x,FloatNum];
--scale: multipy by 2↑k by adding k to the exponent
fl.exp← fl.exp+LowHalf[k];
RETURN[LOOPHOLE[fl,REAL]];
END;
Log: PUBLIC PROCEDURE [base,arg: REAL] RETURNS [REAL] =
BEGIN
--logxy = Ln(y)/Ln(x)
x,y: REAL;
x ← Ln[base];
y ← Ln[arg];
RETURN [y/x];
END;
Ln: PUBLIC PROCEDURE [x: REAL] RETURNS [REAL] =
BEGIN
fl: FloatNum;
Ln2,half: REAL;
i: CARDINAL;
m: INTEGER;
kludge: REAL;
OneMinusX,OneMinusXPowerRep,PreviousVal: REAL;
IF x<0 THEN
BEGIN
RETURN[-1];
END;
fl ← [sign: 0,exp: 200B,m1: 61B,m2: 71027B];
Ln2←LOOPHOLE[fl,REAL];
fl ← [sign: 0,exp: 200B,m1: 0,m2: 0];
half ← LOOPHOLE[fl,REAL];
--if x = z * 2↑m then Ln(x) = m * Ln2 + Ln(z). Here, 1/2 <= z < 1
fl ← LOOPHOLE[x,FloatNum];
m←fl.exp - 200B;
fl.exp ← 200B;--now 1/2 <= z < 1
x←LOOPHOLE[fl,REAL];
--special case for limit condition (x a power of 2)
IF x=half THEN
BEGIN
m←m-1;
kludge←m;
x←Ln2*kludge;
RETURN [x];
END;
OneMinusX ← 1 - x;
x←m*Ln2;
OneMinusXPowerRep ← 1;
FOR i IN [1..100] DO
PreviousVal←x;
OneMinusXPowerRep ← OneMinusXPowerRep*OneMinusX;
x←x-OneMinusXPowerRep/i;
IF x=PreviousVal THEN EXIT;
ENDLOOP;
RETURN[x];
END;
SqRt: PUBLIC PROCEDURE [x: REAL] RETURNS [REAL] =
BEGIN
y,x1,epsilon: REAL;
fl: FloatNum;
IF x <= 0 THEN RETURN[0];
fl ← LOOPHOLE[x,FloatNum];
fl.exp ← (fl.exp - 200B)/2 +200B;
y ← LOOPHOLE[fl,REAL];--initial guess = half original exponent
fl ← [sign: 0,exp: magic,m1: 0,m2: 0];
epsilon ← LOOPHOLE[fl,REAL];
DO
x1 ← y*y;
y←y-(x1-x)/(2*y);--slope of the curve here is 1/2
IF ABS[(x1-x)/x] <= ABS[epsilon] THEN EXIT;
ENDLOOP;
RETURN[y];
END;
Root: PUBLIC PROCEDURE [index,arg: REAL] RETURNS [REAL] =
BEGIN
--yth root of x = e↑(Ln(x)/y)
RETURN [Exp[Ln[arg]/index]];
END;
Power: PUBLIC PROCEDURE [base,exponent: REAL] RETURNS [REAL] =
BEGIN
--x↑y = e↑(y*Ln(x))
x: REAL;
IF base = 0 THEN IF exponent = 0 THEN RETURN[1] ELSE RETURN[0];
x ← Ln[base];
RETURN [Exp[exponent*x]];
END;
Cos: PUBLIC PROCEDURE [radians:REAL] RETURNS [cos:REAL] =
-- This function is good to 7.33 decimal places and is taken from:
-- Computer Approximations, John F. Hart et.al. pp118(top 2nd col).
BEGIN
negsign:REAL←1;
nangle:REAL←ABS[radians-RealDefs.Fix[radians*recpi]*twoPI];
SELECT nangle FROM
IN [PIovr2..PI3ovr2) => BEGIN negsign←-1; nangle←ABS[PI-nangle]; END;
IN [PI3ovr2.. twoPI) => nangle ← twoPI - nangle;
ENDCASE;
nangle←nangle*nangle;
cos←negsign*(cp0+nangle*(cp1+nangle*(cp2+nangle*(cp3+nangle*cp4))));
END;
CosDeg: PUBLIC PROCEDURE [degrees:REAL] RETURNS [cos:REAL] =
BEGIN
radians:REAL←degrees*degtorad;
cos←Cos[radians];
END;
Sin: PUBLIC PROCEDURE [radians:REAL] RETURNS [sin:REAL] =
BEGIN
radians←PIovr2 - radians;
sin ← Cos[radians];
END;
SinDeg: PUBLIC PROCEDURE [degrees:REAL] RETURNS [sin:REAL] =
BEGIN
radians:REAL←degrees*degtorad;
sin←Sin[radians];
END;
Tan: PUBLIC PROCEDURE [radians:REAL] RETURNS [tan:REAL] =
BEGIN
tan ← Sin[radians]/Cos[radians];
END;
TanDeg: PUBLIC PROCEDURE [degrees:REAL] RETURNS [tan:REAL] =
BEGIN
radians:REAL←degrees*degtorad;
tan←Sin[radians]/Cos[radians];
END;
ArcTan: PUBLIC PROCEDURE [y,x:REAL] RETURNS [radians:REAL] =
-- This function is good to 8.7 decimal places and is taken from:
-- Computer Approximations, John F. Hart et.al. pp129(top 2nd col).
BEGIN
s,v,t,t2,c,q:REAL;
IF ABS[x]<=ABS[smallnumber] THEN IF ABS[y]<=ABS[smallnumber] THEN RETURN [0] ELSE IF y<0 THEN RETURN[-PIovr2] ELSE RETURN[PIovr2];
IF x<0 THEN
IF y<0 THEN BEGIN
q←-PI;
s←1;
END ELSE
BEGIN
q←PI;
s←-1;
END
ELSE
IF y<0 THEN BEGIN
q←0;
s←-1;
END ELSE
BEGIN
q←0;
s←1;
END;
v←ABS [y/x];
SELECT v FROM
IN [zero..tanPI16)=>BEGIN t←v;
c←0;
END;
IN [tanPI16..tan3PI16)=>BEGIN t←(x2)-((x22)/(x2+v));
c←PIovr8;
END;
IN [tan3PI16..tan5PI16)=>BEGIN t←(x4)-((x44)/(x4+v));
c←PIovr4;
END;
IN [tan5PI16..tan7PI16)=>BEGIN t←(x6)-((x66)/(x6+v));
c←PI3ovr8;
END;
>= tan7PI16=>BEGIN t←-1/v;
c←PIovr2;
END;
ENDCASE;
t2←t*t;
radians←s*(t*(p0+t2*(p1+t2*(p2+t2*p3)))+c)+q;
END;
ArcTanDeg: PUBLIC PROCEDURE [y,x:REAL] RETURNS [degrees:REAL] =
BEGIN
radians:REAL←ArcTan[y,x];
degrees←radians*radtodeg;
END;
smallnumber:REAL←RealDefs.StringToFloat["0.00000000005"L];
zero:REAL←0;
PI:REAL←RealDefs.StringToFloat["3.1415926535"L];
twoPI:REAL←PI*2;
PI3ovr2:REAL←3*PI/2;
PIovr2:REAL←RealDefs.StringToFloat["1.5707963267"L];
PI3ovr8:REAL←RealDefs.StringToFloat["1.1780972"L];
PIovr4:REAL←RealDefs.StringToFloat[".7853981633"L];
PIovr8:REAL←RealDefs.StringToFloat[".392699"L];
cp0:REAL←RealDefs.StringToFloat[".9999999534"L];
cp1:REAL←-RealDefs.StringToFloat[".4999990534"L];
cp2:REAL←RealDefs.StringToFloat[".04166358467"L];
cp3:REAL←-RealDefs.StringToFloat[".001385370426"L];
cp4:REAL←RealDefs.StringToFloat[".00002315393167"L];
radtodeg:REAL←180/PI;
degtorad:REAL←PI/180;
recpi:REAL←1/twoPI;
tanPI16:REAL←RealDefs.StringToFloat[".1989123673"L];
x2:REAL←1/(RealDefs.StringToFloat[".4142135623"L]);
x22:REAL ← x2*x2+1;
tan3PI16:REAL←RealDefs.StringToFloat[".6681786379"L];
x4:REAL←1/(RealDefs.StringToFloat["1."L]);
x44:REAL ← x4*x4+1;
tan5PI16:REAL←RealDefs.StringToFloat["1.4966057626"L];
x6:REAL←1/(RealDefs.StringToFloat["2.4142135623"L]);
x66:REAL ← x6*x6+1;
tan7PI16:REAL←RealDefs.StringToFloat["5.0273394921"L];
p0:REAL←RealDefs.StringToFloat[".9999999980228"L];
p1:REAL←-RealDefs.StringToFloat[".333331726035"L];
p2:REAL←RealDefs.StringToFloat[".1997952738"L];
p3:REAL←-RealDefs.StringToFloat[".134450639"L];
END.