summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
authorS. Gilles <sgilles@math.umd.edu>2019-05-21 03:51:51 -0400
committerS. Gilles <sgilles@math.umd.edu>2019-05-21 03:51:51 -0400
commit895c8ae58de3ece766a1f9d70094f61adb864f71 (patch)
treeb4e10ee21ae2e6b83ca0e2b930dd5ea12a7aa5a2 /lib
parenta00e33e21435b8a4961d6e0f910e09f9c17f65ba (diff)
downloadmc-895c8ae58de3ece766a1f9d70094f61adb864f71.tar.gz
Correct sign-handling for pown special case
Diffstat (limited to 'lib')
-rw-r--r--lib/math/ancillary/generate-minimax-by-Remez.gp10
-rw-r--r--lib/math/bld.sub2
-rw-r--r--lib/math/pown-impl.myr2
-rw-r--r--lib/math/test/pown-impl.myr50
4 files changed, 45 insertions, 19 deletions
diff --git a/lib/math/ancillary/generate-minimax-by-Remez.gp b/lib/math/ancillary/generate-minimax-by-Remez.gp
index 2863120..14041ff 100644
--- a/lib/math/ancillary/generate-minimax-by-Remez.gp
+++ b/lib/math/ancillary/generate-minimax-by-Remez.gp
@@ -136,6 +136,10 @@
return(1/tanoverx(x));
}
+{ log2xoverx(x) =
+ return(if(x == 1,1,log(x)/(x-1))/log(2));
+}
+
print("\n");
print("Minimaxing sin(x) / x, degree 6, on [-Pi/(4 * 256), Pi/(4 * 256)]:");
find_minimax(sinoverx, 6, -Pi/1024, Pi/1024)
@@ -176,4 +180,10 @@ print("\n");
print("(You'll need to add a 0x0 at the beginning to make a degree 13...\n");
print("\n");
print("---\n");
+print("Minmimaxing log_2(x) / (x - 1), degree 7, on [1, 2^(1/8)]:");
+find_minimax(log2xoverx, 7, 1, 2^(1/8))
+print("\n");
+/* print("(You'll need to add a 0x0 at the beginning to make a degree 13...\n"); */
+/* print("\n"); */
+print("---\n");
print("Remember that there's that extra, ugly E term at the end of the vector that you want to lop off.\n");
diff --git a/lib/math/bld.sub b/lib/math/bld.sub
index 11bceeb..3463e3f 100644
--- a/lib/math/bld.sub
+++ b/lib/math/bld.sub
@@ -23,7 +23,9 @@ lib math =
poly-impl.myr
# x^n and x^1/q
+ log-overkill+posixy-x64-sse2.s
log-overkill.myr
+ pown-impl+posixy-x64-sse2.s
pown-impl.myr
# x^y
diff --git a/lib/math/pown-impl.myr b/lib/math/pown-impl.myr
index 2feecb5..a932541 100644
--- a/lib/math/pown-impl.myr
+++ b/lib/math/pown-impl.myr
@@ -109,7 +109,7 @@ generic powngen = {x : @f, n : @i, d : fltdesc(@f, @u, @i) :: numeric,floating,s
/* Propagate NaN (why doesn't this come first? Ask IEEE.) */
-> d.frombits(d.nan)
elif (x == 0.0 || x == -0.0)
- if n < 0 && (n % 2 == 1) && xn
+ if n < 0 && (n % 2 == -1) && xn
/* (+/- 0)^n = +/- oo */
-> d.frombits(d.neginf)
elif n < 0
diff --git a/lib/math/test/pown-impl.myr b/lib/math/test/pown-impl.myr
index a7dbeb4..c50cb35 100644
--- a/lib/math/test/pown-impl.myr
+++ b/lib/math/test/pown-impl.myr
@@ -21,7 +21,20 @@ const int64fromuint64 = {b; -> (&b : int64#)#}
const pown01 = {c
var inputs : (uint32, uint32, uint32)[:] = [
(0x000000f6, 0x00000000, 0x3f800000),
+ (0x7fc00000, 0x00000001, 0x7fc00000),
+ (0x7fc00000, 0x00000021, 0x7fc00000),
+ (0x7fc00000, 0x00030021, 0x7fc00000),
+ (0x7fc00000, 0xfacecafe, 0x7fc00000),
(0x00000000, 0x3d800000, 0x00000000),
+ (0x80000000, 0x00000124, 0x00000000),
+ (0x80000000, 0x00000123, 0x80000000),
+ (0x00000000, 0x00000124, 0x00000000),
+ (0x00000000, 0x00000123, 0x00000000),
+ (0x00000000, 0xad800001, 0x7f800000),
+ (0x80000000, 0x80000123, 0xff800000),
+ (0x80000000, 0x80000122, 0x7f800000),
+ (0x00000000, 0x80000127, 0x7f800000),
+ (0x00000000, 0x80000128, 0x7f800000),
(0x946fc13b, 0x3b21efc7, 0x80000000),
(0xb76e98b6, 0xdbeb6637, 0xff800000),
(0xc04825b7, 0x53cdd772, 0x7f800000),
@@ -52,11 +65,11 @@ const pown01 = {c
for (x, y, z) : inputs
var xf : flt32 = std.flt32frombits(x)
var yi : int32 = int32fromuint32(y)
- var zf : flt32 = std.flt32frombits(z)
var rf = math.pown(xf, yi)
- testr.check(c, rf == zf,
+ var ru : uint32 = std.flt32bits(rf)
+ testr.check(c, ru == z,
"pown(0x{b=16,w=8,p=0}, {}) should be 0x{b=16,w=8,p=0}, was 0x{b=16,w=8,p=0}",
- x, yi, z, std.flt32bits(rf))
+ x, yi, z, ru)
;;
}
@@ -92,11 +105,11 @@ const pown02 = {c
for (x, y, z) : inputs
var xf : flt64 = std.flt64frombits(x)
var yi : int64 = int64fromuint64(y)
- var zf : flt64 = std.flt64frombits(z)
- var rf = math.pown(xf, yi)
- testr.check(c, rf == zf,
+ var rf : flt64 = math.pown(xf, yi)
+ var ru : uint64 = std.flt64bits(rf)
+ testr.check(c, ru == z,
"pown(0x{b=16,w=16,p=0}, {}) should be 0x{b=16,w=16,p=0}, was 0x{b=16,w=16,p=0}",
- x, yi, z, std.flt64bits(rf))
+ x, yi, z, ru)
;;
}
@@ -112,12 +125,11 @@ const pown03 = {c
for (x, y, z_perfect, z_accepted) : inputs
var xf : flt32 = std.flt32frombits(x)
var yi : int32 = int32fromuint32(y)
- var zf_perfect : flt32 = std.flt32frombits(z_perfect)
- var zf_accepted : flt32 = std.flt32frombits(z_accepted)
- var rf = math.pown(xf, yi)
- testr.check(c, rf == zf_perfect || rf == zf_accepted,
+ var rf : flt32 = math.pown(xf, yi)
+ var ru : uint32 = std.flt32bits(rf)
+ testr.check(c, ru == z_perfect || ru == z_accepted,
"pown(0x{b=16,w=8,p=0}, {}) should be 0x{b=16,w=8,p=0}, will also accept 0x{b=16,w=8,p=0}, was 0x{b=16,w=8,p=0}",
- x, yi, z_perfect, z_accepted, std.flt32bits(rf))
+ x, yi, z_perfect, z_accepted, ru)
;;
}
@@ -279,10 +291,11 @@ const rootn01 = {c
for (x, y, z) : inputs
var xf : flt32 = std.flt32frombits(x)
var zf : flt32 = std.flt32frombits(z)
- var rf = math.rootn(xf, y)
- testr.check(c, rf == zf,
+ var rf : flt32 = math.rootn(xf, y)
+ var ru : uint32 = std.flt32bits(rf)
+ testr.check(c, ru == z,
"rootn(0x{b=16,w=8,p=0}, {}) should be 0x{b=16,w=8,p=0}, was 0x{b=16,w=8,p=0}",
- x, y, z, std.flt32bits(rf))
+ x, y, z, ru)
;;
}
@@ -441,10 +454,11 @@ const rootn02 = {c
for (x, y, z) : inputs
var xf : flt64 = std.flt64frombits(x)
var zf : flt64 = std.flt64frombits(z)
- var rf = math.rootn(xf, y)
- testr.check(c, rf == zf,
+ var rf : flt64 = math.rootn(xf, y)
+ var ru : uint64 = std.flt64bits(rf)
+ testr.check(c, ru == z,
"rootn(0x{b=16,w=16,p=0}, {}) should be 0x{b=16,w=16,p=0}, was 0x{b=16,w=16,p=0}",
- x, y, z, std.flt64bits(rf))
+ x, y, z, ru)
;;
}