Skip to content

Commit 059856a

Browse files
andrewrkandrewrk
authored
Merge pull request #20878 from tiehuis/std-math-complex-fixes
std.math.complex fixes
2 parents 7c5ee3e + 8438855 commit 059856a

19 files changed

+102
-120
lines changed

‎lib/std/math/complex/abs.zig

+2-3
Original file line numberDiff line numberDiff line change
@@ -9,10 +9,9 @@ pub fn abs(z: anytype) @TypeOf(z.re, z.im) {
99
return math.hypot(z.re, z.im);
1010
}
1111

12-
const epsilon = 0.0001;
13-
1412
test abs {
13+
const epsilon = math.floatEps(f32);
1514
const a = Complex(f32).init(5, 3);
1615
const c = abs(a);
17-
try testing.expect(math.approxEqAbs(f32, c, 5.83095, epsilon));
16+
try testing.expectApproxEqAbs(5.8309517, c, epsilon);
1817
}

‎lib/std/math/complex/acos.zig

+3-4
Original file line numberDiff line numberDiff line change
@@ -11,12 +11,11 @@ pub fn acos(z: anytype) Complex(@TypeOf(z.re, z.im)) {
1111
return Complex(T).init(@as(T, math.pi) / 2 - q.re, -q.im);
1212
}
1313

14-
const epsilon = 0.0001;
15-
1614
test acos {
15+
const epsilon = math.floatEps(f32);
1716
const a = Complex(f32).init(5, 3);
1817
const c = acos(a);
1918

20-
try testing.expect(math.approxEqAbs(f32, c.re, 0.546975, epsilon));
21-
try testing.expect(math.approxEqAbs(f32, c.im, -2.452914, epsilon));
19+
try testing.expectApproxEqAbs(0.5469737, c.re, epsilon);
20+
try testing.expectApproxEqAbs(-2.4529128, c.im, epsilon);
2221
}

‎lib/std/math/complex/acosh.zig

+8-5
Original file line numberDiff line numberDiff line change
@@ -8,15 +8,18 @@ const Complex = cmath.Complex;
88
pub fn acosh(z: anytype) Complex(@TypeOf(z.re, z.im)) {
99
const T = @TypeOf(z.re, z.im);
1010
const q = cmath.acos(z);
11-
return Complex(T).init(-q.im, q.re);
12-
}
1311

14-
const epsilon = 0.0001;
12+
return if (math.signbit(z.im))
13+
Complex(T).init(q.im, -q.re)
14+
else
15+
Complex(T).init(-q.im, q.re);
16+
}
1517

1618
test acosh {
19+
const epsilon = math.floatEps(f32);
1720
const a = Complex(f32).init(5, 3);
1821
const c = acosh(a);
1922

20-
try testing.expect(math.approxEqAbs(f32, c.re, 2.452914, epsilon));
21-
try testing.expect(math.approxEqAbs(f32, c.im, 0.546975, epsilon));
23+
try testing.expectApproxEqAbs(2.4529128, c.re, epsilon);
24+
try testing.expectApproxEqAbs(0.5469737, c.im, epsilon);
2225
}

‎lib/std/math/complex/arg.zig

+2-3
Original file line numberDiff line numberDiff line change
@@ -9,10 +9,9 @@ pub fn arg(z: anytype) @TypeOf(z.re, z.im) {
99
return math.atan2(z.im, z.re);
1010
}
1111

12-
const epsilon = 0.0001;
13-
1412
test arg {
13+
const epsilon = math.floatEps(f32);
1514
const a = Complex(f32).init(5, 3);
1615
const c = arg(a);
17-
try testing.expect(math.approxEqAbs(f32, c, 0.540420, epsilon));
16+
try testing.expectApproxEqAbs(0.5404195, c, epsilon);
1817
}

‎lib/std/math/complex/asin.zig

+3-4
Original file line numberDiff line numberDiff line change
@@ -17,12 +17,11 @@ pub fn asin(z: anytype) Complex(@TypeOf(z.re, z.im)) {
1717
return Complex(T).init(r.im, -r.re);
1818
}
1919

