summaryrefslogtreecommitdiff
path: root/lib/math/scale2-impl.myr
diff options
context:
space:
mode:
authorS. Gilles <sgilles@math.umd.edu>2018-05-03 11:13:58 -0400
committerS. Gilles <sgilles@math.umd.edu>2018-05-03 12:56:36 -0400
commitbf035b686f58f139ee0b04746d3215e7bb8413e9 (patch)
tree085f9e79fe6df5766060708fd12baf7f290b1b01 /lib/math/scale2-impl.myr
parent0dd6c9d5124707e4c69116745469e21f7e076429 (diff)
downloadmc-bf035b686f58f139ee0b04746d3215e7bb8413e9.tar.gz
Implment scale2.
This is, perhaps, a conforming implementation of scalB, but it is highly probable that, in the future, it will be desireable for scalB to accept arbitrary, exponents, not simply integers. Furthermore, we find scalB a non-intuitive name. Therefore, calling the function ‘scale2’ serves two purposes.
Diffstat (limited to 'lib/math/scale2-impl.myr')
-rw-r--r--lib/math/scale2-impl.myr73
1 files changed, 73 insertions, 0 deletions
diff --git a/lib/math/scale2-impl.myr b/lib/math/scale2-impl.myr
new file mode 100644
index 0000000..6301e14
--- /dev/null
+++ b/lib/math/scale2-impl.myr
@@ -0,0 +1,73 @@
+use std
+
+use "fpmath"
+use "util"
+
+/*
+ The scaleB function recommended by the 2008 revision of IEEE
+ 754. Since only integer exponents are considered, the naive
+ approach works quite well.
+ */
+pkg math =
+ const scale232 : (f : flt32, m : int32 -> flt32)
+ const scale264 : (f : flt64, m : int64 -> flt64)
+;;
+
+const scale232 = {f : flt32, m : int32
+ var n, e, s
+ (n, e, s) = std.flt32explode(f)
+ (n, e, s) = scale2gen(n, e, s, -126, 127, 24, m)
+ -> std.flt32assem(n, e, s)
+}
+
+const scale264 = {f : flt64, m : int64
+ var n, e, s
+ (n, e, s) = std.flt64explode(f)
+ (n, e, s) = scale2gen(n, e, s, -1022, 1023, 53, m)
+ -> std.flt64assem(n, e, s)
+}
+
+generic scale2gen = {n : bool, e : @i, s : @u, emin : @i, emax : @i, p : @u, m : @i :: numeric, integral @i, numeric, integral @u
+ if e < emin && s == 0
+ -> (n, e, s)
+ elif m == 0
+ -> (n, e, s)
+ elif m < 0
+ var sprime = s
+ var eprime = e
+ if e < emin
+ sprime = s >> (-m)
+ if need_round_away(0, (s : uint64), (-m : int64))
+ sprime++
+ ;;
+ eprime = emin - 1
+ elif e + m < emin
+ sprime = s >> (emin - m - e)
+ if need_round_away(0, (s : uint64), ((emin - m - e) : int64))
+ sprime++
+ ;;
+ eprime = emin - 1
+ else
+ eprime = e + m
+ ;;
+ -> (n, eprime, sprime)
+ ;;
+
+ if e < emin
+ var first_1 = find_first1_64((s : uint64), (p : int64))
+ var shift = p - (first_1 : @u) - 1
+ if m >= (shift : @i)
+ s = s << shift
+ m -= (shift : @i)
+ e = emin
+ else
+ -> (n, e, s << m)
+ ;;
+ ;;
+
+ if e + m > emax && s != 0
+ -> (n, emax + 1, 1 << (p - 1))
+ ;;
+
+ -> (n, e + m, s)
+}