-- RealFnsImpl.mesa
-- Last modified: Stewart on August 27, 1982 1:07 pm
DIRECTORY
Ieee USING [CVExtended, ExponentBias, SingleReal, Unpack],
Real USING [Extended, Fix, RealError, RealException],
Inline USING [LowHalf],
RealFns;
RealFnsImpl: CEDAR PROGRAM IMPORTS Ieee, Inline, Real EXPORTS RealFns =
BEGIN OPEN Ieee;
Precision: CARDINAL ← 5;
LoopCount: INTEGER ← 3;
Ln2: REAL = 0.693147181;
LogBase2ofE: REAL = 1/Ln2;
--compute e↑n using continued fractions
--See ACM Algorithm 229, Morelock, James C., "Elementary functions
-- by continued fractions"
Exp: PUBLIC PROC [x: REAL] RETURNS [REAL] = {
r, f, kludge: REAL;
k: LONG INTEGER;
i, ikludge: INTEGER;
IF x = 0 THEN RETURN[1];
x ← x*LogBase2ofE; --exponent of 2 (instead of e)
k ← Real.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); --scale: multipy by 2↑k by adding k to the exponent
x ← Pow2[x, LOOPHOLE[Inline.LowHalf[LOOPHOLE[k]]]];
RETURN[x];
};
Pow2: PROC [x: REAL, exp: INTEGER] RETURNS [REAL] = INLINE {
fl: SingleReal ← LOOPHOLE[1.0, SingleReal];
IF exp NOT IN [-126..127] THEN ERROR;
fl.exp ← exp + ExponentBias;
--scale: multipy by 2↑exp by adding exp to the exponent
x ← x*LOOPHOLE[fl, REAL];
RETURN[x];
};
Log: PUBLIC PROC [base, arg: REAL] RETURNS [REAL] = {
--logxy = Ln(y)/Ln(x)
x, y: REAL;
x ← Ln[base];
y ← Ln[arg];
RETURN[y/x];
};
Ln: PUBLIC PROC [x: REAL] RETURNS [REAL] = {
fl: SingleReal;
i: CARDINAL;
m: INTEGER;
kludge: REAL;
OneMinusX, OneMinusXPowerRep, PreviousVal: REAL;
IF x < 0 THEN {
[] ← SIGNAL Real.RealException[flags: [invalidOperation: TRUE], vp: NEW[Real.Extended ← Ieee.CVExtended[Ieee.Unpack[x]]]];
ERROR Real.RealError[];
};
--if x = z * 2↑m then Ln(x) = m * Ln2 + Ln(z). Here, 1/2 <= z < 1
fl ← LOOPHOLE[x, SingleReal];
m ← fl.exp - (ExponentBias - 1);
fl.exp ← ExponentBias - 1; --now 1/2 <= z < 1
x ← LOOPHOLE[fl, REAL];
--special case for limit condition (x a power of 2)
IF x = 0.5 THEN {m ← m - 1; kludge ← m; x ← Ln2*kludge; RETURN[x]; };
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];
};
SqRt: PUBLIC PROC [x: REAL] RETURNS [REAL] = {
epsilon: REAL = Pow2[1, -23];
y, x1: REAL;
fl: SingleReal;
IF x <= 0 THEN RETURN[0];
fl ← LOOPHOLE[x, SingleReal];
fl.exp ← (fl.exp - ExponentBias)/2 + ExponentBias;
y ← LOOPHOLE[fl, REAL]; --initial guess = half original exponent
DO
x1 ← y*y;
y ← y - (x1 - x)/(2*y); --slope of the curve here is 1/2
IF ABS[(x1 - x)/x] <= epsilon THEN EXIT; -- may need ABS[]
ENDLOOP;
RETURN[y];
};
Root: PUBLIC PROC [index, arg: REAL] RETURNS [REAL] = {
--yth root of x = e↑(Ln(x)/y)
RETURN[Exp[Ln[arg]/index]];
};
Power: PUBLIC PROC [base, exponent: REAL] RETURNS [REAL] = {
--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]];
};
Cos: PUBLIC PROC [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).
negsign: REAL ← 1;
nangle: REAL ← ABS[radians - Real.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))));
};
CosDeg: PUBLIC PROC [degrees: REAL] RETURNS [cos: REAL] = {
radians: REAL ← degrees*degtorad; cos ← Cos[radians]; };
Sin: PUBLIC PROC [radians: REAL] RETURNS [sin: REAL] = {
radians ← PIovr2 - radians; sin ← Cos[radians]; };
SinDeg: PUBLIC PROC [degrees: REAL] RETURNS [sin: REAL] = {
radians: REAL ← degrees*degtorad; sin ← Sin[radians]; };
Tan: PUBLIC PROC [radians: REAL] RETURNS [tan: REAL] = {
tan ← Sin[radians]/Cos[radians]; };
TanDeg: PUBLIC PROC [degrees: REAL] RETURNS [tan: REAL] = {
radians: REAL ← degrees*degtorad; tan ← Sin[radians]/Cos[radians]; };
ArcTan: PUBLIC PROC [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).
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 {q ← -PI; s ← 1; } ELSE {q ← PI; s ← -1; }
ELSE IF y < 0 THEN {q ← 0; s ← -1; } ELSE {q ← 0; s ← 1; };
v ← ABS[y/x];
SELECT v FROM
IN [0..tanPI16) => {t ← v; c ← 0; };
IN [tanPI16..tan3PI16) => {t ← (x2) - ((x22)/(x2 + v)); c ← PIovr8; };
IN [tan3PI16..tan5PI16) => {t ← (x4) - ((x44)/(x4 + v)); c ← PIovr4; };
IN [tan5PI16..tan7PI16) => {t ← (x6) - ((x66)/(x6 + v)); c ← PI3ovr8; };
>= tan7PI16 => {t ← -1/v; c ← PIovr2; };
ENDCASE;
t2 ← t*t;
radians ← s*(t*(p0 + t2*(p1 + t2*(p2 + t2*p3))) + c) + q;
};
ArcTanDeg: PUBLIC PROC [y, x: REAL] RETURNS [degrees: REAL] = {
radians: REAL ← ArcTan[y, x]; degrees ← radians*radtodeg; };
AlmostZero: PUBLIC PROC [x: REAL, distance: INTEGER [-126..127]]
RETURNS [BOOLEAN] = {
fl: Ieee.SingleReal ← LOOPHOLE[x];
RETURN[fl.exp < (distance + Ieee.ExponentBias)];
};
AlmostEqual: PUBLIC PROC [y, x: REAL, distance: INTEGER [-126..0]]
RETURNS [BOOLEAN] = {
fl: Ieee.SingleReal ← LOOPHOLE[(x - y)];
oldexp: INTEGER ← MAX[
LOOPHOLE[x, Ieee.SingleReal].exp, LOOPHOLE[y, Ieee.SingleReal].exp];
RETURN[fl.exp < (distance + oldexp)];
};
smallnumber: REAL = 0.00000000005;
PI: REAL = 3.14159265;
twoPI: REAL = PI*2;
PI3ovr2: REAL = 3*PI/2;
PIovr2: REAL = 1.570796327;
PI3ovr8: REAL = 1.1780972;
PIovr4: REAL = .785398163;
PIovr8: REAL = .392699;
cp0: REAL = .999999953;
cp1: REAL = -.499999053;
cp2: REAL = .0416635847;
cp3: REAL = -.001385370426;
cp4: REAL = .0000231539317;
radtodeg: REAL = 180/PI;
degtorad: REAL = PI/180;
recpi: REAL = 1/twoPI;
tanPI16: REAL = .1989123673;
x2: REAL = 1/(.414213562);
x22: REAL = x2*x2 + 1;
tan3PI16: REAL = .668178638;
x4: REAL = 1/(1.0);
x44: REAL = x4*x4 + 1;
tan5PI16: REAL = 1.496605763;
x6: REAL = 1/(2.41421356);
x66: REAL = x6*x6 + 1;
tan7PI16: REAL = 5.02733949;
p0: REAL = .999999998;
p1: REAL = -.333331726;
p2: REAL = .1997952738;
p3: REAL = -.134450639;
END.
September 16, 1980 3:13 PM, Stewart; Bug in Exp
September 28, 1980 8:28 PM, Stewart; Add AlmostEqual and AlmostZero, format
November 7, 1980 3:33 PM, Stewart; fix AlmostZero
August 27, 1982 1:07 pm, Stewart; CEDAR