20-
const epsilon = 0.0001;
21-
2220
test asin {
21+
const epsilon = math.floatEps(f32);
2322
const a = Complex(f32).init(5, 3);
2423
const c = asin(a);
2524

26-
try testing.expect(math.approxEqAbs(f32, c.re, 1.023822, epsilon));
27-
try testing.expect(math.approxEqAbs(f32, c.im, 2.452914, epsilon));
25+
try testing.expectApproxEqAbs(1.0238227, c.re, epsilon);
26+
try testing.expectApproxEqAbs(2.4529128, c.im, epsilon);
2827
}

‎lib/std/math/complex/asinh.zig

+3-4
Original file line numberDiff line numberDiff line change
@@ -12,12 +12,11 @@ pub fn asinh(z: anytype) Complex(@TypeOf(z.re, z.im)) {
1212
return Complex(T).init(r.im, -r.re);
1313
}
1414

15-
const epsilon = 0.0001;
16-
1715
test asinh {
16+
const epsilon = math.floatEps(f32);
1817
const a = Complex(f32).init(5, 3);
1918
const c = asinh(a);
2019

21-
try testing.expect(math.approxEqAbs(f32, c.re, 2.459831, epsilon));
22-
try testing.expect(math.approxEqAbs(f32, c.im, 0.533999, epsilon));
20+
try testing.expectApproxEqAbs(2.4598298, c.re, epsilon);
21+
try testing.expectApproxEqAbs(0.5339993, c.im, epsilon);
2322
}

‎lib/std/math/complex/atan.zig

+10-40
Original file line numberDiff line numberDiff line change
@@ -32,37 +32,22 @@ fn redupif32(x: f32) f32 {
3232
t -= 0.5;
3333
}
3434

35-
const u = @as(f32, @floatFromInt(@as(i32, @intFromFloat(t))));
36-
return ((x - u * DP1) - u * DP2) - t * DP3;
35+
const u: f32 = @trunc(t);
36+
return ((x - u * DP1) - u * DP2) - u * DP3;
3737
}
3838

3939
fn atan32(z: Complex(f32)) Complex(f32) {
40-
const maxnum = 1.0e38;
41-
4240
const x = z.re;
4341
const y = z.im;
4442

45-
if ((x == 0.0) and (y > 1.0)) {
46-
// overflow
47-
return Complex(f32).init(maxnum, maxnum);
48-
}
49-
5043
const x2 = x * x;
5144
var a = 1.0 - x2 - (y * y);
52-
if (a == 0.0) {
53-
// overflow
54-
return Complex(f32).init(maxnum, maxnum);
55-
}
5645

5746
var t = 0.5 * math.atan2(2.0 * x, a);
5847
const w = redupif32(t);
5948

6049
t = y - 1.0;
6150
a = x2 + t * t;
62-
if (a == 0.0) {
63-
// overflow
64-
return Complex(f32).init(maxnum, maxnum);
65-
}
6651

6752
t = y + 1.0;
6853
a = (x2 + (t * t)) / a;
@@ -81,57 +66,42 @@ fn redupif64(x: f64) f64 {
8166
t -= 0.5;
8267
}
8368

84-
const u = @as(f64, @floatFromInt(@as(i64, @intFromFloat(t))));
85-
return ((x - u * DP1) - u * DP2) - t * DP3;
69+
const u: f64 = @trunc(t);
70+
return ((x - u * DP1) - u * DP2) - u * DP3;
8671
}
8772

8873
fn atan64(z: Complex(f64)) Complex(f64) {
89-
const maxnum = 1.0e308;
90-
9174
const x = z.re;
9275
const y = z.im;
9376

94-
if ((x == 0.0) and (y > 1.0)) {
95-
// overflow
96-
return Complex(f64).init(maxnum, maxnum);
97-
}
98-
9977
const x2 = x * x;
10078
var a = 1.0 - x2 - (y * y);
101-
if (a == 0.0) {
102-
// overflow
103-
return Complex(f64).init(maxnum, maxnum);
104-
}
10579

10680
var t = 0.5 * math.atan2(2.0 * x, a);
10781
const w = redupif64(t);
10882

10983
t = y - 1.0;
11084
a = x2 + t * t;
111-
if (a == 0.0) {
112-
// overflow
113-
return Complex(f64).init(maxnum, maxnum);
114-
}
11585

11686
t = y + 1.0;
11787
a = (x2 + (t * t)) / a;
11888
return Complex(f64).init(w, 0.25 * @log(a));
11989
}
12090

