[_CDCSL_93-16_]<1>Cedar>release>FloatingPoint>RealFnsImpl.mesa

RealFnsImpl.mesa

Stewart on August 27, 1982 1:07 pm

Paul Rovner on May 4, 1983 10:13 am

Levin, August 8, 1983 4:39 pm

Plass, August 5, 1991 11:28 am PDT

E. McCreight, November 16, 1984 11:48:22 am PST (RealFnsExtrasImpl)

Russ Atkinson (RRA), July 5, 1989 11:36:56 pm PDT

Doug Wyatt, January 15, 1987 1:59:52 pm PST

Chauser, March 16, 1990 4:06 pm PST

Loner, February 26, 1991 8:33 am PST

Willie-s, August 6, 1991 1:06 pm PDT

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;

the bias, amount to subtract to make the exponent "correct"

DenormExponent: INTEGER = Exponent.FIRST;

biased, just as it appears in SingleReal

NaNExponent: INTEGER = Exponent.LAST;

biased, just as it appears in SingleReal

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)];

This declaration obviously depends on IEEE format 32-bit numbers (SingleReal). The idea is that the index field has the low bit of the exponent and the high 4 bits of the mantissa, which is used to index a table of initial guesses at the square root, each guess is roughly halfway between the lowest and highest root with the given index, which allows convergence with only two iterations. I don't know where the table originated.

Bomb: PROC [x: REAL, exCode: FloatingPointCommon.Exception ¬ invalidOperation] = {

ERROR FloatingPointCommon.Error[exCode];

};

Some useful constants (must be variables)

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] = {

compute e^x using continued fractions

ACM Algorithm 229, Morelock, James C., "Elementary functions by continued fractions"

r, f: REAL;

k: INTEGER;

fl: SingleReal ¬ LOOPHOLE[1.0, SingleReal];

To keep the error small near 0.0, use this Taylor's series

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 = 2k * exp(x) and 0<=x<Ln2

r ¬ x*x;

{

Expand the continued fraction loop a few times

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;

Note, if LoopCount > 3, then use the following comments

f ← (4*(LoopCount-3) - 2) + r/f;

f ← (4*(LoopCount-4) - 2) + r/f;

};

unscaled result is (f+x)/(f-x)

x ¬ (f + x)/(f - x);

scale: multipy by 2^k by adding k to the exponent

fl.exp ¬ k + ExponentBias;

RETURN[x*LOOPHOLE[fl, REAL]];

};

Log: PUBLIC PROC [base, arg: REAL] RETURNS [REAL] = {

logxy = Ln(y)/Ln(x)

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];

To keep the error small near 1.0, use this Taylor's series

SELECT ABS[xm1] FROM

> 3e-3 => {};

> 3e-4 => RETURN [xm1 - 0.5*xm1*xm1];

ENDCASE => RETURN [xm1];

if x = z * 2^m then Ln(x) = m * Ln2 + Ln(z). Here, 1/2 <= z < 1

fl ¬ LOOPHOLE[x, SingleReal];

m ¬ fl.exp - (ExponentBias - 1);

fl.exp ¬ 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 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] = {

RRA: The idea is to lookup a square root guess from the table, adjusting for the exponent, then perform a few iterations to get within a reasonable error. HIGHLY IEEE DEPENDENT!

SELECT TRUE FROM

a <= 0.0 => {

IF a < 0.0 THEN Bomb[a];

b ¬ 0.0;

};

LOOPHOLE[a, SingleReal].exp = 0 =>

Denormalized exponent. The guess is only valid when we have a normalized number, so we force it to be normalized, then correct for the shift.

b ¬ SqRt[a * (LONG[20000B] * LONG[20000B])] / 20000B;

ENDCASE => {

aFmt: SquareReal ¬ LOOPHOLE[a, SquareReal];

bFmt: SingleReal ¬ LOOPHOLE[sqGuesses[aFmt.index], SingleReal];

Get the initial guess based on the top 4 bits of the mantissa and the bottom bit of the exponent.

aExp: CARDINAL ¬ aFmt.expD2;

bFmt.exp ¬ bFmt.exp + aExp - 63;

Adjust the guess exponent based on adding half of the initial exponent (adjusted for the offset inherent in IEEE format).

b ¬ LOOPHOLE[bFmt, REAL];

The result is initially the guess

b ¬ LOOPHOLE[LOOPHOLE[b + a/b, LONG CARDINAL] - 40000000B, REAL];

This is really b ← (b + a/b)/2, but much faster

b ¬ LOOPHOLE[LOOPHOLE[b + a/b, LONG CARDINAL] - 40000000B, REAL];

This is really b ← (b + a/b)/2, but much faster

};

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] = {

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))

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];

