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;
Bessel function 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.
pi: REAL ← 3.1415926536;
pi4: REAL ← 0.7853981633;
J0:
PUBLIC
PROC [x:
REAL]
RETURNS [f:
REAL] =
Bessel function of first kind of order 0
BEGIN
IF x < 0.0 THEN x ← -x;
IF x < 8.0
THEN
BEGIN
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;
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] =
Bessel function of first kind of order 1
BEGIN
sign: REAL ← 1.0;
IF x < 0.0 THEN {sign ← -1.0; x ← -x};
IF x < 8.0
THEN
BEGIN
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;
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] =
Bessel function of first kind of order n
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};
Bessel function of second kind of order 0
Y1:
PUBLIC
PROC [x:
REAL]
RETURNS [f:
REAL] = {
ERROR};
Bessel function of second kind of order 1
Yn:
PUBLIC
PROC [n:
INT, x:
REAL]
RETURNS [f:
REAL] = {
ERROR};
Bessel function of second kind of order n
END.