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 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 LDragonRealFnsImpl.mesa Last modified: Stewart on August 27, 1982 1:07 pm Last modified: Paul Rovner on May 4, 1983 10:13 am Last Edited by: Levin, August 8, 1983 4:39 pm compute e^n using continued fractions See ACM Algorithm 229, Morelock, James C., "Elementary functions by continued fractions" f_4*LoopCount+2; unscaled result is (f+x)/(f-x) scale: multipy by 2^exp by adding exp to the exponent logxy = Ln(y)/Ln(x) if x = z * 2^m then Ln(x) = m * Ln2 + Ln(z). Here, 1/2 <= z < 1 special case for limit condition (x a power of 2) yth root of x = e^(Ln(x)/y) x^y = e^(y*Ln(x)) This function is good to 7.33 decimal places and is taken from: Computer Approximations, John F. Hart et.al. pp118(top 2nd col). This function is good to 8.7 decimal places and is taken from: Computer Approximations, John F. Hart et.al. pp129(top 2nd col). Ê †˜Jšœ™Jšœ1™1Jšœ2™2Jšœ-™-J˜šÏk ˜ Jšœ œ0˜BJšœœ+˜5J˜J˜—š œ œœœœ ˜GJš˜J˜Jšœ œ˜Jšœ œ˜J˜Jšœœ˜Jšœ œ ˜J˜Jšœ%™%Jšœ@™@Jšœ™š Ïnœœœœœœ˜-Jšœœ˜Jšœœœ˜Jšœ œ˜Jšœœœ˜J˜JšœÏc˜1JšœŸ˜Jšœ Ÿ˜J˜Jšœ Ÿ(˜3J˜J˜Jšœ™J˜J˜J˜šœ œœ˜%Jšœ7œ˜?J˜—Jšœ™JšœŸ3˜HJ˜Jšœ˜ J˜J˜—šžœœœœœœœ˜™>Jšœ@™@Jšœœ˜šœœœ˜"Jš œœœœœ˜,Jš œœœœ œœ ˜7—Jšœœœœœ œœ ˜FJšœœœœ˜;Jšœœ˜ šœ˜ Jšœ"˜$JšœD˜FJšœE˜GJšœF˜HJ˜(Jšœ˜—J˜ J˜9J˜J˜—š ž œœœœœ œ˜?Jšœ œ/˜