<> <> DIRECTORY RealFns, RealHypBessel; RealHypBesselImpl: CEDAR PROGRAM IMPORTS RealFns EXPORTS RealHypBessel = 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; <> pi: REAL _ 3.1415926536; pi4: REAL _ 0.7853981633; 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}; <> <<>> Y1: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = {ERROR}; <> <<>> Yn: PUBLIC PROC [n: INT, x: REAL] RETURNS [f: REAL] = {ERROR}; <> <<>> END.