summaryrefslogtreecommitdiff
path: root/lib/math/atan-impl.myr
diff options
context:
space:
mode:
authorS. Gilles <sgilles@math.umd.edu>2018-07-31 13:20:50 -0400
committerS. Gilles <sgilles@math.umd.edu>2018-07-31 13:20:50 -0400
commitf762f96f07cf2dfc7b3bb745d6a767dccda0a01a (patch)
tree9a14f5fb98f84b53634dcbd186d30d690a13461f /lib/math/atan-impl.myr
parent56a659f634f2a591cf641eb875a4e7f457c40321 (diff)
downloadmc-f762f96f07cf2dfc7b3bb745d6a767dccda0a01a.tar.gz
Negate atan properly when atan(1/x) is computed.
Diffstat (limited to 'lib/math/atan-impl.myr')
-rw-r--r--lib/math/atan-impl.myr37
1 files changed, 16 insertions, 21 deletions
diff --git a/lib/math/atan-impl.myr b/lib/math/atan-impl.myr
index 911c99c..c296a75 100644
--- a/lib/math/atan-impl.myr
+++ b/lib/math/atan-impl.myr
@@ -338,9 +338,16 @@ const atan264 = {y, x
var xabs = std.flt64frombits(std.flt64bits(x) & (~(0 : uint64)) >> 1)
var yabs = std.flt64frombits(std.flt64bits(y) & (~(0 : uint64)) >> 1)
+ var then_negate = false
+ if (xabs == x) != (yabs == y)
+ then_negate = true
+ ;;
+ x = xabs
+ y = yabs
+
var xb = std.flt64bits(x)
- if xb == 0x0 || xb == 0x8000000000000000
- if (y >= 0.0) == (x >= 0.0)
+ if xb == 0x0
+ if then_negate
-> std.flt64frombits(0x3ff921fb54442d18)
else
-> std.flt64frombits(0xbff921fb54442d18)
@@ -368,20 +375,17 @@ const atan264 = {y, x
ret = (po2_hi - ret) + po2_lo
;;
+ if then_negate
+ ret = -1.0 * ret
+ ;;
+
-> ret
}
-/* Handle arctan(z) with z in (-1, 1) */
+/* Handle arctan(z) with z in [0, 1] */
const atan_reduced = {u : flt64, du : flt64
- var then_negate : bool = false
var ret : flt64
- if u < 0.0
- then_negate = true
- u = -1.0 * u
- du = -1.0 * du
- ;;
-
/*
If u is less than 1/16, [GB91] uses a single polynomial
approximation. I am not quite sure why they don't employ
@@ -392,8 +396,7 @@ const atan_reduced = {u : flt64, du : flt64
if u < 0.0625
var s = u * u
var p = horner_polyu(s, atan_coeffs[:])
- ret = u*p + du
- goto have_result
+ -> u*p + du
;;
/*
@@ -410,13 +413,5 @@ const atan_reduced = {u : flt64, du : flt64
(xi_u, c[0], c[1], c[2], c[3], c[4], c[5]) = C[j]
var xi : flt64 = std.flt64frombits(xi_u)
- ret = horner_polyu((u - xi) + du, c[:])
-
-:have_result
-
- if then_negate
- ret = -1.0 * ret
- ;;
-
- -> ret
+ -> horner_polyu((u - xi) + du, c[:])
}