<> <> <> <> DIRECTORY IeeeInternal USING [CVExtended, ExponentBias, SingleReal, Unpack], Real USING [Extended, Fix, RealError, RealException], RealFns; RealFnsImpl: CEDAR PROGRAM IMPORTS IeeeInternal, Real EXPORTS RealFns = BEGIN Precision: CARDINAL _ 5; LoopCount: INTEGER _ 3; Ln2: REAL = 0.693147181; LogBase2ofE: REAL = 1/Ln2; <> <> <> 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> i _ 4*LoopCount + 2; f _ i; FOR i DECREASING IN [1..LoopCount] DO ikludge _ 4*i - 2; kludge _ ikludge; f _ kludge + r/f; ENDLOOP; <> 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; <> x _ x*LOOPHOLE[fl, REAL]; RETURN[x]; }; Log: PUBLIC PROC [base, arg: REAL] RETURNS [REAL] = { <> 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[]; }; <> 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]; <> 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] = { <> RETURN[Exp[Ln[arg]/index]]; }; Power: PUBLIC PROC [base, exponent: REAL] RETURNS [REAL] = { <> 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] = { <> <> 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] = { <> <> 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. 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