## math64-intel32.pkg
## *************************************************************************
## *
## Copyright (c) 1985 Regents of the University of California. *
## *
## Use and reproduction of this software are granted in accordance with *
## the terms and conditions specified in the Berkeley Software License *
## Agreement (in particular, this entails acknowledgement of the programs' *
## source, and inclusion of this notice) with the additional understanding *
## that all recipients should regard themselves as participants in an *
## ongoing research project and hence should feel obligated to report *
## their experiences (good or bad) with these elementary function codes, *
## using "sendbug 4bsd-bugs@BERKELEY", to the authors. *
## *
## K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan. *
## Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85. *
## *
## *************************************************************************
# Compiled by:
# src/lib/std/src/standard-core.sublib# The following functions were adapted from the 4.3BSD math library.
# Eventually, each machine supported should have a hand-coded math
# generic with more efficient versions of these functions.
### "The world is full of magical things
### patiently waiting for our wits
### to grow sharper."
###
### -- Bertrand Russell
package math64: (weak) Math { # Math is from src/lib/std/src/math.api # inline_t is from src/lib/core/init/built-in.pkg Float = Float;
infix my 50 ==== ;
my (+) = inline_t::f64::(+);
my (-) = inline_t::f64::(-);
my (*) = inline_t::f64::(*);
my (/) = inline_t::f64::(/);
my (-_) = inline_t::f64::neg;
my neg = inline_t::f64::neg;
my (<) = inline_t::f64::(<);
my (<=) = inline_t::f64::(<=);
my (>) = inline_t::f64::(>);
my (>=) = inline_t::f64::(>=);
my (====) = inline_t::f64::(====);
package i = inline_t::default_int; # inline_t is from src/lib/core/init/built-in.pkg my lessu: (Int, Int) -> Bool = i::ltu;
pi = 3.14159265358979323846;
e = 2.7182818284590452354;
fun is_nan x = bool::not (x====x);
plus_infinity = 1E300 * 1E300;
minus_infinity = -plus_infinity;
na_n = 0.0 / 0.0;
two_to_the_54 = 18014398509481984.0;
two_to_the_minus_54 = 1.0 / 18014398509481984.0;
# This function is IEEE double-precision specific;
# it works correctly on subnormal inputs and outputs;
# we do not apply it to inf's and nan's:
#
fun scalb (x, k)
=
{ j = runtime::asm::logb x;
k' = (i::(+))(k, j);
if (j == -1023)
scalb (x*two_to_the_54, (i::(-))(k, 54)); # 2
elif (lessu((i::(+))(k', 1022), 2046))
runtime::asm::scalb (x, k); # 1
elif ((i::(<))(k', 0))
if ((i::(<))(k', (i::(-))(-1022, 54)))
0.0; # 3
else
scalb (x, (i::(+))(k, 54)) * two_to_the_minus_54; # 4
fi;
else
x * plus_infinity; # 5
fi;
};
# Proof of correctness of scalb: (Appel)
# 1. if x is normal and x*2^k is normal
# then case /*1*/ applies, computes right answer
# 2. if x is subnormal and x*2^k is normal
# then case /*2*/ reduces problem to case 1.
# 3. if x*2^k is sub-subnormal (i.e. 0)
# then case /*3*/ applies, returns 0.0
# 4. if x*2^k is subnormal
# then -1076 < k' <= -1023, case /*4*/ applies,
# computes right answer
# 5. if x*2^k is supernormal (i.e. infinity)
# then case /*5*/ computes right answer
# This function is IEEE double-precision specific;
# it works correctly on subnormal inputs;
# must not be applied to inf's and nan's
#
fun logb x
=
case (runtime::asm::logb x)
-1023 => (i::(-))(runtime::asm::logb (x*two_to_the_54), 54); # Denormalized number
i => i;
esac;
negone = -1.0;
zero = 0.0;
half = 0.5;
one = 1.0;
two = 2.0;
four = 4.0;
# * SHOULD BE INLINE OP *
/* may be applied to inf's and nan's
GETS MINUS-ZERO WRONG! XXX BUGGO FIXME
*/
fun copysign (a, b)
=
case(a < zero, b < zero)
(TRUE, TRUE ) => a;
(FALSE, FALSE ) => a;
_ => -a;
esac;
# May be applied to inf's and nan's
#
fun abs x
=
x < zero ?? -x
:: x;
fun (mod) (a, b)
=
(i::(-))(a, (i::(*))(i::div (a, b), b));
# We will never call floor with an inf or nan
#
fun floor x
=
if (x < 1073741824.0 and x >= -1073741824.0)
runtime::asm::floor x;
elif (is_nan x) raise exception exceptions_guts::DOMAIN; # exceptions_guts is from src/lib/std/src/exceptions-guts.pkg else raise exception exceptions_guts::OVERFLOW;
fi;
real = inline_t::f64::from_tagged_int;
# This is the IEEE double-precision maxint;
# won't work accurately on VAX
#
maxint = 4503599627370496.0;
# realround (x) returns x rounded to some nearby integer, almost always
# the nearest integer.
# May be applied to inf's and nan's.
#
fun realround x
=
x >= 0.0 ?? x+maxint-maxint
:: x-maxint+maxint;
pio4 = 7.853981633974483096E-1;
pio2 = 1.5707963267948966192E0;
pi3o4 = 2.3561944901923449288E0;
# PI = pi;
pi2 = 6.2831853071795864769E0;
one_over2pi = 0.1591549430918953357688837633725143620345;
stipulate
p1 = 1.3887401997267371720E-2;
p2 = 3.3044019718331897649E-5;
q1 = 1.1110813732786649355E-1;
q2 = 9.9176615021572857300E-4;
herein
fun exp__e (x: Float, c: Float)
=
{ z = x*x;
p = z*(p1+z*p2);
q = z*(q1+z*q2);
xp= x*p;
xh= x*half;
w = xh-(q-xp);
c = c+x*((xh*w-(q-(p+p+xp)))/(one-w)+c);
z*half+c;
};
end;
# For exp and ln
ln2hi = 6.9314718036912381649E-1;
ln2lo = 1.9082149292705877000E-10;
sqrt2 = 1.4142135623730951455E0;
lnhuge = 7.1602103751842355450E2;
lntiny = -7.5137154372698068983E2;
invln2 = 1.4426950408889633870E0;
fun exp (x: Float) # Propagates and generates inf's and nan's correctly
=
{ fun exp_norm x
=
{ # Argument reduction: x --> x - k*ln2
k = floor (invln2*x+copysign (half, x)); # k=NINT (x/ln2)
kkk = real k;
# express x-k*ln2 as z+c
hi = x-kkk*ln2hi;
lo = kkk*ln2lo;
z = hi - lo;
c = (hi-z)-lo;
# return 2^k*[expm1 (x) + 1]
z = z + exp__e (z, c);
scalb (z+one, k);
};
if (x <= lnhuge)
x >= lntiny ?? exp_norm x
:: zero;
elif (is_nan x) x;
else plus_infinity;
fi;
};
stipulate
c1 = 6.6666666666667340202E-1;
c2 = 3.9999999999416702146E-1;
c3 = 2.8571428742008753154E-1;
c4 = 2.2222198607186277597E-1;
c5 = 1.8183562745289935658E-1;
c6 = 1.5314087275331442206E-1;
c7 = 1.4795612545334174692E-1;
herein
fun log__l z
=
z*(c1+z*(c2+z*(c3+z*(c4+z*(c5+z*(c6+z*c7))))));
end;
fun ln (x: Float) # handles inf's and nan's correctly
=
if (x>0.0)
if (x < plus_infinity)
k = logb (x);
x = scalb (x, (i::neg) k);
my (k, x)
=
x >= sqrt2 ?? ((i::(+))(k, 1), x*half)
:: (k, x);
kkk = real k;
x = x - one;
# Compute log (1+x)
s = x/(two+x);
t = x*x*half;
z = kkk*ln2lo+s*(t+log__l (s*s));
x = x + (z - t);
kkk*ln2hi+x;
else
x;
fi;
elif ((x ==== 0.0))
minus_infinity;
elif (is_nan x)
x;
else
na_n;
fi;
one_overln10 = 1.0 / ln 10.0;
fun log10 x
=
ln x * one_overln10;
fun is_int y
=
realround (y)-y ==== 0.0;
fun is_odd_int y
=
is_int((y - 1.0)*0.5);
fun intpow (x, 0)
=>
1.0;
intpow (x, y)
=>
{ h = i::rshift (y, 1);
z = intpow (x, h);
zz = z*z;
if (y==(i::(+))(h, h)) zz;
else x*zz;
fi;
};
end;
# We do not properly handle negative zeros. XXX BUGGO FIXME
# Also, the copysign function works incorrectly on negative zero. XXX BUGGO FIXME
# The code for "pow" below should work correctly when these other
# bugs are fixed. Andrew. Appel, 5/8/97 */
#
fun pow (x, y)
=
if (y > 0.0)
if (y < plus_infinity)
if (x > minus_infinity)
if (x > 0.0)
exp (y*ln (x));
else
if (x ==== 0.0)
is_odd_int y ?? x
:: 0.0;
else
is_int y ?? intpow (x, floor (y+0.5))
:: na_n;
fi;
fi;
else
if (is_nan x) x;
elif (is_odd_int y) x;
else plus_infinity;
fi;
fi;
else
ax = abs (x);
if (ax > 1.0) plus_infinity;
elif (ax < 1.0) 0.0;
else na_n;
fi;
fi;
else
if (y < 0.0)
if (y > minus_infinity)
if (x > minus_infinity)
if (x > 0.0) exp (y*ln (x));
elif (x====0.0)
if (is_odd_int y) copysign (plus_infinity, x);
else plus_infinity;
fi;
else
if (is_int y) 1.0 / intpow (x, floor(-y+0.5));
else na_n;
fi;
fi;
else
if (is_nan x) x;
elif (is_odd_int y) -0.0;
else 0.0;
fi;
fi;
else
ax = abs (x);
if (ax > 1.0) 0.0;
elif (ax < 1.0) plus_infinity;
else na_n;
fi;
fi;
else
is_nan y ?? y
:: 1.0;
fi;
fi;
my (**) = pow;
stipulate
athfhi = 4.6364760900080611433E-1;
athflo = 1.0147340032515978826E-18;
at1hi = 0.78539816339744830676;
at1lo = 1.11258708870781088040E-18;
a1 = 3.3333333333333942106E-1;
a2 = -1.9999999999979536924E-1;
a3 = 1.4285714278004377209E-1;
a4 = -1.1111110579344973814E-1;
a5 = 9.0908906105474668324E-2;
a6 = -7.6919217767468239799E-2;
a7 = 6.6614695906082474486E-2;
a8 = -5.8358371008508623523E-2;
a9 = 4.9850617156082015213E-2;
a10 = -3.6700606902093604877E-2;
a11 = 1.6438029044759730479E-2;
fun atn (t, hi, lo) # for -0.4375 <= t <= 0.4375
=
{ z = t*t;
hi+(t+(lo-t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+
z*(a8+z*(a9+z*(a10+z*a11)))))))))))));
};
fun atan (t) # 0 <= t <= 1
=
if (t <= 0.4375 ) atn (t, zero, zero);
elif (t <= 0.6875 ) atn((t-half)/(one+half*t), athfhi, athflo);
else atn((t-one)/(one+t), at1hi, at1lo);
fi;
fun atanpy y # y>=0
=
if (y>one) pio2 - atan (one/y);
else atan y;
fi;
fun atan2pypx (x, y)
=
if (y > x) pio2 - atan (x/y);
else atan (y/x);
fi;
fun atan2py (x, y)
=
if (x > 0.0 ) atan2pypx (x, y);
elif (x ==== 0.0 and y ==== 0.0 ) 0.0;
else pi - atan2pypx(-x, y);
fi;
herein
fun atan y # miraculously handles inf's and nan's correctly
=
if (y <= 0.0) -(atanpy(-y));
else atanpy y;
fi;
fun atan2 (y, x) # miraculously handles inf's and nan's correctly
=
if (y >= 0.0) atan2py (x, y);
else -(atan2py (x,-y));
fi;
end;
sqrt = math_inline_t::sqrt;
sin = math_inline_t::sine;
cos = math_inline_t::cosine;
tan = math_inline_t::tangent;
fun asin x = atan2 (x, sqrt (1.0-x*x));
fun acos x = 2.0 * atan2 (sqrt((1.0-x)/(1.0+x)), 1.0);
fun cosh u
=
{ a = exp u;
if (a====0.0)
plus_infinity;
else 0.5 * (a + 1.0 / a);
fi;
};
fun sinh u
=
{ a = exp u;
if (a ==== 0.0) copysign (plus_infinity, u);
else 0.5 * (a - 1.0 / a);
fi;
};
fun tanh u
=
{ a = exp u;
b = 1.0 / a;
if (a====0.0) copysign (1.0, u);
else (a-b) / (a+b);
fi;
};
};