summaryrefslogtreecommitdiff
path: root/lib/math/exp-impl.myr
diff options
context:
space:
mode:
Diffstat (limited to 'lib/math/exp-impl.myr')
-rw-r--r--lib/math/exp-impl.myr73
1 files changed, 36 insertions, 37 deletions
diff --git a/lib/math/exp-impl.myr b/lib/math/exp-impl.myr
index d7f07d9..65cc868 100644
--- a/lib/math/exp-impl.myr
+++ b/lib/math/exp-impl.myr
@@ -10,15 +10,15 @@ use "util"
enough to be in the same function.
*/
pkg math =
- pkglocal const exp32 : (f : flt32 -> flt32)
- pkglocal const exp64 : (f : flt64 -> flt64)
+ pkglocal const exp32 : (x : flt32 -> flt32)
+ pkglocal const exp64 : (x : flt64 -> flt64)
- pkglocal const expm132 : (f : flt32 -> flt32)
- pkglocal const expm164 : (f : flt64 -> flt64)
+ pkglocal const expm132 : (x : flt32 -> flt32)
+ pkglocal const expm164 : (x : flt64 -> flt64)
;;
-extern const horner_polyu32 : (f : flt32, a : uint32[:] -> flt32)
-extern const horner_polyu64 : (f : flt64, a : uint64[:] -> flt64)
+extern const horner_polyu32 : (x : flt32, a : uint32[:] -> flt32)
+extern const horner_polyu64 : (x : flt64, a : uint64[:] -> flt64)
type fltdesc(@f, @u, @i) = struct
explode : (f : @f -> (bool, @i, @u))
@@ -50,7 +50,6 @@ type fltdesc(@f, @u, @i) = struct
L1 : @u
L2 : @u
S : (@u, @u)[32]
-
;;
const desc32 : fltdesc(flt32, uint32, int32) = [
@@ -194,18 +193,18 @@ const desc64 : fltdesc(flt64, uint64, int64) = [
.precision = 53,
]
-const exp32 = {f : flt32
- -> expgen(f, desc32)
+const exp32 = {x : flt32
+ -> expgen(x, desc32)
}
-const exp64 = {f : flt64
- -> expgen(f, desc64)
+const exp64 = {x : flt64
+ -> expgen(x, desc64)
}
-generic expgen = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i
- var b = d.tobits(f)
+generic expgen = {x : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i
+ var b = d.tobits(x)
var n, e, s
- (n, e, s) = d.explode(f)
+ (n, e, s) = d.explode(x)
/*
Detect if exp(f) would round to outside representability.
@@ -226,7 +225,7 @@ generic expgen = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equata
/* Argument reduction to [ -ln(2)/64, ln(2)/64 ] */
var inv_L = d.frombits(d.inv_L)
- var N = rn(f * inv_L)
+ var N = rn(x * inv_L)
var N2 = N % (32 : @i)
if N2 < 0
N2 += (32 : @i)
@@ -242,9 +241,9 @@ generic expgen = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equata
(very well) f reduced into [ -ln(2)/64, ln(2)/64 ]
*/
if std.abs(N) >= (1 << d.nabs)
- R1 = (f - (N1 : @f) * d.frombits(d.L1)) - ((N2 : @f) * d.frombits(d.L1))
+ R1 = (x - (N1 : @f) * d.frombits(d.L1)) - ((N2 : @f) * d.frombits(d.L1))
else
- R1 = f - (N : @f) * d.frombits(d.L1)
+ R1 = x - (N : @f) * d.frombits(d.L1)
;;
R2 = -1.0 * (N : @f) * d.frombits(d.L2)
@@ -275,33 +274,33 @@ generic expgen = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equata
-> exp
}
-const expm132 = {f : flt32
- -> expm1gen(f, desc32)
+const expm132 = {x : flt32
+ -> expm1gen(x, desc32)
}
-const expm164 = {f : flt64
- -> expm1gen(f, desc64)
+const expm164 = {x : flt64
+ -> expm1gen(x, desc64)
}
-generic expm1gen = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i
- var b = d.tobits(f)
+generic expm1gen = {x : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i
+ var b = d.tobits(x)
var n, e, s
- (n, e, s) = d.explode(f)
+ (n, e, s) = d.explode(x)
/* Special cases: +/- 0, inf, NaN, tiny, and huge */
if (b & ~d.sgnmask == 0)
- -> f
+ -> x
elif n && (b & ~d.sgnmask == d.inf)
-> (-1.0 : @f)
elif (b & ~d.sgnmask == d.inf)
- -> f
- elif std.isnan(f)
+ -> x
+ elif std.isnan(x)
-> d.frombits(d.nan)
elif (b & ~d.sgnmask) <= d.thresh_tiny
var two_to_large = d.assem(false, 100, 0)
var two_to_small = d.assem(false, -100, 0)
- var abs_f = d.assem(false, e, s)
- -> (two_to_large * f + abs_f) * two_to_small
+ var abs_x = d.assem(false, e, s)
+ -> (two_to_large * x + abs_x) * two_to_small
elif !n && b >= d.thresh_1_max /* exp(x) = oo <=> expm1(x) = oo, as it turns out */
-> d.frombits(d.inf)
elif n && b >= d.thresh_huge_neg
@@ -312,25 +311,25 @@ generic expm1gen = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equa
/* Procedure 2 */
/* compute x^2 / 2 with extra precision */
- var u = round(f, d)
- var v = f - u
+ var u = round(x, d)
+ var v = x - u
var y = u * u * (0.5 : @f)
- var z = v * (f + u) * (0.5 : @f)
- var q = f * f * f * d.horner(f, d.Bi)
+ var z = v * (x + u) * (0.5 : @f)
+ var q = x * x * x * d.horner(x, d.Bi)
var yn, ye, ys
(yn, ye, ys) = d.explode(y)
if (ye >= -7)
-> (u + y) + (q + (v + z))
else
- -> f + (y + (q + z))
+ -> x + (y + (q + z))
;;
;;
/* Procedure 1 */
var inv_L = d.frombits(d.inv_L)
- var N = rn(f * inv_L)
+ var N = rn(x * inv_L)
var N2 = N % (32 : @i)
if N2 < 0
N2 += (32 : @i)
@@ -344,9 +343,9 @@ generic expm1gen = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equa
reduced into [ -ln(2)/64, ln(2)/64 ]
*/
if std.abs(N) >= (1 << d.nabs)
- R1 = (f - (N1 : @f) * d.frombits(d.L1)) - ((N2 : @f) * d.frombits(d.L1))
+ R1 = (x - (N1 : @f) * d.frombits(d.L1)) - ((N2 : @f) * d.frombits(d.L1))
else
- R1 = f - (N : @f) * d.frombits(d.L1)
+ R1 = x - (N : @f) * d.frombits(d.L1)
;;
R2 = -1.0 * (N : @f) * d.frombits(d.L2)