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