diff options
author | S. Gilles <sgilles@math.umd.edu> | 2019-05-21 03:51:51 -0400 |
---|---|---|
committer | S. Gilles <sgilles@math.umd.edu> | 2019-05-21 03:51:51 -0400 |
commit | 895c8ae58de3ece766a1f9d70094f61adb864f71 (patch) | |
tree | b4e10ee21ae2e6b83ca0e2b930dd5ea12a7aa5a2 | |
parent | a00e33e21435b8a4961d6e0f910e09f9c17f65ba (diff) | |
download | mc-895c8ae58de3ece766a1f9d70094f61adb864f71.tar.gz |
Correct sign-handling for pown special case
-rw-r--r-- | lib/math/ancillary/generate-minimax-by-Remez.gp | 10 | ||||
-rw-r--r-- | lib/math/bld.sub | 2 | ||||
-rw-r--r-- | lib/math/pown-impl.myr | 2 | ||||
-rw-r--r-- | lib/math/test/pown-impl.myr | 50 |
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) ;; } |