<> <> <> <> <> <> <> DIRECTORY Ieee USING [CVExtended, ExponentBias, SingleReal, Unpack], Real USING [Extended, Fix, FixC, FixI, RealError, RealException], RealFns; RealFnsImpl: CEDAR PROGRAM IMPORTS Ieee, Real EXPORTS RealFns = BEGIN <> <> Ln2: REAL = 0.693147181; LogBase2ofE: REAL _ 1/Ln2; 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; radiansPerDegree: REAL _ PI/180; degreesPerRadian: REAL _ 180/PI; rec2pi: REAL _ 1.0/twoPI; rec360: REAL _ 1.0/360.0; fourOverPI: REAL _ 4.0/PI; Exp: PUBLIC PROC [x: REAL] RETURNS [REAL] = { <> <> r, f: REAL; k: INTEGER; fl: Ieee.SingleReal _ LOOPHOLE[1.0, Ieee.SingleReal]; <> SELECT ABS[x] FROM < 3e-4 => RETURN [1+x]; < 3e-3 => RETURN [1+x+x*x*0.5]; ENDCASE; x _ x*LogBase2ofE; --exponent of 2 (instead of e) IF x < -Ieee.ExponentBias THEN RETURN [0.0]; IF x > Ieee.ExponentBias THEN ERROR Real.RealError[]; k _ Real.FixI[x]; --integer part x _ x - k; --fractional part x _ x*Ln2; --now result = 2^k * exp(x) and 0<=x> LoopCount: NAT = 3; f _ 4*LoopCount + 2; f _ (4*(LoopCount-0) - 2) + r/f; f _ (4*(LoopCount-1) - 2) + r/f; f _ (4*(LoopCount-2) - 2) + r/f; < 3, then use the following comments>> <> <> }; <> x _ (f + x)/(f - x); <> fl.exp _ k + Ieee.ExponentBias; RETURN[x*LOOPHOLE[fl, REAL]]; }; Log: PUBLIC PROC [base, arg: REAL] RETURNS [REAL] = { <> RETURN[Ln[arg]/Ln[base]]; }; Ln: PUBLIC PROC [x: REAL] RETURNS [REAL] = { fl: Ieee.SingleReal; m: INTEGER; xm1: REAL _ x-1; OneMinusX, OneMinusXPowerRep: REAL; IF x <= 0 THEN Bomb[x]; <> SELECT ABS[xm1] FROM > 3e-3 => {}; > 3e-4 => RETURN [xm1 - 0.5*xm1*xm1]; ENDCASE => RETURN [xm1]; <> fl _ LOOPHOLE[x, Ieee.SingleReal]; m _ fl.exp - (Ieee.ExponentBias - 1); fl.exp _ Ieee.ExponentBias - 1; --now 1/2 <= z < 1 x _ LOOPHOLE[fl, REAL]; <> IF x = 0.5 THEN RETURN[Ln2*(m-1)]; OneMinusX _ 1 - x; x _ m*Ln2; OneMinusXPowerRep _ 1; FOR i: CARDINAL _ 1, i+1 DO PreviousVal: REAL _ x; OneMinusXPowerRep _ OneMinusXPowerRep*OneMinusX; x _ x - OneMinusXPowerRep/i; IF x = PreviousVal THEN EXIT; ENDLOOP; RETURN[x]; }; Bomb: PROC [x: REAL] = { [] _ SIGNAL Real.RealException[flags: [invalidOperation: TRUE], vp: NEW[Real.Extended _ Ieee.CVExtended[Ieee.Unpack[x]]]]; ERROR Real.RealError[]; }; SquareReal: TYPE = MACHINE DEPENDENT RECORD [ m2: CARDINAL, sign: BOOL, expD2: [0..128), index: [0..32), m1: [0..8)]; <> SqRt: PUBLIC PROC [a: REAL] RETURNS [b: REAL _ 0.0] = { <> SELECT TRUE FROM a <= 0.0 => { <> IF a < 0.0 THEN Bomb[a]; b _ 0.0; }; LOOPHOLE[a, Ieee.SingleReal].exp = 0 => <> b _ SqRt[a * (LONG[20000B] * LONG[20000B])] / 20000B; ENDCASE => { aFmt: SquareReal _ LOOPHOLE[a, SquareReal]; bFmt: Ieee.SingleReal _ LOOPHOLE[sqGuesses[aFmt.index], Ieee.SingleReal]; <> bFmt.exp _ bFmt.exp + aFmt.expD2 - 63; <> b _ LOOPHOLE[bFmt, REAL]; <> b _ LOOPHOLE[LOOPHOLE[b + a/b, LONG CARDINAL] - 40000000B, REAL]; <> b _ LOOPHOLE[LOOPHOLE[b + a/b, LONG CARDINAL] - 40000000B, REAL]; <> }; }; sqGuesses: ARRAY[0..32) OF REAL _ [ LOOPHOLE[7715741440B, REAL], -- 0.7178211 LOOPHOLE[7717240700B, REAL], LOOPHOLE[7720514140B, REAL], LOOPHOLE[7721745140B, REAL], LOOPHOLE[7723155200B, REAL], LOOPHOLE[7724346000B, REAL], LOOPHOLE[7725517500B, REAL], LOOPHOLE[7726654000B, REAL], -- 0.8568115 LOOPHOLE[7727773400B, REAL], -- 0.8748627 LOOPHOLE[7731077100B, REAL], LOOPHOLE[7732167300B, REAL], LOOPHOLE[7733245000B, REAL], LOOPHOLE[7734310400B, REAL], LOOPHOLE[7735342500B, REAL], LOOPHOLE[7736363400B, REAL], LOOPHOLE[7737373700B, REAL], -- 0.9920616 LOOPHOLE[7740370240B, REAL], -- 1.015156 LOOPHOLE[7741351440B, REAL], LOOPHOLE[7742314540B, REAL], LOOPHOLE[7743243000B, REAL], LOOPHOLE[7744155200B, REAL], LOOPHOLE[7745054400B, REAL], LOOPHOLE[7745741300B, REAL], LOOPHOLE[7746614600B, REAL], -- 1.211716 LOOPHOLE[7747457040B, REAL], -- 1.237247 LOOPHOLE[7750310640B, REAL], LOOPHOLE[7751132500B, REAL], LOOPHOLE[7751745000B, REAL], LOOPHOLE[7752550100B, REAL], LOOPHOLE[7753344400B, REAL], LOOPHOLE[7754132600B, REAL], LOOPHOLE[7754712500B, REAL] -- 1.402992 ]; Root: PUBLIC PROC [index, arg: REAL] RETURNS [REAL] = { <> RETURN[Exp[Ln[arg]/index]]; }; Power: PUBLIC PROC [base, exponent: REAL] RETURNS [REAL] = { <> IF base = 0 THEN IF exponent = 0 THEN RETURN[1] ELSE RETURN[0]; RETURN[Exp[exponent*Ln[base]]]; }; Sin: PUBLIC PROC [radians: REAL] RETURNS [sin: REAL] = { rem: REAL _ ABS[radians]; IF rem > twoPI THEN rem _ GetWithinTwoPI[rem]; <> { <> SELECT Real.FixC[rem * fourOverPI] FROM 0 => {GO TO posSin}; 1 => {rem _ PIovr2 - rem; GO TO posCos}; 2 => {rem _ rem - PIovr2; GO TO posCos}; 3 => {rem _ PI - rem; GO TO posSin}; 4 => {rem _ rem - PI; GO TO negSin}; 5 => {rem _ PI3ovr2 - rem; GO TO negCos}; 6 => {rem _ rem - PI3ovr2; GO TO negCos}; 7 => {rem _ twoPI - rem; GO TO negSin}; ENDCASE => RETURN [0.0]; EXITS posSin => {sin _ SinFirstOctant[rem]; IF radians < 0.0 THEN GO TO neg}; negSin => {sin _ SinFirstOctant[rem]; IF radians > 0.0 THEN GO TO neg}; posCos => {sin _ CosFirstOctant[rem]; IF radians < 0.0 THEN GO TO neg}; negCos => {sin _ CosFirstOctant[rem]; IF radians > 0.0 THEN GO TO neg}; }; EXITS neg => sin _ - sin; }; Cos: PUBLIC PROC [radians: REAL] RETURNS [cos: REAL] = { rem: REAL _ ABS[radians]; IF rem > twoPI THEN rem _ GetWithinTwoPI[rem]; <> { <> SELECT Real.FixC[rem * fourOverPI] FROM 0 => {GO TO posCos}; 1 => {rem _ PIovr2 - rem; GO TO posSin}; 2 => {rem _ rem - PIovr2; GO TO negSin}; 3 => {rem _ PI - rem; GO TO negCos}; 4 => {rem _ rem - PI; GO TO negCos}; 5 => {rem _ PI3ovr2 - rem; GO TO negSin}; 6 => {rem _ rem - PI3ovr2; GO TO posSin}; 7 => {rem _ twoPI - rem; GO TO posCos}; ENDCASE => cos _ 1.0; EXITS posSin => {cos _ SinFirstOctant[rem]}; negSin => {cos _ SinFirstOctant[rem]; cos _ -cos}; posCos => {cos _ CosFirstOctant[rem]}; negCos => {cos _ CosFirstOctant[rem]; cos _ -cos}; }; }; SinDeg: PUBLIC PROC [degrees: REAL] RETURNS [sin: REAL] = { IF ABS[degrees] > 360 THEN degrees _ GetWithin360[degrees]; <> sin _ Sin[degrees*radiansPerDegree]; }; CosDeg: PUBLIC PROC [degrees: REAL] RETURNS [cos: REAL] = { degrees _ ABS[degrees]; IF degrees > 360 THEN degrees _ GetWithin360[degrees]; <> cos _ Cos[degrees*radiansPerDegree]; }; Tan: PUBLIC PROC [radians: REAL] RETURNS [tan: REAL] = { IF ABS[radians] > twoPI THEN radians _ GetWithinTwoPI[radians]; tan _ Sin[radians] / Cos[radians]; }; TanDeg: PUBLIC PROC [degrees: REAL] RETURNS [tan: REAL] = { IF ABS[degrees] > 360 THEN degrees _ GetWithin360[degrees]; tan _ SinDeg[degrees] / CosDeg[degrees]; }; <> tanPI16: REAL = .1989123673; x2: REAL _ 1/(.414213562); x22: REAL _ x2*x2 + 1; tan3PI16: REAL = .668178638; x4: REAL = 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; 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 < tanPI16 => {t _ v; c _ 0}; < tan3PI16 => {t _ x2 - x22/(x2 + v); c _ PIovr8}; < tan5PI16 => {t _ x4 - x44/(x4 + v); c _ PIovr4}; < tan7PI16 => {t _ x6 - x66/(x6 + v); c _ PI3ovr8}; ENDCASE => {t _ -1/v; c _ PIovr2}; t2 _ t*t; radians _ s*(t*(p0 + t2*(p1 + t2*(p2 + t2*p3))) + c) + q; }; ArcTanDeg: PUBLIC PROC [y, x: REAL] RETURNS [degrees: REAL] = { degrees _ ArcTan[y, x]*degreesPerRadian; }; AlmostZero: PUBLIC PROC [x: REAL, distance: INTEGER [-126..127]] RETURNS [BOOL] = { RETURN[LOOPHOLE[x, Ieee.SingleReal].exp < (distance + Ieee.ExponentBias)]; }; AlmostEqual: PUBLIC PROC [y, x: REAL, distance: INTEGER [-126..0]] RETURNS [BOOL] = { fl: Ieee.SingleReal _ LOOPHOLE[(x - y)]; oldexp: INTEGER _ MAX[ LOOPHOLE[x, Ieee.SingleReal].exp, LOOPHOLE[y, Ieee.SingleReal].exp]; RETURN[fl.exp < (distance + oldexp)]; }; <> <> fact9recip: REAL _ 1.0/(INT[1]*2*3*4*5*6*7*8*9); fact8recip: REAL _ 1.0/(INT[1]*2*3*4*5*6*7*8); fact7recip: REAL _ 1.0/(INT[1]*2*3*4*5*6*7); fact6recip: REAL _ 1.0/(INT[1]*2*3*4*5*6); fact5recip: REAL _ 1.0/(INT[1]*2*3*4*5); fact4recip: REAL _ 1.0/(INT[1]*2*3*4); fact3recip: REAL _ 1.0/(INT[1]*2*3); fact2recip: REAL _ 1.0/(INT[1]*2); minEps: REAL _ 1.0 - LOOPHOLE[LOOPHOLE[1.0, LONG CARDINAL]-1, REAL]; <> sinCosEps: REAL _ SqRt[minEps]; <> SinFirstOctant: PROC [radians: REAL] RETURNS [REAL] = { <> IF radians <= sinCosEps THEN RETURN [radians] ELSE { r2: REAL ~ radians*radians; RETURN [(((fact9recip*r2-fact7recip)*r2+fact5recip)*r2-fact3recip)*r2*radians+radians]; }; }; CosFirstOctant: PROC [radians: REAL] RETURNS [REAL] = { <> IF radians <= sinCosEps THEN RETURN [1.0] ELSE { r2: REAL ~ radians*radians; RETURN [((((fact8recip*r2-fact6recip)*r2+fact4recip)*r2-fact2recip)*r2+1)]; }; }; <> PIc1: REAL _ 6.28125; -- 402.0/64 (exact representation near twoPI) PIc2: REAL = 1.9353072e-3; -- 2*PI - c1 GetWithinTwoPI: PROC [radians: REAL] RETURNS [REAL] = { <> <> <= 0 >> xn: REAL _ Real.Fix[radians * rec2pi]; RETURN [(radians - xn*PIc1) - xn*PIc2]; }; GetWithin360: PROC [degrees: REAL] RETURNS [REAL] = { <> <> <= 0 >> xn: REAL _ Real.Fix[degrees * rec360]; RETURN [degrees - xn*360]; }; 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 October 24, 1984 4:37:17 pm PDT, Plass; New impls for Sin, Cos, SinDeg, CosDeg November 2, 1984, Atkinson; Faster & more accurate SqRt, Sin, Cos, SinDeg, CosDeg