DRealFnsImpl.mesa
Copyright Ó 1989, 1991 by Xerox Corporation. All rights reserved.
Russ Atkinson (RRA) June 22, 1989 10:14:23 pm PDT
Mna, February 18, 1991 4:14 pm PST
DIRECTORY
DRealFns,
DRealSupport,
FloatingPointCommon;
DRealFnsImpl: CEDAR PROGRAM
IMPORTS DRealSupport, FloatingPointCommon
EXPORTS DRealFns
= BEGIN OPEN DRealSupport;
Rational: TYPE = DRealFns.Rational;
Exponent and logarithm functions
Exp: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
For an input argument n, returns e^n (e=2.718...).
DRealExpI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealExpI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = exp(DRealPtr(x));\n";
" };\n";
".XR𡤍RealExpI";
};
DRealExpI[@ret, @x];
};
Ln: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
Computes the natural logarithm (base e) of the input argument.
DRealLogI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealLogI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = log(DRealPtr(x));\n";
" };\n";
".XR𡤍RealLogI";
};
DRealLogI[@ret, @x];
};
Log: PUBLIC PROC [base, arg: DREAL] RETURNS [DREAL] = {
Computes logarithm to the base base of arg.
RETURN [Ln[arg]/Ln[base]];
};
Power: PUBLIC PROC [base, exponent: DREAL] RETURNS [DREAL] = {
Calculates base to the exponent power by e(exponent*Ln(base)).
RETURN [Exp[exponent*Ln[base]]];
};
Root: PUBLIC PROC [index, arg: DREAL] RETURNS [DREAL] = {
Calculates the index root of arg by e(Ln(arg)/index).
RETURN [Exp[Ln[arg]/index]];
};
SqRt: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
DRealSqrtI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealSqrtI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = sqrt(DRealPtr(x));\n";
" };\n";
".XR𡤍RealSqrtI";
};
DRealSqrtI[@ret, @x];
};
Trigonometric functions
Sin: PUBLIC PROC [radians: DREAL] RETURNS [ret: DREAL] = TRUSTED {
DRealSinI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealSinI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = sin(DRealPtr(x));\n";
" };\n";
".XR𡤍RealSinI";
};
DRealSinI[@ret, @radians];
};
SinDeg: PUBLIC PROC [degrees: DREAL] RETURNS [DREAL] = {
RETURN [Sin[degreesToRadians*degrees]];
};
Cos: PUBLIC PROC [radians: DREAL] RETURNS [ret: DREAL] = TRUSTED {
DRealCosI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealCosI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = cos(DRealPtr(x));\n";
" };\n";
".XR𡤍RealCosI";
};
DRealCosI[@ret, @radians];
};
CosDeg: PUBLIC PROC [degrees: DREAL] RETURNS [DREAL] = {
RETURN [Cos[degreesToRadians*degrees]];
};
Tan: PUBLIC PROC [radians: DREAL] RETURNS [ret: DREAL] = TRUSTED {
DRealTanI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealTanI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = tan(DRealPtr(x));\n";
" };\n";
".XR𡤍RealTanI";
};
DRealTanI[@ret, @radians];
};
TanDeg: PUBLIC PROC [degrees: DREAL] RETURNS [DREAL] = {
RETURN [Tan[degreesToRadians*degrees]];
};
ArcTan: PUBLIC PROC [y, x: DREAL] RETURNS [radians: DREAL] = TRUSTED {
DRealArcTanI: UNSAFE PROC [ret, y, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealArcTanI (ret, y, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = atan2(DRealPtr(y), DRealPtr(x));\n";
" };\n";
".XR𡤍RealArcTanI";
};
DRealArcTanI[@radians, @y, @x];
};
ArcTanDeg: PUBLIC PROC [y, x: DREAL] RETURNS [degrees: DREAL] = {
radians: DREAL = ArcTan[y, x];
RETURN [radiansToDegrees*radians];
};
Hyperbolic circle functions
SinH: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
Hyperbolic sine
DRealSinHI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealSinHI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = sinh(DRealPtr(x));\n";
" };\n";
".XR𡤍RealSinHI";
};
DRealSinHI[@ret, @x];
};
CosH: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
Hyperbolic cosine
DRealCosHI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealCosHI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = cosh(DRealPtr(x));\n";
" };\n";
".XR𡤍RealCosHI";
};
DRealCosHI[@ret, @x];
};
TanH: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
Hyperbolic tangent
DRealTanHI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealTanHI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = tanh(DRealPtr(x));\n";
" };\n";
".XR𡤍RealTanHI";
};
DRealTanHI[@ret, @x];
};
CotH: PUBLIC PROC [x: DREAL] RETURNS [DREAL] = {
Hyperbolic cotangent
RETURN [1.0/TanH[x]];
};
InvSinH: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
DRealInvSinHI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealInvSinHI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = asinh(DRealPtr(x));\n";
" };\n";
".XR𡤍RealInvSinHI";
};
DRealInvSinHI[@ret, @x];
};
InvCosH: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
...defined for x >= 1.0, returns non-negative result
DRealInvCosHI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealInvCosHI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = acosh(DRealPtr(x));\n";
" };\n";
".XR𡤍RealInvCosHI";
};
DRealInvCosHI[@ret, @x];
};
InvTanH: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
...defined for x IN (-1.0..1.0)
DRealInvTanHI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealInvTanHI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = atanh(DRealPtr(x));\n";
" };\n";
".XR𡤍RealInvTanHI";
};
DRealInvTanHI[@ret, @x];
};
InvCotH: PUBLIC PROC [x: DREAL] RETURNS [DREAL] = {
...defined for x NOT IN [-1.0..1.0]
RETURN [1.0/InvTanH[x]];
};
Gamma functions
LnGamma: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
.. defined for x>0
DRealLnGammaI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealLnGammaI (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = lgamma(DRealPtr(x));\n";
" };\n";
".XR𡤍RealLnGammaI";
};
DRealLnGammaI[@ret, @x];
};
Gamma: PUBLIC PROC [x: DREAL] RETURNS [DREAL] = {
.. defined for x>0; computed by Exp[LnGamma[x]]
Gamma[x+1] = x! for integer x.
RETURN [Exp[LnGamma[x]]];
};
Bessel functions
J0: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
Bessel function of first kind of order 0
DRealJ0I: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealJ0I (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = j0(DRealPtr(x));\n";
" };\n";
".XR𡤍RealJ0I";
};
DRealJ0I[@ret, @x];
};
J1: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
Bessel function of first kind of order 1
DRealJ1I: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealJ1I (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = j1(DRealPtr(x));\n";
" };\n";
".XR𡤍RealJ1I";
};
DRealJ1I[@ret, @x];
};
Jn: PUBLIC PROC [n: INT, x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
Bessel function of first kind of order n
DRealJnI: UNSAFE PROC [ret: PDREAL, n: INTEGER, x: PDREAL]
= UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealJnI (ret, n, x) W2 *ret; word n; W2 *x; {\n";
" DRealPtr(ret) = jn(((int) n), DRealPtr(x));\n";
" };\n";
".XR𡤍RealJnI";
};
DRealJnI[@ret, n, @x];
};
Y0: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
Bessel function of second kind of order 0
DRealY0I: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealY0I (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = y0(DRealPtr(x));\n";
" };\n";
".XR𡤍RealY0I";
};
DRealY0I[@ret, @x];
};
Y1: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
Bessel function of second kind of order 1
DRealY1I: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealY1I (ret, x) W2 *ret, *x; {\n";
" DRealPtr(ret) = y1(DRealPtr(x));\n";
" };\n";
".XR𡤍RealY1I";
};
DRealY1I[@ret, @x];
};
Yn: PUBLIC PROC [n: INT, x: DREAL] RETURNS [ret: DREAL] = TRUSTED {
Bessel function of second kind of order n
DRealYnI: UNSAFE PROC [ret: PDREAL, n: INTEGER, x: PDREAL]
= UNCHECKED MACHINE CODE {
"+extern void XR𡤍RealYnI (ret, n, x) W2 *ret; word n; W2 *x; {\n";
" DRealPtr(ret) = yn(((int) n), DRealPtr(x));\n";
" };\n";
".XR𡤍RealYnI";
};
DRealYnI[@ret, n, @x];
};
Rational approximation
Rational: TYPE ~ RECORD [numerator: INT, denominator: INT];
RationalFromDReal: PUBLIC PROC [x: DREAL, limit: INT ¬ INT.LAST] RETURNS [Rational] = {
... computes the best rational approximation to a real using continued fractions.
|numerator| <= limit, 0 < denominator <= limit.
sign: INT ¬ 1;
inverse: BOOL ¬ FALSE;
p0, q0: INT ¬ 0;
lambda: INT ¬ 0;
q1, p2: INT ¬ 0;
p1, q2: INT ¬ 1;
IF limit < 1 THEN ERROR FloatingPointCommon.Error[other];
IF x > limit THEN RETURN [[limit, 1]];
IF -x > limit THEN RETURN [[-limit, 1]];
IF x < 0.0 THEN { x ¬ -x; sign ¬ -1 };
IF x = 1.0 THEN RETURN [[numerator: sign, denominator: 1]];
IF x > 1.0 THEN { x ¬ 1.0/x; inverse ¬ TRUE };
Here x lies in [0..1), lambda = Fix[x]
DO
oneOverX: DREAL ¬ (x-lambda);
IF oneOverX*(limit-q1) <= q2 THEN EXIT;
This test ensures that q2 never exceeds limit.
Because the initial x was <1, p2<=q2.
The best approximation so far is p2/q2, and q2 > p2.
x ¬ 1.0/oneOverX;
lambda ¬ DRealSupport.Fix[x];
q0 ¬ q1;
q1 ¬ q2;
q2 ¬ lambda*q1+q0;
p0 ¬ p1;
p1 ¬ p2;
p2 ¬ lambda*p1+p0;
ENDLOOP;
IF NOT inverse
THEN RETURN [[numerator: sign*p2, denominator: q2]]
ELSE RETURN [[numerator: sign*q2, denominator: p2]];
};
Approximate equality
AlmostZero: PUBLIC PROC [x: DREAL, distance: INTEGER] RETURNS [BOOL] = {
AlmostZero[x, distance] == ABS[x] < 2^distance
IF x = 0.0 THEN RETURN [TRUE];
RETURN [ABS[x] < DRealSupport.FScale[1.0, distance]];
};
AlmostEqual: PUBLIC PROC [y, x: DREAL, distance: INTEGER] RETURNS [BOOL] = {
AlmostEqual[y, x, distance] ==
NotNaN[x] AND NotNaN[y] AND
ABS[x-y] <= MAX[ABS[x], ABS[y]] * 2^distance
Notes:
When either number is a Nan (not-a-number), then equality is never guaranteed, so FALSE is a conservative answer.
Comparing anything against 0.0 can give surprising results using the above definition, since ABS[x-0] = MAX[ABS[x], ABS[0]], which implies that AlmostEqual[0, x, distance] is only true for distance = 0.
IF DRealSupport.Classify[x] > normal THEN RETURN [FALSE];
IF DRealSupport.Classify[y] > normal THEN RETURN [FALSE];
{
delta: DREAL = ABS[x-y];
absX: DREAL = ABS[x];
absY: DREAL = ABS[y];
max: DREAL = IF absX < absY THEN absY ELSE absX;
IF delta = max THEN RETURN [TRUE];
RETURN [delta <= DRealSupport.FScale[max, distance]];
};
};
Initialization
PREAL: TYPE = POINTER TO REAL;
PDREAL: TYPE = POINTER TO DREAL;
pi: DREAL = 3.14159265358979323846264338327950288419716939937510;
degreesToRadians: DREAL = pi/180;
radiansToDegrees: DREAL = 180/pi;
DefInclude: PROC = TRUSTED MACHINE CODE {
"*#include <math.h>\n";
"."
};
DefDRealPtr: PROC = TRUSTED MACHINE CODE {
"+#define DRealPtr(x) (*((double *) (x)))\n."
};
DefInclude[];
DefDRealPtr[];
END.
Russ Atkinson (RRA) May 29, 1989 5:05:29 pm PDT
Adapted from RealFns