-- 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
--log
xy = 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
--y
th 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;
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
END;

Sin
: PUBLIC PROCEDURE [radians:REAL] RETURNS [sin:REAL] =
BEGIN
END;

SinDeg
: PUBLIC PROCEDURE [degrees:REAL] RETURNS [sin:REAL] =
BEGIN
END;

Tan
: PUBLIC PROCEDURE [radians:REAL] RETURNS [tan:REAL] =
BEGIN
END;

TanDeg
: PUBLIC PROCEDURE [degrees:REAL] RETURNS [tan:REAL] =
BEGIN
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;
END;

ArcTanDeg
: PUBLIC PROCEDURE [y,x:REAL] RETURNS [degrees:REAL] =
BEGIN
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];
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.