121-
const epsilon = 0.0001;
122-
12391
test atan32 {
92+
const epsilon = math.floatEps(f32);
12493
const a = Complex(f32).init(5, 3);
12594
const c = atan(a);
12695

127-
try testing.expect(math.approxEqAbs(f32, c.re, 1.423679, epsilon));
128-
try testing.expect(math.approxEqAbs(f32, c.im, 0.086569, epsilon));
96+
try testing.expectApproxEqAbs(1.423679, c.re, epsilon);
97+
try testing.expectApproxEqAbs(0.086569, c.im, epsilon);
12998
}
13099

131100
test atan64 {
101+
const epsilon = math.floatEps(f64);
132102
const a = Complex(f64).init(5, 3);
133103
const c = atan(a);
134104

135-
try testing.expect(math.approxEqAbs(f64, c.re, 1.423679, epsilon));
136-
try testing.expect(math.approxEqAbs(f64, c.im, 0.086569, epsilon));
105+
try testing.expectApproxEqAbs(1.4236790442393028, c.re, epsilon);
106+
try testing.expectApproxEqAbs(0.08656905917945844, c.im, epsilon);
137107
}

‎lib/std/math/complex/atanh.zig

+3-4
Original file line numberDiff line numberDiff line change
@@ -12,12 +12,11 @@ pub fn atanh(z: anytype) Complex(@TypeOf(z.re, z.im)) {
1212
return Complex(T).init(r.im, -r.re);
1313
}
1414

15-
const epsilon = 0.0001;
16-
1715
test atanh {
16+
const epsilon = math.floatEps(f32);
1817
const a = Complex(f32).init(5, 3);
1918
const c = atanh(a);
2019

21-
try testing.expect(math.approxEqAbs(f32, c.re, 0.146947, epsilon));
22-
try testing.expect(math.approxEqAbs(f32, c.im, 1.480870, epsilon));
20+
try testing.expectApproxEqAbs(0.14694665, c.re, epsilon);
21+
try testing.expectApproxEqAbs(1.4808695, c.im, epsilon);
2322
}

‎lib/std/math/complex/conj.zig

+2-1
Original file line numberDiff line numberDiff line change
@@ -14,5 +14,6 @@ test conj {
1414
const a = Complex(f32).init(5, 3);
1515
const c = a.conjugate();
1616

17-
try testing.expect(c.re == 5 and c.im == -3);
17+
try testing.expectEqual(5, c.re);
18+
try testing.expectEqual(-3, c.im);
1819
}

‎lib/std/math/complex/cos.zig

+3-4
Original file line numberDiff line numberDiff line change
@@ -11,12 +11,11 @@ pub fn cos(z: anytype) Complex(@TypeOf(z.re, z.im)) {
1111
return cmath.cosh(p);
1212
}
1313

14-
const epsilon = 0.0001;
15-
1614
test cos {
15+
const epsilon = math.floatEps(f32);
1716
const a = Complex(f32).init(5, 3);
1817
const c = cos(a);
1918

20-
try testing.expect(math.approxEqAbs(f32, c.re, 2.855815, epsilon));
21-
try testing.expect(math.approxEqAbs(f32, c.im, 9.606383, epsilon));
19+
try testing.expectApproxEqAbs(2.8558152, c.re, epsilon);
20+
try testing.expectApproxEqAbs(9.606383, c.im, epsilon);
2221
}

‎lib/std/math/complex/cosh.zig

