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;
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.
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] =
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 Unimplemented};
Bessel function of second kind of order 0
Y1:
PUBLIC
PROC [x:
REAL]
RETURNS [f:
REAL] = {
ERROR Unimplemented};
Bessel function of second kind of order 1
Yn:
PUBLIC
PROC [n:
INT, x:
REAL]
RETURNS [f:
REAL] = {
ERROR Unimplemented};
Bessel function of second kind of order n
K:
PUBLIC
PROC [m:
REAL]
RETURNS [
REAL] = {
ERROR Unimplemented};
Complete elliptic integral of the first kind
E:
PUBLIC
PROC [m:
REAL]
RETURNS [
REAL] = {
ERROR Unimplemented};
Complete elliptic integral of the second kind
Unimplemented: PUBLIC ERROR = CODE;
BadArgument: PUBLIC ERROR = CODE;
END.