authorgravatar for marc@tiehu.isMarc Tiehuis <marc@tiehu.is> 2018-09-13 20:33:05+12:00
committergravatar for marc@tiehu.isMarc Tiehuis <marc@tiehu.is> 2018-09-13 20:33:05+12:00
loge70c543bc4c097d7c3287feb0f233a13d4dc3829
treeeed3f73ee9e52cbeb801fa85b407b5342d0e26f6
parentafe6316d32d0e61687c58d152207252bdc33dad7

math/complex: cexp test correction and ldexp usage fix

3 files changed, 8 insertions(+), 9 deletions(-)

std/math/complex/cosh.zig+2-2
...@@ -44,7 +44,7 @@ fn cosh32(z: *const Complex(f32)) Complex(f32) {...@@ -44,7 +44,7 @@ fn cosh32(z: *const Complex(f32)) Complex(f32) {
44 else if (ix < 0x4340b1e7) {44 else if (ix < 0x4340b1e7) {
45 const v = Complex(f32).new(math.fabs(x), y);45 const v = Complex(f32).new(math.fabs(x), y);
46 const r = ldexp_cexp(v, -1);46 const r = ldexp_cexp(v, -1);
47 return Complex(f32).new(x, y * math.copysign(f32, 1, x));47 return Complex(f32).new(r.re, r.im * math.copysign(f32, 1, x));
48 }48 }
49 // x >= 192.7: result always overflows49 // x >= 192.7: result always overflows
50 else {50 else {
...@@ -112,7 +112,7 @@ fn cosh64(z: *const Complex(f64)) Complex(f64) {...@@ -112,7 +112,7 @@ fn cosh64(z: *const Complex(f64)) Complex(f64) {
112 else if (ix < 0x4096bbaa) {112 else if (ix < 0x4096bbaa) {
113 const v = Complex(f64).new(math.fabs(x), y);113 const v = Complex(f64).new(math.fabs(x), y);
114 const r = ldexp_cexp(v, -1);114 const r = ldexp_cexp(v, -1);
115 return Complex(f64).new(x, y * math.copysign(f64, 1, x));115 return Complex(f64).new(r.re, r.im * math.copysign(f64, 1, x));
116 }116 }
117 // x >= 1455: result always overflows117 // x >= 1455: result always overflows
118 else {118 else {
std/math/complex/exp.zig+4-5
...@@ -69,7 +69,7 @@ fn exp64(z: Complex(f64)) Complex(f64) {...@@ -69,7 +69,7 @@ fn exp64(z: Complex(f64)) Complex(f64) {
69 const y = z.im;69 const y = z.im;
7070
71 const fy = @bitCast(u64, y);71 const fy = @bitCast(u64, y);
72 const hy = u32(fy >> 32) & 0x7fffffff;72 const hy = @intCast(u32, (fy >> 32) & 0x7fffffff);
73 const ly = @truncate(u32, fy);73 const ly = @truncate(u32, fy);
7474
75 // cexp(x + i0) = exp(x) + i075 // cexp(x + i0) = exp(x) + i0
...@@ -78,7 +78,7 @@ fn exp64(z: Complex(f64)) Complex(f64) {...@@ -78,7 +78,7 @@ fn exp64(z: Complex(f64)) Complex(f64) {
78 }78 }
7979
80 const fx = @bitCast(u64, x);80 const fx = @bitCast(u64, x);
81 const hx = u32(fx >> 32);81 const hx = @intCast(u32, fx >> 32);
82 const lx = @truncate(u32, fx);82 const lx = @truncate(u32, fx);
8383
84 // cexp(0 + iy) = cos(y) + isin(y)84 // cexp(0 + iy) = cos(y) + isin(y)
...@@ -101,8 +101,7 @@ fn exp64(z: Complex(f64)) Complex(f64) {...@@ -101,8 +101,7 @@ fn exp64(z: Complex(f64)) Complex(f64) {
101101
102 // 709.7 <= x <= 1454.3 so must scale102 // 709.7 <= x <= 1454.3 so must scale
103 if (hx >= exp_overflow and hx <= cexp_overflow) {103 if (hx >= exp_overflow and hx <= cexp_overflow) {
104 const r = ldexp_cexp(z, 0);104 return ldexp_cexp(z, 0);
105 return r.*;
106 } // - x < exp_overflow => exp(x) won't overflow (common)105 } // - x < exp_overflow => exp(x) won't overflow (common)
107 // - x > cexp_overflow, so exp(x) * s overflows for s > 0106 // - x > cexp_overflow, so exp(x) * s overflows for s > 0
108 // - x = +-inf107 // - x = +-inf
...@@ -124,7 +123,7 @@ test "complex.cexp32" {...@@ -124,7 +123,7 @@ test "complex.cexp32" {
124}123}
125124
126test "complex.cexp64" {125test "complex.cexp64" {
127 const a = Complex(f32).new(5, 3);126 const a = Complex(f64).new(5, 3);
128 const c = exp(a);127 const c = exp(a);
129128
130 debug.assert(math.approxEq(f64, c.re, -146.927917, epsilon));129 debug.assert(math.approxEq(f64, c.re, -146.927917, epsilon));
std/math/complex/sinh.zig+2-2
...@@ -44,7 +44,7 @@ fn sinh32(z: Complex(f32)) Complex(f32) {...@@ -44,7 +44,7 @@ fn sinh32(z: Complex(f32)) Complex(f32) {
44 else if (ix < 0x4340b1e7) {44 else if (ix < 0x4340b1e7) {
45 const v = Complex(f32).new(math.fabs(x), y);45 const v = Complex(f32).new(math.fabs(x), y);
46 const r = ldexp_cexp(v, -1);46 const r = ldexp_cexp(v, -1);
47 return Complex(f32).new(x * math.copysign(f32, 1, x), y);47 return Complex(f32).new(r.re * math.copysign(f32, 1, x), r.im);
48 }48 }
49 // x >= 192.7: result always overflows49 // x >= 192.7: result always overflows
50 else {50 else {
...@@ -111,7 +111,7 @@ fn sinh64(z: Complex(f64)) Complex(f64) {...@@ -111,7 +111,7 @@ fn sinh64(z: Complex(f64)) Complex(f64) {
111 else if (ix < 0x4096bbaa) {111 else if (ix < 0x4096bbaa) {
112 const v = Complex(f64).new(math.fabs(x), y);112 const v = Complex(f64).new(math.fabs(x), y);
113 const r = ldexp_cexp(v, -1);113 const r = ldexp_cexp(v, -1);
114 return Complex(f64).new(x * math.copysign(f64, 1, x), y);114 return Complex(f64).new(r.re * math.copysign(f64, 1, x), r.im);
115 }115 }
116 // x >= 1455: result always overflows116 // x >= 1455: result always overflows
117 else {117 else {