| author |  LemonBoy <thatlemon@gmail.com> 2021-06-12 15:35:52+02:00 |
| committer |  Veikka Tuominen <git@vexu.eu> 2021-06-12 19:56:08+03:00 |
| log | 44cdafd9d46bed36f08effb99aa0bb7fc2145efa |
| tree | a4c988e36af2aac626644b1c080d652e943bbdc1 |
| parent | 0d022eb1d6333e6e027d35e9588b5a10035664b7 |
std: Fix complex ldexp implementation
Two bugs in the implementation ported from musl made all the complex
functions relying on ldexp return incorrect results in some cases.
Spotted in #90472 files changed, 40 insertions(+), 14 deletions(-)
lib/std/math/complex/exp.zig+32-10
| ... | @@ -124,20 +124,42 @@ fn exp64(z: Complex(f64)) Complex(f64) { | ... | @@ -124,20 +124,42 @@ fn exp64(z: Complex(f64)) Complex(f64) { |
| 124 | } | 124 | } |
| 125 | } | 125 | } |
| 126 | 126 | ||
| 127 | const epsilon = 0.0001; | ||
| 128 | |||
| 129 | test "complex.cexp32" { | 127 | test "complex.cexp32" { |
| 130 | const a = Complex(f32).init(5, 3); | 128 | const tolerance_f32 = math.sqrt(math.epsilon(f32)); |
| 131 | const c = exp(a); | 129 | |
| 130 | { | ||
| 131 | const a = Complex(f32).init(5, 3); | ||
| 132 | const c = exp(a); | ||
| 133 | |||
| 134 | try testing.expectApproxEqRel(@as(f32, -1.46927917e+02), c.re, tolerance_f32); | ||
| 135 | try testing.expectApproxEqRel(@as(f32, 2.0944065e+01), c.im, tolerance_f32); | ||
| 136 | } | ||
| 137 | |||
| 138 | { | ||
| 139 | const a = Complex(f32).init(88.8, 0x1p-149); | ||
| 140 | const c = exp(a); | ||
| 132 | 141 | ||
| 133 | try testing.expect(math.approxEqAbs(f32, c.re, -146.927917, epsilon)); | 142 | try testing.expectApproxEqAbs(math.inf(f32), c.re, tolerance_f32); |
| 134 | try testing.expect(math.approxEqAbs(f32, c.im, 20.944065, epsilon)); | 143 | try testing.expectApproxEqAbs(@as(f32, 5.15088629e-07), c.im, tolerance_f32); |
| 144 | } | ||
| 135 | } | 145 | } |
| 136 | 146 | ||
| 137 | test "complex.cexp64" { | 147 | test "complex.cexp64" { |
| 138 | const a = Complex(f64).init(5, 3); | 148 | const tolerance_f64 = math.sqrt(math.epsilon(f64)); |
| 139 | const c = exp(a); | ||
| 140 | 149 | ||
| 141 | try testing.expect(math.approxEqAbs(f64, c.re, -146.927917, epsilon)); | 150 | { |
| 142 | try testing.expect(math.approxEqAbs(f64, c.im, 20.944065, epsilon)); | 151 | const a = Complex(f64).init(5, 3); |
| 152 | const c = exp(a); | ||
| 153 | |||
| 154 | try testing.expectApproxEqRel(@as(f64, -1.469279139083189e+02), c.re, tolerance_f64); | ||
| 155 | try testing.expectApproxEqRel(@as(f64, 2.094406620874596e+01), c.im, tolerance_f64); | ||
| 156 | } | ||
| 157 | |||
| 158 | { | ||
| 159 | const a = Complex(f64).init(709.8, 0x1p-1074); | ||
| 160 | const c = exp(a); | ||
| 161 | |||
| 162 | try testing.expectApproxEqAbs(math.inf(f64), c.re, tolerance_f64); | ||
| 163 | try testing.expectApproxEqAbs(@as(f64, 9.036659362159884e-16), c.im, tolerance_f64); | ||
| 164 | } | ||
| 143 | } | 165 | } |
lib/std/math/complex/ldexp.zig+8-4
| ... | @@ -13,6 +13,7 @@ const std = @import("../../std.zig"); | ... | @@ -13,6 +13,7 @@ const std = @import("../../std.zig"); |
| 13 | const debug = std.debug; | 13 | const debug = std.debug; |
| 14 | const math = std.math; | 14 | const math = std.math; |
| 15 | const cmath = math.complex; | 15 | const cmath = math.complex; |
| 16 | const testing = std.testing; | ||
| 16 | const Complex = cmath.Complex; | 17 | const Complex = cmath.Complex; |
| 17 | 18 | ||
| 18 | /// Returns exp(z) scaled to avoid overflow. | 19 | /// Returns exp(z) scaled to avoid overflow. |
| ... | @@ -48,7 +49,10 @@ fn ldexp_cexp32(z: Complex(f32), expt: i32) Complex(f32) { | ... | @@ -48,7 +49,10 @@ fn ldexp_cexp32(z: Complex(f32), expt: i32) Complex(f32) { |
| 48 | const half_expt2 = exptf - half_expt1; | 49 | const half_expt2 = exptf - half_expt1; |
| 49 | const scale2 = @bitCast(f32, (0x7f + half_expt2) << 23); | 50 | const scale2 = @bitCast(f32, (0x7f + half_expt2) << 23); |
| 50 | 51 | ||
| 51 | return Complex(f32).init(math.cos(z.im) * exp_x * scale1 * scale2, math.sin(z.im) * exp_x * scale1 * scale2); | 52 | return Complex(f32).init( |
| 53 | math.cos(z.im) * exp_x * scale1 * scale2, | ||
| 54 | math.sin(z.im) * exp_x * scale1 * scale2, | ||
| 55 | ); | ||
| 52 | } | 56 | } |
| 53 | 57 | ||
| 54 | fn frexp_exp64(x: f64, expt: *i32) f64 { | 58 | fn frexp_exp64(x: f64, expt: *i32) f64 { |
| ... | @@ -57,7 +61,7 @@ fn frexp_exp64(x: f64, expt: *i32) f64 { | ... | @@ -57,7 +61,7 @@ fn frexp_exp64(x: f64, expt: *i32) f64 { |
| 57 | 61 | ||
| 58 | const exp_x = math.exp(x - kln2); | 62 | const exp_x = math.exp(x - kln2); |
| 59 | 63 | ||
| 60 | const fx = @bitCast(u64, x); | 64 | const fx = @bitCast(u64, exp_x); |
| 61 | const hx = @intCast(u32, fx >> 32); | 65 | const hx = @intCast(u32, fx >> 32); |
| 62 | const lx = @truncate(u32, fx); | 66 | const lx = @truncate(u32, fx); |
| 63 | 67 | ||
| ... | @@ -73,10 +77,10 @@ fn ldexp_cexp64(z: Complex(f64), expt: i32) Complex(f64) { | ... | @@ -73,10 +77,10 @@ fn ldexp_cexp64(z: Complex(f64), expt: i32) Complex(f64) { |
| 73 | const exptf = @as(i64, expt + ex_expt); | 77 | const exptf = @as(i64, expt + ex_expt); |
| 74 | 78 | ||
| 75 | const half_expt1 = @divTrunc(exptf, 2); | 79 | const half_expt1 = @divTrunc(exptf, 2); |
| 76 | const scale1 = @bitCast(f64, (0x3ff + half_expt1) << 20); | 80 | const scale1 = @bitCast(f64, (0x3ff + half_expt1) << (20 + 32)); |
| 77 | 81 | ||
| 78 | const half_expt2 = exptf - half_expt1; | 82 | const half_expt2 = exptf - half_expt1; |
| 79 | const scale2 = @bitCast(f64, (0x3ff + half_expt2) << 20); | 83 | const scale2 = @bitCast(f64, (0x3ff + half_expt2) << (20 + 32)); |
| 80 | 84 | ||
| 81 | return Complex(f64).init( | 85 | return Complex(f64).init( |
| 82 | math.cos(z.im) * exp_x * scale1 * scale2, | 86 | math.cos(z.im) * exp_x * scale1 * scale2, |