+19-10
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@ fn cosh32(z: Complex(f32)) Complex(f32) {
3434

3535
if (ix < 0x7f800000 and iy < 0x7f800000) {
3636
if (iy == 0) {
37-
return Complex(f32).init(math.cosh(x), y);
37+
return Complex(f32).init(math.cosh(x), x * y);
3838
}
3939
// small x: normal case
4040
if (ix < 0x41100000) {
@@ -45,7 +45,7 @@ fn cosh32(z: Complex(f32)) Complex(f32) {
4545
if (ix < 0x42b17218) {
4646
// x < 88.7: exp(|x|) won't overflow
4747
const h = @exp(@abs(x)) * 0.5;
48-
return Complex(f32).init(math.copysign(h, x) * @cos(y), h * @sin(y));
48+
return Complex(f32).init(h * @cos(y), math.copysign(h, x) * @sin(y));
4949
}
5050
// x < 192.7: scale to avoid overflow
5151
else if (ix < 0x4340b1e7) {
@@ -68,7 +68,7 @@ fn cosh32(z: Complex(f32)) Complex(f32) {
6868
if (hx & 0x7fffff == 0) {
6969
return Complex(f32).init(x * x, math.copysign(@as(f32, 0.0), x) * y);
7070
}
71-
return Complex(f32).init(x, math.copysign(@as(f32, 0.0), (x + x) * y));
71+
return Complex(f32).init(x * x, math.copysign(@as(f32, 0.0), (x + x) * y));
7272
}
7373

7474
if (ix < 0x7f800000 and iy >= 0x7f800000) {
@@ -123,7 +123,7 @@ fn cosh64(z: Complex(f64)) Complex(f64) {
123123
}
124124
// x >= 1455: result always overflows
125125
else {
126-
const h = 0x1p1023;
126+
const h = 0x1p1023 * x;
127127
return Complex(f64).init(h * h * @cos(y), h * @sin(y));
128128
}
129129
}
@@ -153,20 +153,29 @@ fn cosh64(z: Complex(f64)) Complex(f64) {
153153
return Complex(f64).init((x * x) * (y - y), (x + x) * (y - y));
154154
}
155155

156-
const epsilon = 0.0001;
157-
158156
test cosh32 {
157+
const epsilon = math.floatEps(f32);
159158
const a = Complex(f32).init(5, 3);
160159
const c = cosh(a);
161160

162-
try testing.expect(math.approxEqAbs(f32, c.re, -73.467300, epsilon));
163-
try testing.expect(math.approxEqAbs(f32, c.im, 10.471557, epsilon));
161+
try testing.expectApproxEqAbs(-73.467300, c.re, epsilon);
162+
try testing.expectApproxEqAbs(10.471557, c.im, epsilon);
164163
}
165164

166165
test cosh64 {
166+
const epsilon = math.floatEps(f64);
167167
const a = Complex(f64).init(5, 3);
168168
const c = cosh(a);
169169

170-
try testing.expect(math.approxEqAbs(f64, c.re, -73.467300, epsilon));
171-
try testing.expect(math.approxEqAbs(f64, c.im, 10.471557, epsilon));
170+
try testing.expectApproxEqAbs(-73.46729221264526, c.re, epsilon);
171+
try testing.expectApproxEqAbs(10.471557674805572, c.im, epsilon);
172+
}
173+
174+
test "cosh64 musl" {
175+
const epsilon = math.floatEps(f64);
176+
const a = Complex(f64).init(7.44648873421389e17, 1.6008058402057622e19);
177+
const c = cosh(a);
178+
179+
try testing.expectApproxEqAbs(std.math.inf(f64), c.re, epsilon);
180+
try testing.expectApproxEqAbs(std.math.inf(f64), c.im, epsilon);
172181
}

‎lib/std/math/complex/log.zig

+3-4
Original file line numberDiff line numberDiff line change
@@ -13,12 +13,11 @@ pub fn log(z: anytype) Complex(@TypeOf(z.re, z.im)) {
1313
return Complex(T).init(@log(r), phi);
1414
}
1515

16-
const epsilon = 0.0001;
17-
1816
test log {
17+
const epsilon = math.floatEps(f32);
1918
const a = Complex(f32).init(5, 3);
2019
const c = log(a);
2120

22-
try testing.expect(math.approxEqAbs(f32, c.re, 1.763180, epsilon));
23-
try testing.expect(math.approxEqAbs(f32, c.im, 0.540419, epsilon));
21+
try testing.expectApproxEqAbs(1.7631803, c.re, epsilon);
22+
try testing.expectApproxEqAbs(0.5404195, c.im, epsilon);
2423
}

‎lib/std/math/complex/pow.zig

+3-4
Original file line numberDiff line numberDiff line change
@@ -9,13 +9,12 @@ pub fn pow(z: anytype, s: anytype) Complex(@TypeOf(z.re, z.im, s.re, s.im)) {
99
return cmath.exp(cmath.log(z).mul(s));
1010
}
1111

12-
const epsilon = 0.0001;
13-
1412
test pow {
13+
const epsilon = math.floatEps(f32);
1514
const a = Complex(f32).init(5, 3);
1615
const b = Complex(f32).init(2.3, -1.3);
1716
const c = pow(a, b);
1817

19-
try testing.expect(math.approxEqAbs(f32, c.re, 58.049110, epsilon));
20-
try testing.expect(math.approxEqAbs(f32, c.im, -101.003433, epsilon));
18+
try testing.expectApproxEqAbs(58.049110, c.re, epsilon);
19+
try testing.expectApproxEqAbs(-101.003433, c.im, epsilon);
2120
}

‎lib/std/math/complex/proj.zig

+2-1
Original file line numberDiff line numberDiff line change
@@ -19,5 +19,6 @@ test proj {
1919
const a = Complex(f32).init(5, 3);
2020
const c = proj(a);
2121

22-
try testing.expect(c.re == 5 and c.im == 3);
22+
try testing.expectEqual(5, c.re);
23+
try testing.expectEqual(3, c.im);
2324
}

‎lib/std/math/complex/sin.zig

+3-4
Original file line numberDiff line numberDiff line change
@@ -12,12 +12,11 @@ pub fn sin(z: anytype) Complex(@TypeOf(z.re, z.im)) {
1212
return Complex(T).init(q.im, -q.re);
1313
}
1414

15-
const epsilon = 0.0001;
16-
1715
test sin {
16+
const epsilon = math.floatEps(f32);
1817
const a = Complex(f32).init(5, 3);
1918
const c = sin(a);
2019

21-
try testing.expect(math.approxEqAbs(f32, c.re, -9.654126, epsilon));
22-
try testing.expect(math.approxEqAbs(f32, c.im, 2.841692, epsilon));
20+
try testing.expectApproxEqAbs(-9.654126, c.re, epsilon);
21+
try testing.expectApproxEqAbs(2.8416924, c.im, epsilon);
2322
}

‎lib/std/math/complex/sinh.zig

+6-6
Original file line numberDiff line numberDiff line change
@@ -152,20 +152,20 @@ fn sinh64(z: Complex(f64)) Complex(f64) {
152152
return Complex(f64).init((x * x) * (y - y), (x + x) * (y - y));
153153
}
154154

155-
const epsilon = 0.0001;
156-
157155
test sinh32 {
156+
const epsilon = math.floatEps(f32);
158157
const a = Complex(f32).init(5, 3);
159158
const c = sinh(a);
160159

161-
try testing.expect(math.approxEqAbs(f32, c.re, -73.460617, epsilon));
162-
try testing.expect(math.approxEqAbs(f32, c.im, 10.472508, epsilon));
160+
try testing.expectApproxEqAbs(-73.460617, c.re, epsilon);
161+
try testing.expectApproxEqAbs(10.472508, c.im, epsilon);
163162
}
164163

165164
test sinh64 {
165+
const epsilon = math.floatEps(f64);
166166
const a = Complex(f64).init(5, 3);
167167
const c = sinh(a);
168168

169-
try testing.expect(math.approxEqAbs(f64, c.re, -73.460617, epsilon));
170-
try testing.expect(math.approxEqAbs(f64, c.im, 10.472508, epsilon));
169+
try testing.expectApproxEqAbs(-73.46062169567367, c.re, epsilon);
170+
try testing.expectApproxEqAbs(10.472508533940392, c.im, epsilon);
171171
}

0 commit comments

Comments
 (0)