<<>> <> <> <> <> <> <> <> <> <> <> <> <> DIRECTORY Basics, FloatingPointCommon, Real, RealSupport USING [Fix], RealFns; RealFnsImpl: CEDAR PROGRAM IMPORTS Basics, FloatingPointCommon, Real, RealSupport EXPORTS RealFns = BEGIN SingleReal: TYPE = MACHINE DEPENDENT RECORD [ sign (0: 0..0): BOOL, exp (0: 1..8): Exponent, m (0: 9..31): Mantissa]; Exponent: TYPE = CARDINAL [0..256); Mantissa: TYPE = CARDINAL [0..HiddenBit); HiddenBit: CARDINAL = 800000h; ExponentBias: INTEGER = 127; <> DenormExponent: INTEGER = Exponent.FIRST; <> NaNExponent: INTEGER = Exponent.LAST; <> SquareReal: TYPE = MACHINE DEPENDENT RECORD [ sign (0: 0..0): BOOL, expD2 (0: 1..7): CARDINAL [0..128), index (0: 8..12): CARDINAL [0..32), m (0: 13..31): CARDINAL [0..2000000B)]; <> Bomb: PROC [x: REAL, exCode: FloatingPointCommon.Exception ¬ invalidOperation] = { ERROR FloatingPointCommon.Error[exCode]; }; <> Ln2: REAL = 0.693147181; LogBase2ofE: REAL ¬ 1/Ln2; smallnumber: REAL = 0.00000000005; PI: REAL = 3.1415926535; twoPI: REAL ¬ PI*2; piHalves: REAL ¬ PI/2; pi3ovr2: REAL ¬ 3*piHalves; pi3Eighths: REAL = pi3ovr2/4; piFourths: REAL = PI/4; piEighths: REAL ¬ PI/8; radiansPerDegree: REAL ¬ PI/180; degreesPerRadian: REAL ¬ 180/PI; rec2pi: REAL ¬ 1.0/twoPI; rec360: REAL ¬ 1.0/360.0; fourOverPI: REAL ¬ 4.0/PI; twoOverPI: REAL ¬ 2.0/PI; invRoot2: REAL ¬ 0.70710678119; lnRoot2pi: REAL ¬ 0.9189385332; sqrtEps: REAL ¬ Real.FScale[1.0, -12]; -- 2­-12 = square root of machine epsilon Exp: PUBLIC PROC [x: REAL] RETURNS [REAL] = { <> <> r, f: REAL; k: INTEGER; fl: SingleReal ¬ LOOPHOLE[1.0, 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 < -ExponentBias THEN RETURN [0.0]; IF x > ExponentBias THEN Bomb[x]; k ¬ RealSupport.Fix[x]; --integer part (was InlineFixI) 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 + 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: 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, SingleReal]; m ¬ fl.exp - (ExponentBias - 1); fl.exp ¬ 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]; }; 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, SingleReal].exp = 0 => <> b ¬ SqRt[a * (LONG[20000B] * LONG[20000B])] / 20000B; ENDCASE => { aFmt: SquareReal ¬ LOOPHOLE[a, SquareReal]; bFmt: SingleReal ¬ LOOPHOLE[sqGuesses[aFmt.index], SingleReal]; <> aExp: CARDINAL ¬ aFmt.expD2; bFmt.exp ¬ bFmt.exp + aExp - 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 CARD = [ 7715741440B, 7717240700B, 7720514140B, 7721745140B, -- 0.71782 .. 0.78043 7723155200B, 7724346000B, 7725517500B, 7726654000B, -- 0.80021 .. 0.85681 7727773400B, 7731077100B, 7732167300B, 7733245000B, -- 0.87486 .. 0.92691 7734310400B, 7735342500B, 7736363400B, 7737373700B, -- 0.94362 .. 0.99206 7740370240B, 7741351440B, 7742314540B, 7743243000B, -- 1.01516 .. 1.10370 7744155200B, 7745054400B, 7745741300B, 7746614600B, -- 1.13167 .. 1.21172 7747457040B, 7750310640B, 7751132500B, 7751745000B, -- 1.23725 .. 1.31085 7752550100B, 7753344400B, 7754132600B, 7754712500B -- 1.33448 .. 1.40299 ]; 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 RealSupport.Fix[rem * fourOverPI] FROM -- was InlineFixC 0 => {GO TO posSin}; 1 => {rem ¬ piHalves - rem; GO TO posCos}; 2 => {rem ¬ rem - piHalves; 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 RealSupport.Fix[rem * fourOverPI] FROM -- was InlineFixC 0 => {GO TO posCos}; 1 => {rem ¬ piHalves - rem; GO TO posSin}; 2 => {rem ¬ rem - piHalves; 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] = { RETURN[TanCot[radians, FALSE]]; }; CoTan: PUBLIC PROC [radians: REAL] RETURNS [cotan: REAL] = { RETURN[TanCot[radians, TRUE]]; }; TanCot: PUBLIC PROC [x: REAL, cot: BOOLEAN] RETURNS [tan: REAL] = { c1: REAL = 201.0/128.0; -- exactly representable c2: REAL = .000483826794897; -- c1 + c2 = pi/2 xn, f, g: REAL; num, den: REAL; n: INT; p1: REAL = -.0958017723; q1: REAL = -0.429135777; q2: REAL = .00971685835; n ¬ Real.Round[twoOverPI*x]; xn ¬ REAL[n]; f ¬ (x - xn*c1) - xn*c2; <> IF ABS[f] < sqrtEps THEN { num ¬ f; den ¬ 1.0; } ELSE { g ¬ f*f; num ¬ f + f*g*p1; den ¬ (q2*g + q1)*g + 1; }; IF cot THEN {IF Basics.OddInt[n] THEN RETURN[-num/den] ELSE RETURN[den/num]} ELSE {IF Basics.OddInt[n] THEN RETURN[-den/num] ELSE RETURN[num/den]}; }; TanDeg: PUBLIC PROC [degrees: REAL] RETURNS [tan: REAL] = { IF ABS[degrees] > 360 THEN degrees ¬ GetWithin360[degrees]; tan ¬ SinDeg[degrees] / CosDeg[degrees]; }; CoTanDeg: PUBLIC PROC [degrees: REAL] RETURNS [cotan: REAL] = { IF ABS[degrees] > 360 THEN degrees ¬ GetWithin360[degrees]; cotan ¬ CoTan[degrees*radiansPerDegree]; }; ArcSin: PUBLIC PROC [x: REAL] RETURNS [ REAL] = { z: REAL; SELECT (z ¬ ABS[x]) FROM < sqrtEps => RETURN[x]; -- so small that arcsin(x) = x + x­3/6 + x­5/40 + ... < 0.5 => { RETURN[ArcTan[x, SqRt[1 - z*z]]]; }; ENDCASE => { RETURN[ArcTan[x, SqRt[2*(1-z) - (1-z)*(1-z)]]]; }; }; ArcCos: PUBLIC PROC [x: REAL] RETURNS [REAL] = { RETURN[IF x = -1.0 THEN PI ELSE 2.0*ArcTan[SqRt[(1.0-x)/(1.0+x)], 1.0]]; <> }; ArcSinDeg: PUBLIC PROC [x: REAL] RETURNS [degrees: REAL] = { degrees ¬ ArcSin[x]*degreesPerRadian; }; ArcCosDeg: PUBLIC PROC [x: REAL] RETURNS [degrees: REAL] = { degrees ¬ ArcCos[x]*degreesPerRadian; }; << David Goldberg notes regarding the accuracy of two methods for computing ArcCos, where Method 1 is ArcTan[SqRt[1.0-x*x], x], and Method 2 is the one above: When x ~ y, then x-y is dangerous if x and/or y have suffered roundoff error, because all the significant digits may be cancelled, leaving only roundoff error. That's the trouble with 1-x*x; there will be roundoff error in x*x, and if x ~ 1, then 1-x*x can have a large relative error. However, if x ~ 1, (in fact, if 1/2 < x < 2), then 1-x is exact because there is no roundoff error. Thus 1/(1-x) is very accurate. That's the theory according to floating point experts. Method 1 is worse when x = 1-n*2^(-24), n = sqrt(2^(23)), i.e., x = 0.9998274. The correct value is 1.857983e-02; Method 1 gives 1.858144e-2 (error .000161e-2) while Method 2 gives 1.858064e-2 (error .000081e-2). So, in this case, Method 2 has half the error of Method 1.>> << [Artwork node; type 'Artwork on' to command tool] >> <> <<>> <> 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[-piHalves] ELSE RETURN[piHalves]; 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 ¬ piEighths}; < tan5PI16 => {t ¬ x4 - x44/(x4 + v); c ¬ piFourths}; < tan7PI16 => {t ¬ x6 - x66/(x6 + v); c ¬ pi3Eighths}; ENDCASE => {t ¬ -1/v; c ¬ piHalves}; 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; }; SinH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { ex: REAL = Exp[x]; emx: REAL = 1/ex; f ¬ (ex-emx)/2.0 }; CosH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { ex: REAL = Exp[x]; emx: REAL = 1/ex; f ¬ (ex+emx)/2.0 }; TanH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { IF ABS[x] > 20. THEN f ¬ (IF x<0. THEN -1. ELSE 1.) ELSE { ex: REAL = Exp[x]; emx: REAL = 1/ex; f ¬ (ex-emx)/(ex+emx); }; }; CotH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { IF ABS[x] > 20. THEN f ¬ (IF x<0. THEN -1. ELSE 1.) ELSE { ex: REAL = Exp[x]; emx: REAL = 1/ex; f ¬ (ex+emx)/(ex-emx); }; }; InvSinH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { f ¬ Ln[x+SqRt[x*x+1.0]]; }; InvCosH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { IF x < 1.0 THEN Bomb[x]; f ¬ Ln[x+SqRt[x*x-1.0]]; }; InvTanH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { IF ABS[x] >= 1.0 THEN Bomb[x]; f ¬ 0.5*Ln[(1.0+x)/(1.0-x)]; }; InvCotH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { IF ABS[x] <= 1.0 THEN Bomb[x]; f ¬ 0.5*Ln[(x+1.0)/(x-1.0)]; }; <> <<>> LnGamma: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { phi: REAL; lnAdjFact: REAL ¬ 0; IF x <= 0.0 THEN Bomb[x]; IF x < 8 THEN { -- Gamma[x+1] = x! = x*(x-1)! = x*Gamma[x] adjFact: REAL ¬ 1.0; WHILE x < 8 DO adjFact ¬ adjFact*x; x ¬ x+1.0; ENDLOOP; lnAdjFact ¬ Ln[adjFact]; }; phi ¬ (-0.2762472E-2/(x*x)+0.8333327385E-1)/x; -- LGAM 5400 - 8.82 digits f ¬ (x-0.5)*Ln[x]-x+lnRoot2pi+phi-lnAdjFact; -- Stirling's approx. }; Gamma: PUBLIC PROC [x: REAL] RETURNS [REAL] = { RETURN[Exp[LnGamma[x]]]; }; J0: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { <> IF x < 0.0 THEN x ¬ -x; IF x < 8.0 THEN { <> x2: REAL = x*x; P: REAL = (((( (-0.1849052456E3)*x2+ (0.7739233017E5))*x2+ (-0.1121442418E8))*x2+ (0.6516196407E9))*x2+ (-0.1336259035E11))*x2+ (0.5756849057E11); Q: REAL = (((( x2+ (0.2678532712E3))*x2+ (0.5927264853E5))*x2+ (0.9494680718E7))*x2+ (0.1029532985E10))*x2+ (0.5756849041E11); f ¬ P/Q; } ELSE { X0: REAL = x-piFourths; p0, q0: REAL; [p0, q0] ¬ PQ0[x]; f ¬ SqRt[2.0/(PI*x)]*(p0*Cos[radians: X0]-q0*Sin[radians: X0]); }; }; <<>> PQ0: PROC [x: REAL] RETURNS [P0, Q0: REAL] = { z: REAL = 8.0/x; z2: REAL = z*z; P0 ¬ (( -- PZERO 6502 - 8.78 digits (-0.165711693E-5)*z2+ (0.270871648E-4))*z2+ (-0.1098577711E-2))*z2+ (0.9999999983); Q0 ¬ ((( -- QZERO 6901 - 9.13 digits (0.576476505E-6)*z2+ (-0.6796319617E-5))*z2+ (0.1430262718E-3))*z2+ (-0.1562499927E-1))*z; }; pi34: REAL ¬ 3*piFourths; J1: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = { <> sign: REAL ¬ 1.0; IF x < 0.0 THEN {sign ¬ -1.0; x ¬ -x}; IF x < 8.0 THEN { <> <<..actually, JONE 6038 - 8.19 digits would be better, but that page is missing from my copy.>> x2: REAL = x*x; P: REAL = ((( (0.6543001487E1)*x2+ (-0.1988185996E4))*x2+ (0.1978615055E6))*x2+ (-0.7064812917E7))*x2+ (0.6706259146E8); Q: REAL = ((( x2+ (0.1766217266E3))*x2+ (0.2668508067E5))*x2+ (0.2635953855E7))*x2+ (0.1341252501E9); f ¬ sign*x*P/Q; } ELSE { X1: REAL = x-pi34; p1, q1: REAL; [p1, q1] ¬ PQ1[x]; f ¬ sign*SqRt[2.0/(PI*x)]*(p1*Cos[radians: X1]-q1*Sin[radians: X1]); }; }; <<>> PQ1: PROC [x: REAL] RETURNS [P1, Q1: REAL] = { z: REAL = 8.0/x; z2: REAL = z*z; P1 ¬ (( -- PONE 6702 - 8.72 digits (0.19796752E-5)*z2+ (-0.348677998E-4))*z2+ (0.1830991520E-2))*z2+ (0.1000000001E1); Q1 ¬ ((( -- QONE 7102 - 9.08 digits (-0.67222066E-6)*z2+ (0.831923005E-5))*z2+ (-0.2002434943E-3))*z2+ (0.4687499917E-1))*z; }; <<>> OldJn: PUBLIC PROC [n: INT, x: REAL] RETURNS [f: REAL] = { <> sign: REAL ¬ 1.0; IF n<0 THEN {n ¬ -n; sign ¬ IF n MOD 2 # 0 THEN -1.0 ELSE 1.0}; SELECT n FROM 0 => f ¬ J0[x]; 1 => f ¬ sign*J1[x]; ENDCASE => { v0: REAL ¬ J0[x]; v1: REAL ¬ J1[x]; FOR m: INT IN [2..n] DO v2: REAL ¬ 2*(m-1)*v1/x-v0; v0 ¬ v1; v1 ¬ v2; ENDLOOP; f ¬ sign*v1; }; }; lagestNumber: REAL ¬ Real.LargestNumber; factLimit: NAT ¬ 8; Jn: PUBLIC PROC [n: INT, x: REAL] RETURNS [f: REAL] = { <> sign: REAL ¬ 1.0; IF n<0 THEN {n ¬ -n; sign ¬ IF n MOD 2 # 0 THEN -1.0 ELSE 1.0}; SELECT TRUE FROM n=0 => RETURN [J0[x]]; n=1 => f ¬ sign*J1[x]; x=0.0 => RETURN [0.0]; ENDCASE => { <=0] (-1)k(x/2)n+2k / k! (n+k)!>> halfX: REAL ~ x*0.5; halfX2: REAL ~ halfX * halfX; term, sum: REAL; IF n <= factLimit THEN { term ¬ halfX; FOR j: NAT IN [2..n] DO term ¬ term * (halfX / j) ENDLOOP; } ELSE term ¬ Exp[n*Ln[halfX] - LnGamma[n+1] ]; sum ¬ term; FOR k: REAL ¬ 1, k+1 DO newSum: REAL; term ¬ -term * (halfX2 / (k*(n+k))); newSum ¬ sum + term; IF newSum=sum THEN EXIT ELSE sum ¬ newSum; ENDLOOP; f ¬ sum*sign; }; }; RationalFromReal: PUBLIC PROC [x: REAL, limit: INT ¬ LAST[INT]] RETURNS [RealFns.Rational] ~ { <<-- All best rational approximations of a real x may be obtained from x's continued fraction representation, x = c0 + 1/(c1 + 1/(c2 + 1/(...))) = [c0, c1, c2, ... ], by truncation to k terms and possibly "interpolation" of the last term. The continued fraction expansion itself is obtained by a variant of the standard GCD algorithm, which is folded into the recursions generating successive numerators and denominators. These recursions both have the same form: fk _ ck*fk-1 + fk-2. For further information, see CSL-Notebook entry 88CSLN-0045, Nov. 28, 1988; cited there is Fundamentals of Mathematics, Volume I, MIT Press, 1983. On return, |numer| <= limit, 0 < denom <= limit.>> tooLargeToFix: REAL ~ Real.FScale[1.0, BITS[INT] - 1]; -- 2147483648.0 tooSmallToFix: REAL ~ Real.FScale[1.0, 1 - BITS[INT]]; -- 4.65661280e-10 halfTooSmallToFix: REAL ~ Real.FScale[tooSmallToFix, -1]; -- 4.65661280e-10 / 2 expForInt: INT ~ 150; -- This exponent in SingleReal makes the significand an INT ratZero: RealFns.Rational ~ [0, 1]; sign: INT ¬ 1; flip: BOOL ¬ FALSE; -- If TRUE, nk and dk are swapped scale: INT; -- Power of 2 to get x into integer domain ak2, ak1, ak: CARD; -- GCD arguments, initially 1 and x ck, climit: CARD; -- ck is GCD quotient and c.f. term k nk, dk: INT; -- Result num. and den., recursively found nk1, dk2: INT ¬ 0; -- History terms for recursion nk2, dk1: INT ¬ 1; IF limit <= 0 THEN RETURN[ratZero]; -- Insist limit > 0 IF x < 0.0 THEN { x ¬ -x; sign ¬ -1; }; <<-- Handle first non-zero term of continued fraction, rest setup for integer GCD, sure to fit.>> IF x >= 1.0 THEN { -- First continued fraction term is non-zero rest: REAL; IF x >= tooLargeToFix OR (ck ¬ Real.Fix[x]) >= CARD[limit] THEN RETURN [[sign*limit, 1]]; flip ¬ TRUE; -- Keep denominator larger, for fast loop test nk ¬ 1; dk ¬ ck; -- Make new numerator and denominator rest ¬ x - ck; scale ¬ expForInt - LOOPHOLE[1.0, SingleReal].exp; ak ¬ Real.Fix[Real.FScale[rest, scale]]; ak1 ¬ Real.Fix[Real.FScale[1.0, scale]]; } ELSE { -- First continued fraction term is zero n, num: CARD; IF x <= tooSmallToFix -- Is x too tiny to be 1/INT ? THEN { IF x <= halfTooSmallToFix THEN RETURN[ratZero]; IF CARD[limit] > CARD[Real.Fix[0.5/x]] THEN RETURN [[sign, limit]] ELSE RETURN[ratZero]; }; <<-- Treating 1.0 and x as integers, divide 1/x in a peculiar way to get accurate remainder>> scale ¬ expForInt - LOOPHOLE[x, SingleReal].exp; ak1 ¬ Real.Fix[Real.FScale[x, scale]]; n ¬ MIN[BITS[CARD]-1, scale]; -- Stay within CARD arithmetic <> num ¬ Basics.BITLSHIFT[1, n]; ck ¬ num/ak1; -- First attempt at 1/x ak ¬ num MOD ak1; -- First attempt at remainder WHILE (scale ¬ scale-n) > 0 DO -- Shift quotient and remainder until done n ¬ MIN[8, scale]; -- That 8 is from 24 bits of x in a 32 bit CARD <> <> num ¬ Basics.BITLSHIFT[ak, n]; ck ¬ Basics.BITLSHIFT[ck, n] + num/ak1; ak ¬ num MOD ak1; -- Reduce remainder ENDLOOP; <> IF ck >= CARD[limit] -- Is x too tiny to be 1/limit ? THEN { IF CARD[2*limit] > ck THEN RETURN [[sign, limit]] ELSE RETURN[ratZero]; }; nk ¬ 1; dk ¬ ck; -- Make new numerator and denominator }; WHILE ak # 0 DO -- If possible, quit when have exact result ak2 ¬ ak1; ak1 ¬ ak; -- Prepare for next term nk2 ¬ nk1; nk1 ¬ nk; -- (This loop does essentially all the work) dk2 ¬ dk1; dk1 ¬ dk; ck ¬ ak2/ak1; -- Get next term of continued fraction ak ¬ ak2 - ck*ak1; -- Get remainder (GCD step) climit ¬ (limit - dk2)/dk1; -- Anticipate result of recursion on denominator IF climit <= ck THEN GOTO Interpolate; -- Do not let denominator exceed limit nk ¬ ck*nk1 + nk2; -- Make new result numerator and denominator dk ¬ ck*dk1 + dk2; REPEAT Interpolate => { twoClimit: CARD ¬ 2*climit; IF twoClimit >= ck THEN { -- If climit < ck/2 no improvement possible nk ¬ climit*nk1 + nk2; -- Make limited numerator and denominator dk ¬ climit*dk1 + dk2; IF twoClimit = ck THEN TRUSTED {-- If climit = ck improvement not sure -- Using climit is better only when dk2/dk1 > ak/ak1 dk2ak1Hi, dk2ak1Lo, dk1akHi, dk1akLo: CARD32; <<[dk2ak1Hi, dk2ak1Lo] _ Mul32[[lc[IF flip THEN nk2 ELSE dk2]],[lc[ak1]]];>> <<[dk1akHi, dk1akLo] _ Mul32[[lc[IF flip THEN nk1 ELSE dk1]],[lc[ak]]];>> [dk2ak1Hi, dk2ak1Lo] ¬ Mul32[IF flip THEN nk2 ELSE dk2, ak1]; [dk1akHi, dk1akLo] ¬ Mul32[IF flip THEN nk1 ELSE dk1, ak]; IF dk2ak1Hi < dk1akHi OR (dk2ak1Hi = dk1akHi AND dk2ak1Lo <= dk1akLo) THEN { nk ¬ nk1; dk ¬ dk1; }; -- Not an improvement, so undo step } }; }; ENDLOOP; IF flip THEN RETURN [[sign*dk, nk]] ELSE RETURN [[sign*nk, dk]]; }; OldRationalFromReal: PROC [x: REAL, limit: INT] RETURNS [RealFns.Rational] = { <> oneOverX: REAL; sign: INT ¬ 1; inverse: BOOL ¬ FALSE; p0, q0, lambda: INT; q1, p2: INT ¬ 0; p1, q2: INT ¬ 1; IF x < 0.0 THEN { x ¬ -x; sign ¬ -1 }; IF x = 1.0 THEN RETURN [[numerator: sign, denominator: 1]]; IF x > 1.0 THEN { x ¬ 1.0/x; inverse ¬ TRUE }; <> lambda ¬ 0; -- = FLOOR[x] WHILE (oneOverX ¬ (x-lambda))*REAL[limit-q1] > REAL[q2] DO <> <> < p2.>> x ¬ 1.0/oneOverX; lambda ¬ RealSupport.Fix[x]; -- was InlineFix q0 ¬ q1; q1 ¬ q2; q2 ¬ lambda*q1+q0; p0 ¬ p1; p1 ¬ p2; p2 ¬ lambda*p1+p0; ENDLOOP; IF NOT inverse THEN RETURN [[numerator: sign*p2, denominator: q2]] ELSE RETURN [[numerator: sign*q2, denominator: p2]]; }; AlmostZero: PUBLIC PROC [x: REAL, distance: INTEGER [-126..127]] RETURNS [BOOL] = { RETURN[LOOPHOLE[x, SingleReal].exp < (distance + ExponentBias)]; <> <<1. make zero and all denormalized values of x return TRUE>> <<2. make all NaN values of x return FALSE>> }; AlmostEqual: PUBLIC PROC [y, x: REAL, distance: INTEGER [-126..0]] RETURNS [BOOL] = { <> <> <> <> < TRUE>> < FALSE>> < TRUE>> <> xe: SingleReal ¬ LOOPHOLE[x]; ye: SingleReal ¬ LOOPHOLE[y]; SELECT TRUE FROM xe.exp = NaNExponent => RETURN [FALSE]; < FALSE>> ye.exp = NaNExponent => RETURN [FALSE]; < FALSE>> x = y => RETURN [TRUE]; <> <<(use REAL equality to check with -0.0)>> xe.sign # ye.sign => RETURN [FALSE]; < ABS[x-y] > MAX[ABS[x],ABS[y]]>> distance = 0 => RETURN [TRUE]; < ABS[x-y] < MAX[ABS[x],ABS[y]]>> xe.exp > ye.exp+1, ye.exp > xe.exp+1 => RETURN [FALSE]; < MAX[ABS[x],ABS[y]]*0.5 >> ENDCASE => { IF xe.exp = 0 OR ye.exp = 0 THEN { <> reNorm: REAL ¬ 1.0; LOOPHOLE[reNorm, SingleReal].exp ¬ LOOPHOLE[reNorm, SingleReal].exp + 24; xe ¬ LOOPHOLE[x ¬ x * reNorm]; ye ¬ LOOPHOLE[y ¬ y * reNorm]; }; { delta: SingleReal ¬ LOOPHOLE[(x - y)]; keyExp: INTEGER ¬ MAX[xe.exp, ye.exp] + distance; SELECT INTEGER[delta.exp] FROM < keyExp => RETURN [TRUE]; > keyExp => RETURN [FALSE]; ENDCASE => RETURN [delta.m = 0]; }; }; }; <> <> 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 ¬ twoPI-PIc1; -- 2*PI - c1 GetWithinTwoPI: PROC [radians: REAL] RETURNS [REAL] = { <> <> <= 0 >> xn: REAL ¬ RealSupport.Fix[radians * rec2pi]; -- was InlineFix RETURN [(radians - xn*PIc1) - xn*PIc2]; }; GetWithin360: PROC [degrees: REAL] RETURNS [REAL] = { <> <> <= 0 >> xn: REAL ¬ RealSupport.Fix[degrees * rec360]; -- was InlineFix RETURN [degrees - xn*360]; }; Mul32: PROC [x, y: CARD32] RETURNS [CARD32, CARD32] ~ { a: Basics.LongNumber ¬ [card[x]]; b: Basics.LongNumber ¬ [card[y]]; hi, lo, t1, t2: Basics.LongNumber; cy: CARDINAL; lo.card ¬ CARD32[a.lo]*CARD32[b.lo]; hi.card ¬ CARD32[a.hi]*CARD32[b.hi]; t1.card ¬ CARD32[a.hi]*CARD32[b.lo]; t2.card ¬ CARD32[a.lo]*CARD32[b.hi]; [cy, lo.hi] ¬ ADC3[lo.hi, t1.lo, t2.lo]; hi.card ¬ hi.card + t1.hi + t2.hi + cy; RETURN[hi.card, lo.card]; }; <> ADC3: PROC [a, b, c: CARD16] RETURNS [CARD16, CARD16] = INLINE { s: Basics.LongNumber; s.card ¬ CARD32[a] + CARD32[b] + CARD32[c]; RETURN[s.hi, s.lo]; }; 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 <> <> <> <<>>