DIRECTORY RealFns, RealFnsExtras; RealFnsExtrasImpl: CEDAR PROGRAM IMPORTS RealFns EXPORTS RealFnsExtras = BEGIN SinH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN ex: REAL = RealFns.Exp[x]; emx: REAL = 1/ex; f _ (ex-emx)/2.0 END; CosH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN ex: REAL = RealFns.Exp[x]; emx: REAL = 1/ex; f _ (ex+emx)/2.0 END; TanH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN IF ABS[x] > 20. THEN f _ (IF x<0. THEN -1. ELSE 1.) ELSE BEGIN ex: REAL = RealFns.Exp[x]; emx: REAL = 1/ex; f _ (ex-emx)/(ex+emx); END; END; CotH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN IF ABS[x] > 20. THEN f _ (IF x<0. THEN -1. ELSE 1.) ELSE BEGIN ex: REAL = RealFns.Exp[x]; emx: REAL = 1/ex; f _ (ex+emx)/(ex-emx); END; END; InvSinH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN f _ RealFns.Ln[x+RealFns.SqRt[x*x+1.0]]; END; InvCosH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN IF x < 1.0 THEN ERROR BadArgument; f _ RealFns.Ln[x+RealFns.SqRt[x*x-1.0]]; END; InvTanH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN IF ABS[x] >= 1.0 THEN ERROR BadArgument; f _ 0.5*RealFns.Ln[(1.0+x)/(1.0-x)]; END; InvCotH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN IF ABS[x] <= 1.0 THEN ERROR BadArgument; f _ 0.5*RealFns.Ln[(x+1.0)/(x-1.0)]; END; pi: REAL _ RealFnsExtras.pi; pi4: REAL _ 0.7853981633; lnRoot2pi: REAL _ 0.9189385332; LnGamma: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN phi: REAL; lnAdjFact: REAL _ 0; IF x <= 0.0 THEN ERROR BadArgument; IF x < 8 THEN BEGIN -- 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 _ RealFns.Ln[adjFact]; END; phi _ (-0.2762472E-2/(x*x)+0.8333327385E-1)/x; -- LGAM 5400 - 8.82 digits f _ (x-0.5)*RealFns.Ln[x]-x+lnRoot2pi+phi-lnAdjFact; -- Stirling's approx. END; Gamma: PUBLIC PROC [x: REAL] RETURNS [REAL] = {RETURN[RealFns.Exp[LnGamma[x]]]}; Erf: PUBLIC PROC [x: REAL] RETURNS [REAL] = {RETURN[1.0-Erfc[x]]}; Erfc: PUBLIC PROC [x: REAL] RETURNS [REAL] = {ERROR Unimplemented}; J0: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN IF x < 0.0 THEN x _ -x; IF x < 8.0 THEN BEGIN 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; END ELSE BEGIN X0: REAL = x-pi4; p0, q0: REAL; [p0, q0] _ PQ0[x]; f _ RealFns.SqRt[2.0/(pi*x)]*(p0*RealFns.Cos[radians: X0]-q0*RealFns.Sin[radians: X0]); END; END; PQ0: PROC [x: REAL] RETURNS [P0, Q0: REAL] = BEGIN 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; END; pi34: REAL _ 3*pi4; J1: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN sign: REAL _ 1.0; IF x < 0.0 THEN {sign _ -1.0; x _ -x}; IF x < 8.0 THEN BEGIN 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; END ELSE BEGIN X1: REAL = x-pi34; p1, q1: REAL; [p1, q1] _ PQ1[x]; f _ sign*RealFns.SqRt[2.0/(pi*x)]*(p1*RealFns.Cos[radians: X1]-q1*RealFns.Sin[radians: X1]); END; END; PQ1: PROC [x: REAL] RETURNS [P1, Q1: REAL] = BEGIN 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; END; Jn: PUBLIC PROC [n: INT, x: REAL] RETURNS [f: REAL] = BEGIN 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 => BEGIN 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; END; END; Y0: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = {ERROR Unimplemented}; Y1: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = {ERROR Unimplemented}; Yn: PUBLIC PROC [n: INT, x: REAL] RETURNS [f: REAL] = {ERROR Unimplemented}; K: PUBLIC PROC [m: REAL] RETURNS [REAL] = {ERROR Unimplemented}; E: PUBLIC PROC [m: REAL] RETURNS [REAL] = {ERROR Unimplemented}; Unimplemented: PUBLIC ERROR = CODE; BadArgument: PUBLIC ERROR = CODE; END. "RealFnsExtrasImpl.mesa Last Update by E. McCreight, November 16, 1984 11:48:22 am PST 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. Bessel function of first kind of order 0 JZERO 5838 - 8.54 digits Bessel function of first kind of order 1 JONE 6036 - 7.08 digits ..actually, JONE 6038 - 8.19 digits would be better, but that page is missing from my copy. Bessel function of first kind of order n Bessel function of second kind of order 0 Bessel function of second kind of order 1 Bessel function of second kind of order n Complete elliptic integral of the first kind Complete elliptic integral of the second kind Ê ÿ˜Jšœ™Jšœ>™>J˜šÏk ˜ Jšœ˜Jšœ˜J˜—š œœœœ œ˜HJš˜J˜š Ïnœœœœœœ˜/Jš˜Jšœœ˜Jšœœ˜Jšœ˜Jšœ˜J˜—š žœœœœœœ˜/Jš˜Jšœœ˜Jšœœ˜Jšœ˜Jšœ˜J˜—š žœœœœœœ˜/Jš˜Jš œœ œœœœ˜3š˜Jš˜Jšœœ˜Jšœœ˜Jšœ˜Jšœ˜—Jšœ˜J˜—š žœœœœœœ˜/Jš˜Jš œœ œœœœ˜3š˜Jš˜Jšœœ˜Jšœœ˜Jšœ˜Jšœ˜—Jšœ˜J˜—š žœœœœœœ˜2Jš˜Jšœ(˜(Jšœ˜J˜—š žœœœœœœ˜2Jš˜Jšœ œœ ˜"Jšœ(˜(Jšœ˜J˜—š žœœœœœœ˜2Jš˜Jšœœ œœ ˜(Jšœ$˜$Jšœ˜J˜—š žœœœœœœ˜2Jš˜Jšœœ œœ ˜(Jšœ$˜$Jšœ˜J˜—J˜JšœÜ™ÜJ˜Jšœœ˜Jšœœ˜Jšœ œ˜J˜š žœœœœœœ˜2Jš˜Jšœœ˜ Jšœ œ˜Jšœ œœ ˜#šœ˜ JšœÏc*˜0Jšœ œ˜šœ˜Jšœ˜Jšœ ˜ Jšœ˜—Jšœ ˜ Jšœ˜—Jšœ/Ÿ˜IJšœ5Ÿ˜JJšœ˜—J˜JšÐbnœœœœœœœ˜PJ˜Jšžœœœœœœœ˜BJ˜Jšžœœœœœœœ˜CJ˜J˜š žœœœœœœ˜-Jšœ(™(Jš˜Jšœ œ˜šœ ˜Jš˜Jšœ™Jšœœ˜šœœ˜Jšœ˜Jšœ˜Jšœ˜Jšœ˜Jšœ˜Jšœ˜—šœœ˜Jšœ˜Jšœ˜Jšœ˜Jšœ˜Jšœ˜Jšœ˜—Jšœœœ˜Jš˜—š˜Jš˜Jšœœ ˜Jšœœ˜ Jšœ˜JšœW˜WJšœ˜—Jšœ˜J™—š žœœœœ œ˜,Jš˜Jšœœ ˜Jšœœ˜šœŸ˜#Jšœ˜Jšœ˜Jšœ˜Jšœ˜—šœ Ÿ˜$Jšœ˜Jšœ˜Jšœ˜Jšœ˜—Jšœ˜J˜—Jšœœ ˜J˜š žœœœœœœ˜-Jšœ(™(Jš˜Jšœœ˜Jšœ œ˜&šœ ˜Jš˜šœ™J™[—Jšœœ˜šœœ˜ Jšœ˜Jšœ˜Jšœ˜Jšœ˜Jšœ˜—šœœ˜ Jšœ˜Jšœ˜Jšœ˜Jšœ˜Jšœ˜—Jšœ œœ˜Jš˜—š˜Jš˜Jšœœ ˜Jšœœ˜ Jšœ˜Jšœ\˜\Jšœ˜—Jšœ˜J™—š žœœœœ œ˜,Jš˜Jšœœ ˜Jšœœ˜šœŸ˜"Jšœ˜Jšœ˜Jšœ˜Jšœ˜—šœ Ÿ˜#Jšœ˜Jšœ˜Jšœ˜Jšœ˜—Jšœ˜J™—šžœœœœœœœ˜5Jšœ(™(Jš˜Jšœœ˜Jš œœœœœœ˜?šœ˜ J˜J˜šœ˜ Jš˜Jšœœ ˜Jšœœ ˜šœœœ˜Jšœœ˜Jšœ˜Jšœ˜Jšœ˜—J˜ Jšœ˜——Jšœ˜J˜J™—šžœœœœœœœ˜DJšœ)™)J™—šžœœœœœœœ˜DJšœ)™)J™—šžœœœœœœœœ˜LJšœ)™)J™J™—šœœœœœœœ˜@J™,—J˜šœœœœœœœ˜@J™-—J˜Jšœœœœ˜#Jšœ œœœ˜!J˜Jšœ˜—J˜—…—5