We need to have rem IN [0..twoPI]

{

Now select the appropriate formula given the octant.

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];

We need to have rem IN [0..twoPI]

{

Now select the appropriate formula given the octant.

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];

We get better results if we compute in the degrees domain for a while.

sin ¬ Sin[degrees*radiansPerDegree];

};

CosDeg: PUBLIC PROC [degrees: REAL] RETURNS [cos: REAL] = {

degrees ¬ ABS[degrees];

IF degrees > 360 THEN degrees ¬ GetWithin360[degrees];

We get better results if we compute in the degrees domain for a while.

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 f is small, then tan(f) H f, otherwise use ratl approximation from Cody & Waite

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 + x3/6 + x5/40 + ... H x

< 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]];

RETURN[2*ArcTan[SqRt[(1-x)/(1+x)], 1]];

};

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]

The ArcCos Geometry for Method 2

Constants for ArcTan

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] = {

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[-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)];

};

The following approximations are taken from the book "Computer Approximations", by Hart, Cheney, Lawson, Maehly, Mesztenyi, Rice, Thacher, and Witzgall, published by Wiley. Function names and indexes refer to that book.

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] = {

Bessel function of first kind of order 0

IF x < 0.0 THEN x ¬ -x;

IF x < 8.0 THEN {

JZERO 5838 - 8.54 digits

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] = {

Bessel function of first kind of order 1

sign: REAL ¬ 1.0;

IF x < 0.0 THEN {sign ¬ -1.0; x ¬ -x};

IF x < 8.0 THEN {

JONE 6036 - 7.08 digits

..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] = {

Bessel function of first kind of order n

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] = {

Bessel function of first kind of order n

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 => {

Jn(x) = Sum[k>=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[[lc[1]], n].lc;

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[[li[ak]], n].lc;

ck ← Basics.BITLSHIFT[[lc[ck]], n].lc + num/ak1;

num ¬ Basics.BITLSHIFT[ak, n];

ck ¬ Basics.BITLSHIFT[ck, n] + num/ak1;

ak ¬ num MOD ak1; -- Reduce remainder

ENDLOOP;

All done with divide

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] = {

The best rational approximation to a real x is obtained using continued fractions. The following algorithm is adapted from Wm. J. LeVeque's book, Fundamentals of Number Theory, Addison Wesley, 1977, page 229. |num| <= limit, 0 < denom <= limit.

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 };

Here x lies in [0..1)

lambda ¬ 0; -- = FLOOR[x]

WHILE (oneOverX ¬ (x-lambda))*REAL[limit-q1] > REAL[q2] DO

This test ensures that q2 never exceeds limit.

Because the initial x was <1, p2<=q2.

The best approximation so far is p2/q2, and q2 > 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)];

The range of distance has been carefully chosen to

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] = {

AlmostEqual[y, x, distance] ==

NotNaN[x] AND NotNaN[y] AND

ABS[x-y] <= MAX[ABS[x], ABS[y]]*2^distance

for valid x & y:

x = y => TRUE

Sign[x] # Sign[y] => FALSE

Sign[x] = Sign[y] AND distance = 0 => TRUE

Rewritten by RRA on November 11, 1986, previous implementation returned FALSE when x = y = 0.0, which was counterintuitive, to say the least.

xe: SingleReal ¬ LOOPHOLE[x];

ye: SingleReal ¬ LOOPHOLE[y];

SELECT TRUE FROM

xe.exp = NaNExponent => RETURN [FALSE];

NotNaN[x] => FALSE

ye.exp = NaNExponent => RETURN [FALSE];

NotNaN[y] => FALSE

x = y => RETURN [TRUE];

exactly equal valid numbers are always AlmostEqual

(use REAL equality to check with -0.0)

xe.sign # ye.sign => RETURN [FALSE];

xe.sign # ye.sign => ABS[x-y] > MAX[ABS[x],ABS[y]]

distance = 0 => RETURN [TRUE];

xe.sign = ye.sign AND x # y => ABS[x-y] < MAX[ABS[x],ABS[y]]

xe.exp > ye.exp+1, ye.exp > xe.exp+1 => RETURN [FALSE];

ABS[x-y] > MAX[ABS[x],ABS[y]]*0.5

ENDCASE => {

IF xe.exp = 0 OR ye.exp = 0 THEN {

Sigh, one of them is denormalized!

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];

};

Helpful internal functions

Coefficients for the taylor series for sin and cos.

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];

roughly 5.96e-8, the minimum relative error between two different reals

sinCosEps: REAL ¬ SqRt[minEps];

RRA: IF radians <= sinEps THEN SinFirstOctant[radians] = radians AND CosFirstOctant[radians] = 1.0 to the precision of our number system. The contribution of the smaller term of two consecutive terms in the Sin or Cos polynomials is down by the square from the larger (and a factor of at least two). So sinCosEps roughly equals 1.726e-4. Testing for this case improves accuracy at small values.

SinFirstOctant: PROC [radians: REAL] RETURNS [REAL] = {

This version of Sin is only guaranteed to be accurate within the first octant: [0.0..PI/4], which is approximately [0.0..0.785+]. When radians = PI/4 we have the maximum contribution of the term after the fact9recip term, for a truncation error of roughly 1.76e-9, which is within minEps*Sin[radians] (= 4.2e-8).

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] = {

This version of Cos is only guaranteed to be accurate within the first octant: [0.0..PI/4]. When radians = PI/4 we have the maximum contribution of the terms after the fact8recip term, for a truncation error of roughly 2.46e-8, which is within minEps*Cos[PI/4] (= 4.2e-8).

IF radians <= sinCosEps THEN RETURN [1.0] ELSE {

r2: REAL ~ radians*radians;

RETURN [((((fact8recip*r2-fact6recip)*r2+fact4recip)*r2-fact2recip)*r2+1)];

};

Two constants used to add N*PI with more than normal precision.

PIc1: REAL ¬ 6.28125; -- 402.0/64 (exact representation near twoPI)

PIc2: REAL ¬ twoPI-PIc1; -- 2*PI - c1

GetWithinTwoPI: PROC [radians: REAL] RETURNS [REAL] = {

Let the returned value be R. Then

radians <= 0 ƒ radians H xn*twoPI + R ' -twoPI d R d 0

radians >= 0 ƒ radians H xn*twoPI + R ' 0 d R d twoPI

xn: REAL ¬ RealSupport.Fix[radians * rec2pi]; -- was InlineFix

RETURN [(radians - xn*PIc1) - xn*PIc2];

};

GetWithin360: PROC [degrees: REAL] RETURNS [REAL] = {

Let the returned value be R. Then

degrees <= 0 ƒ degrees H xn*360 + R ' -360 d R d 0

degrees >= 0 ƒ degrees H xn*360 + R ' 0 d R d 360

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];

};

Add with Carry

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

Carl Hauser, January 21, 1988 2:08:18 pm PST

For PortableCedar

Removed instances of Real.InlineFix* , replacing with RealSupport.Fix and labelling what was replaced.