ref: f762f96f07cf2dfc7b3bb745d6a767dccda0a01a
parent: 56a659f634f2a591cf641eb875a4e7f457c40321
author: S. Gilles <sgilles@math.umd.edu>
date: Tue Jul 31 09:20:50 EDT 2018
Negate atan properly when atan(1/x) is computed.
--- a/lib/math/atan-impl.myr
+++ b/lib/math/atan-impl.myr
@@ -338,9 +338,16 @@
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 @@
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 @@
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 @@
(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[:])
}
--- /dev/null
+++ b/lib/math/test/atan-impl.myr
@@ -1,0 +1,88 @@
+use std
+use math
+use testr
+
+const main = {+ math.fptrap(false)
+ testr.run([
+ [.name="atan-01", .fn = atan01], /* atan, flt32 */
+ [.name="atan-02", .fn = atan02], /* atan, flt64 */
+ [.name="atan-03", .fn = atan03], /* atan2, flt32 */
+ [.name="atan-04", .fn = atan04], /* atan2, flt64 */
+ [.name="atan-05", .fn = atan05], /* off-by-1-ulp quarantine */
+ [.name="atan-06", .fn = atan06], /* exhaustively test C */
+ [.name="atan-07", .fn = atan07], /* NaN handling */
+ ][:])
+}
+
+const same32 = {a, b+ if a == b
+ -> true
+ ;;
+
+ if std.isnan(std.flt32frombits(a)) && std.isnan(std.flt32frombits(b))
+ -> true
+ ;;
+
+ -> false
+}
+
+const same64 = {a, b+ if a == b
+ -> true
+ ;;
+
+ if std.isnan(std.flt64frombits(a)) && std.isnan(std.flt64frombits(b))
+ -> true
+ ;;
+
+ -> false
+}
+
+const atan01 = {c+ var inputs : (uint32, uint32)[:] = [
+ (0x00000000, 0x00000000),
+ (0xec0c0000, 0xbfc90fdb),
+ ][:]
+
+ for (x, y) : inputs
+ var xf : flt32 = std.flt32frombits(x)
+ var yf : flt32 = std.flt32frombits(y)
+ var a1f : flt32 = math.atan(xf)
+ var a2f : flt32 = math.atan2(xf, 1.0)
+ var a1u : uint32 = std.flt32bits(a1f)
+ var a2u : uint32 = std.flt32bits(a2f)
+ testr.check(c, same32(a1u, a2u),
+ "atan(0x{b=16,w=8,p=0}) = 0x{b=16,w=8,p=0}, but atan2(0x{b=16,w=8,p=0}, 1.0) = 0x{b=16,w=8,p=0}",+ x, a1u, x, a2u)
+
+ testr.check(c, same32(a1u, y),
+ "atan(0x{b=16,w=8,p=0}) = 0x{b=16,w=8,p=0}, should be 0x{b=16,w=8,p=0}",+ x, a1u, y)
+ ;;
+}
+
+const atan02 = {c+}
+
+const atan03 = {c+}
+
+const atan04 = {c+}
+
+const atan05 = {c+}
+
+const atan06 = {c+}
+
+const atan07 = {c+ testr.check(c, std.isnan(math.atan(std.flt64nan())), "atan(NaN64) should be NaN")
+ testr.check(c, std.isnan(math.atan(std.flt32nan())), "atan(NaN32) should be NaN")
+
+ testr.check(c, std.isnan(math.atan2(std.flt64nan(), 3.0)), "atan2(NaN64, 3.0) should be NaN")
+ testr.check(c, std.isnan(math.atan2(std.flt32nan(), -2.0)), "atan2(NaN32, -2.0) should be NaN")
+ testr.check(c, std.isnan(math.atan2(6.0, std.flt64nan())), "atan2(6.0, NaN64) should be NaN")
+ testr.check(c, std.isnan(math.atan2(4.0, std.flt32nan())), "atan2(4.0, NaN32) should be NaN")
+}