<<>> <> <> <> <> DIRECTORY DRealFns, DRealSupport, FloatingPointCommon; DRealFnsImpl: CEDAR PROGRAM IMPORTS DRealSupport, FloatingPointCommon EXPORTS DRealFns = BEGIN OPEN DRealSupport; Rational: TYPE = DRealFns.Rational; <> Exp: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealExpI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealExpI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = exp(DRealPtr(x));\n"; " };\n"; ".XR_DRealExpI"; }; DRealExpI[@ret, @x]; }; <<>> Ln: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealLogI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealLogI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = log(DRealPtr(x));\n"; " };\n"; ".XR_DRealLogI"; }; DRealLogI[@ret, @x]; }; <<>> Log: PUBLIC PROC [base, arg: DREAL] RETURNS [DREAL] = { <> RETURN [Ln[arg]/Ln[base]]; }; <<>> Power: PUBLIC PROC [base, exponent: DREAL] RETURNS [DREAL] = { <> RETURN [Exp[exponent*Ln[base]]]; }; <<>> Root: PUBLIC PROC [index, arg: DREAL] RETURNS [DREAL] = { <> 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_DRealSqrtI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = sqrt(DRealPtr(x));\n"; " };\n"; ".XR_DRealSqrtI"; }; DRealSqrtI[@ret, @x]; }; <<>> <> Sin: PUBLIC PROC [radians: DREAL] RETURNS [ret: DREAL] = TRUSTED { DRealSinI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealSinI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = sin(DRealPtr(x));\n"; " };\n"; ".XR_DRealSinI"; }; 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_DRealCosI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = cos(DRealPtr(x));\n"; " };\n"; ".XR_DRealCosI"; }; 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_DRealTanI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = tan(DRealPtr(x));\n"; " };\n"; ".XR_DRealTanI"; }; 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_DRealArcTanI (ret, y, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = atan2(DRealPtr(y), DRealPtr(x));\n"; " };\n"; ".XR_DRealArcTanI"; }; DRealArcTanI[@radians, @y, @x]; }; ArcTanDeg: PUBLIC PROC [y, x: DREAL] RETURNS [degrees: DREAL] = { radians: DREAL = ArcTan[y, x]; RETURN [radiansToDegrees*radians]; }; <> SinH: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealSinHI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealSinHI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = sinh(DRealPtr(x));\n"; " };\n"; ".XR_DRealSinHI"; }; DRealSinHI[@ret, @x]; }; CosH: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealCosHI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealCosHI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = cosh(DRealPtr(x));\n"; " };\n"; ".XR_DRealCosHI"; }; DRealCosHI[@ret, @x]; }; TanH: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealTanHI: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealTanHI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = tanh(DRealPtr(x));\n"; " };\n"; ".XR_DRealTanHI"; }; DRealTanHI[@ret, @x]; }; CotH: PUBLIC PROC [x: DREAL] RETURNS [DREAL] = { <> 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_DRealInvSinHI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = asinh(DRealPtr(x));\n"; " };\n"; ".XR_DRealInvSinHI"; }; 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_DRealInvCosHI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = acosh(DRealPtr(x));\n"; " };\n"; ".XR_DRealInvCosHI"; }; 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_DRealInvTanHI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = atanh(DRealPtr(x));\n"; " };\n"; ".XR_DRealInvTanHI"; }; DRealInvTanHI[@ret, @x]; }; <<>> InvCotH: PUBLIC PROC [x: DREAL] RETURNS [DREAL] = { <<...defined for x NOT IN [-1.0..1.0]>> RETURN [1.0/InvTanH[x]]; }; <<>> <> 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_DRealLnGammaI (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = lgamma(DRealPtr(x));\n"; " };\n"; ".XR_DRealLnGammaI"; }; 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]]]; }; <<>> <> J0: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealJ0I: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealJ0I (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = j0(DRealPtr(x));\n"; " };\n"; ".XR_DRealJ0I"; }; DRealJ0I[@ret, @x]; }; <<>> J1: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealJ1I: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealJ1I (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = j1(DRealPtr(x));\n"; " };\n"; ".XR_DRealJ1I"; }; DRealJ1I[@ret, @x]; }; <<>> Jn: PUBLIC PROC [n: INT, x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealJnI: UNSAFE PROC [ret: PDREAL, n: INTEGER, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealJnI (ret, n, x) W2 *ret; word n; W2 *x; {\n"; " DRealPtr(ret) = jn(((int) n), DRealPtr(x));\n"; " };\n"; ".XR_DRealJnI"; }; DRealJnI[@ret, n, @x]; }; Y0: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealY0I: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealY0I (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = y0(DRealPtr(x));\n"; " };\n"; ".XR_DRealY0I"; }; DRealY0I[@ret, @x]; }; <<>> Y1: PUBLIC PROC [x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealY1I: UNSAFE PROC [ret, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealY1I (ret, x) W2 *ret, *x; {\n"; " DRealPtr(ret) = y1(DRealPtr(x));\n"; " };\n"; ".XR_DRealY1I"; }; DRealY1I[@ret, @x]; }; <<>> Yn: PUBLIC PROC [n: INT, x: DREAL] RETURNS [ret: DREAL] = TRUSTED { <> DRealYnI: UNSAFE PROC [ret: PDREAL, n: INTEGER, x: PDREAL] = UNCHECKED MACHINE CODE { "+extern void XR_DRealYnI (ret, n, x) W2 *ret; word n; W2 *x; {\n"; " DRealPtr(ret) = yn(((int) n), DRealPtr(x));\n"; " };\n"; ".XR_DRealYnI"; }; DRealYnI[@ret, n, @x]; }; <<>> <> <> <<>> 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 }; <> DO oneOverX: DREAL ¬ (x-lambda); IF oneOverX*(limit-q1) <= q2 THEN EXIT; <> <> < 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]]; }; <<>> <> AlmostZero: PUBLIC PROC [x: DREAL, distance: INTEGER] RETURNS [BOOL] = { <> 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] = { <> <> <> <> <> <> 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]]; }; }; <<>> <> 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 \n"; "." }; DefDRealPtr: PROC = TRUSTED MACHINE CODE { "+#define DRealPtr(x) (*((double *) (x)))\n." }; DefInclude[]; DefDRealPtr[]; END. <> <> <<>>