RealFnsImpl:
CEDAR
PROGRAM
IMPORTS IeeeInternal, Real
EXPORTS RealFns =
BEGIN
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𡤄*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, k];
RETURN[x];
};
Pow2:
PROC [x:
REAL, exp:
INTEGER]
RETURNS [
REAL] =
INLINE {
fl: IeeeInternal.SingleReal ← LOOPHOLE[1.0, IeeeInternal.SingleReal];
IF exp NOT IN [-126..127] THEN ERROR;
fl.exp ← exp + IeeeInternal.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: IeeeInternal.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 ← IeeeInternal.CVExtended[IeeeInternal.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, IeeeInternal.SingleReal];
m ← fl.exp - (IeeeInternal.ExponentBias - 1);
fl.exp ← IeeeInternal.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: IeeeInternal.SingleReal;
IF x <= 0 THEN RETURN[0];
fl ← LOOPHOLE[x, IeeeInternal.SingleReal];
fl.exp ← (fl.exp - IeeeInternal.ExponentBias)/2 + IeeeInternal.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: IeeeInternal.SingleReal ← LOOPHOLE[x];
RETURN[fl.exp < (distance + IeeeInternal.ExponentBias)];
};
AlmostEqual:
PUBLIC
PROC [y, x:
REAL, distance:
INTEGER [-126..0]]
RETURNS [BOOLEAN] = {
fl: IeeeInternal.SingleReal ← LOOPHOLE[(x - y)];
oldexp:
INTEGER ←
MAX[
LOOPHOLE[x, IeeeInternal.SingleReal].exp, LOOPHOLE[y, IeeeInternal.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.