## math64.sml
## *************************************************************************
## *
## 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. *
## *
## *************************************************************************
# 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 true spirit of delight, the exaltation,
### the sense of being more than Man, which is the
### touchstone of the highest excellence, is to be
### found in mathematics as surely as poetry."
###
### -- Bertrand Russell
package math64: Math {
type real = real
infix 50 ====
my op + = inline_t::f64::(+)
my op - = inline_t::f64::(-)
my op * = inline_t::f64::(*)
my op / = inline_t::f64::(/)
my op (-_) = inline_t::f64::(-_)
my op neg = inline_t::f64::(-_)
my op < = inline_t::f64::(<)
my op <= = inline_t::f64::(<=)
my op > = inline_t::f64::(>)
my op >= = inline_t::f64::(>=)
op == = inline_t::f64::(====)
package i = inline_t::DfltInt
my lessu: Int * Int -> Bool = i::ltu
pi = 3.14159265358979323846
e = 2.7182818284590452354
fun isNan x
=
bool::not (x====x)
plusInfinity = 1E300 * 1E300
minusInfinity = -plusInfinity
NaN = 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) =
let j = assembly::a::logb x
k' = I.+(k, j)
in if j == -1023
then scalb (x*two_to_the_54, I.-(k, 54)) # 2
else if lessu (I.+(k', 1022), 2046)
then assembly::a::scalb (x, k) # 1
else if I.<(k', 0)
then if I.<(k', I.-(-1022, 54))
then 0.0 # 3
else scalb (x, I.+(k, 54)) *
two_to_the_minus_54 # 4
else x * plusInfinity # 5
end
/* 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 assembly::a::logb x
of -1023 => # Denormalized number
I.-(assembly::a::logb (x*two_to_the_54), 54)
| i => i
) # end case
negone = -1.0
zero = 0.0
half = 0.5
one = 1.0
two = 2.0
four = 4.0
# * SHOULD BE INLINE OP * XXX BUGGO FIXME
/* 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)
of (TRUE, true) => a
| (FALSE, false) => a
| _ => -a
) # end case
# May be applied to inf's and nan's
fun abs x
=
if x < zero then -x else x
fun op 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
then assembly::a::floor x
else if isNan x then raise exception exceptions::DOMAIN
else raise exception exceptions::OVERFLOW
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 = if x>=0.0 then x+maxint-maxint else x-maxint+maxint
# sin/cos
local
S0 = -1.6666666666666463126E-1
S1 = 8.3333333332992771264E-3
S2 = -1.9841269816180999116E-4
S3 = 2.7557309793219876880E-6
S4 = -2.5050225177523807003E-8
S5 = 1.5868926979889205164E-10
in
fun sin__S z = (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*S5))))))
end
local
C0 = 4.1666666666666504759E-2
C1 = -1.3888888888865301516E-3
C2 = 2.4801587269650015769E-5
C3 = -2.7557304623183959811E-7
C4 = 2.0873958177697780076E-9
C5 = -1.1250289076471311557E-11
in
fun cos__C z = (z*z*(C0+z*(C1+z*(C2+z*(C3+z*(C4+z*C5))))))
end
PIo4 = 7.853981633974483096E-1
PIo2 = 1.5707963267948966192E0
PI3o4 = 2.3561944901923449288E0
PI = pi
PI2 = 6.2831853071795864769E0
oneOver2Pi = 0.1591549430918953357688837633725143620345
local
thresh = 2.6117239648121182150E-1
in fun sss y = y + y * sin__S (y*y)
fun ccc y =
let yy = y*y
c = cos__C yy
Y = yy/two
in if Y < thresh then one - (Y - c)
else half - (Y - half - c)
end
end
fun sin x = # This function propagages Inf's and Nan's correctly.
let # x may be finite, inf, or nan at this point.
xover2pi = x * oneOver2Pi
x = PI2*(xover2pi - realround (xover2pi))
# now, probably, -pi <= x <= pi, except on vaxes.
# x may be a nan, but cannot now be an inf.
fun lessPIo2 x = if x>PIo4 then ccc (PIo2-x) else sss x
fun lessPI x = if x>PIo2 then lessPIo2 (PI-x) else lessPIo2 x
fun positive x = if x>PI then sin (x-PI2) # exceedingly rare
else lessPI x
in if x>=0.0
then positive x
else -(positive(-x))
end
fun cos x = sin (PIo2-x)
fun tan x = sin x / cos x
local
p1 = 1.3887401997267371720E-2
p2 = 3.3044019718331897649E-5
q1 = 1.1110813732786649355E-1
q2 = 9.9176615021572857300E-4
in fun exp__E (x: real, c: real) =
let 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)
in z*half+c
end
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: real) = # propagates and generates inf's and nan's correctly
let fun exp_norm x =
let # Argument reduction: x --> x - k*ln2
k = floor (invln2*x+copysign (half, x)) # k=NINT (x/ln2)
K = real k
# express x-k*ln2 as z+c
hi = x-K*ln2hi
lo = K*ln2lo
z = hi - lo
c = (hi-z)-lo
# return 2^k*[expm1 (x) + 1]
z = z + exp__E (z, c)
in scalb (z+one, k)
end
in if x <= lnhuge
then if x >= lntiny
then exp_norm x
else zero
else if isNan x then x else plusInfinity
end
local
L1 = 6.6666666666667340202E-1
L2 = 3.9999999999416702146E-1
L3 = 2.8571428742008753154E-1
L4 = 2.2222198607186277597E-1
L5 = 1.8183562745289935658E-1
L6 = 1.5314087275331442206E-1
L7 = 1.4795612545334174692E-1
in fun log__L (z) = z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7))))))
end
fun ln (x: real) = # handles inf's and nan's correctly
if x>0.0
then if x < plusInfinity
then let
k = logb (x)
x = scalb (x, i::(-_) k)
my (k, x) = if x >= sqrt2 then (I.+(k, 1), x*half) else (k, x)
K = real k
x = x - one
# Compute log (1+x)
s = x/(two+x)
t = x*x*half
z = K*ln2lo+s*(t+log__L (s*s))
x = x + (z - t)
in
K*ln2hi+x
end
else x
else if (x ==== 0.0)
then minusInfinity
else if isNan x then x else NaN
oneOverln10 = 1.0 / ln 10.0
fun log10 x = ln x * oneOverln10
fun isInt y = realround (y)-y ==== 0.0
fun isOddInt (y) = isInt((y - 1.0)*0.5)
fun intpow (x, 0) = 1.0
| intpow (x, y) = let h = i::rshift (y, 1)
z = intpow (x, h)
zz = z*z
in if y==I::(+) (h, h) then zz else x*zz
end
# Mythryl doesn't properly handle negative zeros. XXX BUGGO FIXME
# Also, the copysign function works incorrectly on negative zero.
# The code for "pow" below should work correctly when these other
# bugs are fixed. A. Appel, 5/8/97 */
#
fun pow (x, y) = if y>0.0
then if y<plusInfinity
then if x > minusInfinity
then if x > 0.0
then exp (y*ln (x))
else if x ==== 0.0
then if isOddInt (y)
then x
else 0.0
else if isInt (y)
then intpow (x, floor (y+0.5))
else NaN
else if isNan x
then x
else if isOddInt (y)
then x
else plusInfinity
else let ax = abs (x)
in if ax>1.0 then plusInfinity
else if ax<1.0 then 0.0
else NaN
end
else if y < 0.0
then if y>minusInfinity
then if x > minusInfinity
then if x > 0.0
then exp (y*ln (x))
else if x====0.0
then if isOddInt (y)
then copysign (plusInfinity, x)
else plusInfinity
else if isInt (y)
then 1.0 / intpow (x, floor(-y+0.5))
else NaN
else if isNan x
then x
else if isOddInt (y)
then -0.0
else 0.0
else let ax = abs (x)
in if ax>1.0 then 0.0
else if ax<1.0 then plusInfinity
else NaN
end
else if isNan y
then y
else 1.0
local
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
=
let z = t*t
in 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)))))))))))))
end
fun atan (t) = # 0 <= t <= 1
if t <= 0.4375 then atn (t, zero, zero)
else if t <= 0.6875 then atn((t-half)/(one+half*t), athfhi, athflo)
else atn((t-one)/(one+t), at1hi, at1lo)
fun atanpy y = # y>=0
if y>one then PIo2 - atan (one/y) else atan (y)
fun atan2pypx (x, y) =
if y>x then PIo2 - atan (x/y) else atan (y/x)
fun atan2py (x, y) =
if x > 0.0 then atan2pypx (x, y)
else if x ==== 0.0 and y ==== 0.0 then 0.0
else PI - atan2pypx(-x, y)
in fun atan y = # miraculously handles inf's and nan's correctly
if y<=0.0 then -(atanpy(-y)) else atanpy y
fun atan2 (y, x) = # miraculously handles inf's and nan's correctly
if y>=0.0 then atan2py (x, y) else -(atan2py (x,-y))
end
sqrt = math_inline_t::sqrt
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 = let a = exp u in if a====0.0
then plusInfinity
else 0.5 * (a + 1.0 / a)
end
fun sinh u = let a = exp u
in if a====0.0
then copysign (plusInfinity, u)
else 0.5 * (a - 1.0 / a)
end
fun tanh u = let a = exp u
b = 1.0 / a
in if a====0.0 then copysign (1.0, u)
else (a-b) / (a+b)
end
}
package math64 = Math64