authorgravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2024-07-31 19:19:27-07:00
committergravatar for noreply@github.comGitHub <noreply@github.com> 2024-07-31 19:19:27-07:00
log059856acfc9f87d723a90af6a4214e128b8cae2e
tree851538dd757a467513f7061db4d2ec9b94a2c628
parent7c5ee3efde6d948205c6f6eaa7ab52bda3715fea
parent843885512dd69e083ec9163e6f57822487b46639
signaturebadge-check Signed by PGP key B5690EEEBB952194

Merge pull request #20878 from tiehuis/std-math-complex-fixes

std.math.complex fixes

19 files changed, 102 insertions(+), 120 deletions(-)

lib/std/math/complex/abs.zig+2-3
......@@ -9,10 +9,9 @@ pub fn abs(z: anytype) @TypeOf(z.re, z.im) {
99 return math.hypot(z.re, z.im);
1010}
1111
12const epsilon = 0.0001;
13
1412test 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
......@@ -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
14const epsilon = 0.0001;
15
1614test 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
......@@ -8,15 +8,18 @@ const Complex = cmath.Complex;
88pub 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
14const 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
1618test 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
......@@ -9,10 +9,9 @@ pub fn arg(z: anytype) @TypeOf(z.re, z.im) {
99 return math.atan2(z.im, z.re);
1010}
1111
12const epsilon = 0.0001;
13
1412test 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
......@@ -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
20const epsilon = 0.0001;
21
2220test 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
......@@ -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
15const epsilon = 0.0001;
16
1715test 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
......@@ -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
3939fn 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
8873fn 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
121const epsilon = 0.0001;
122
12391test 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
131100test 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
......@@ -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
15const epsilon = 0.0001;
16
1715test 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
......@@ -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
......@@ -11,12 +11,11 @@ pub fn cos(z: anytype) Complex(@TypeOf(z.re, z.im)) {
1111 return cmath.cosh(p);
1212}
1313
14const epsilon = 0.0001;
15
1614test 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
......@@ -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
156const epsilon = 0.0001;
157
158156test 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
166165test 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
174test "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
......@@ -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
16const epsilon = 0.0001;
17
1816test 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
......@@ -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
12const epsilon = 0.0001;
13
1412test 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
......@@ -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
......@@ -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
15const epsilon = 0.0001;
16
1715test 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
......@@ -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
155const epsilon = 0.0001;
156
157155test 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
165164test 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}
lib/std/math/complex/sqrt.zig+8-8
......@@ -43,7 +43,7 @@ fn sqrt32(z: Complex(f32)) Complex(f32) {
4343 // sqrt(-inf + i nan) = nan +- inf i
4444 // sqrt(-inf + iy) = 0 + inf i
4545 if (math.signbit(x)) {
46 return Complex(f32).init(@abs(x - y), math.copysign(x, y));
46 return Complex(f32).init(@abs(y - y), math.copysign(x, y));
4747 } else {
4848 return Complex(f32).init(x, math.copysign(y - y, y));
4949 }
......@@ -94,7 +94,7 @@ fn sqrt64(z: Complex(f64)) Complex(f64) {
9494 // sqrt(-inf + i nan) = nan +- inf i
9595 // sqrt(-inf + iy) = 0 + inf i
9696 if (math.signbit(x)) {
97 return Complex(f64).init(@abs(x - y), math.copysign(x, y));
97 return Complex(f64).init(@abs(y - y), math.copysign(x, y));
9898 } else {
9999 return Complex(f64).init(x, math.copysign(y - y, y));
100100 }
......@@ -127,20 +127,20 @@ fn sqrt64(z: Complex(f64)) Complex(f64) {
127127 return result;
128128}
129129
130const epsilon = 0.0001;
131
132130test sqrt32 {
131 const epsilon = math.floatEps(f32);
133132 const a = Complex(f32).init(5, 3);
134133 const c = sqrt(a);
135134
136 try testing.expect(math.approxEqAbs(f32, c.re, 2.327117, epsilon));
137 try testing.expect(math.approxEqAbs(f32, c.im, 0.644574, epsilon));
135 try testing.expectApproxEqAbs(2.3271174, c.re, epsilon);
136 try testing.expectApproxEqAbs(0.6445742, c.im, epsilon);
138137}
139138
140139test sqrt64 {
140 const epsilon = math.floatEps(f64);
141141 const a = Complex(f64).init(5, 3);
142142 const c = sqrt(a);
143143
144 try testing.expect(math.approxEqAbs(f64, c.re, 2.3271175190399496, epsilon));
145 try testing.expect(math.approxEqAbs(f64, c.im, 0.6445742373246469, epsilon));
144 try testing.expectApproxEqAbs(2.3271175190399496, c.re, epsilon);
145 try testing.expectApproxEqAbs(0.6445742373246469, c.im, epsilon);
146146}
lib/std/math/complex/tan.zig+3-4
......@@ -12,12 +12,11 @@ pub fn tan(z: anytype) Complex(@TypeOf(z.re, z.im)) {
1212 return Complex(T).init(r.im, -r.re);
1313}
1414
15const epsilon = 0.0001;
16
1715test tan {
16 const epsilon = math.floatEps(f32);
1817 const a = Complex(f32).init(5, 3);
1918 const c = tan(a);
2019
21 try testing.expect(math.approxEqAbs(f32, c.re, -0.002708233, epsilon));
22 try testing.expect(math.approxEqAbs(f32, c.im, 1.004165, epsilon));
20 try testing.expectApproxEqAbs(-0.002708233, c.re, epsilon);
21 try testing.expectApproxEqAbs(1.0041647, c.im, epsilon);
2322}
lib/std/math/complex/tanh.zig+16-7
......@@ -70,7 +70,7 @@ fn tanh64(z: Complex(f64)) Complex(f64) {
7070 const ix = hx & 0x7fffffff;
7171
7272 if (ix >= 0x7ff00000) {
73 if ((ix & 0x7fffff) | lx != 0) {
73 if ((ix & 0xfffff) | lx != 0) {
7474 const r = if (y == 0) y else x * y;
7575 return Complex(f64).init(x, r);
7676 }
......@@ -101,20 +101,29 @@ fn tanh64(z: Complex(f64)) Complex(f64) {
101101 return Complex(f64).init((beta * rho * s) / den, t / den);
102102}
103103
104const epsilon = 0.0001;
105
106104test tanh32 {
105 const epsilon = math.floatEps(f32);
107106 const a = Complex(f32).init(5, 3);
108107 const c = tanh(a);
109108
110 try testing.expect(math.approxEqAbs(f32, c.re, 0.999913, epsilon));
111 try testing.expect(math.approxEqAbs(f32, c.im, -0.000025, epsilon));
109 try testing.expectApproxEqAbs(0.99991274, c.re, epsilon);
110 try testing.expectApproxEqAbs(-0.00002536878, c.im, epsilon);
112111}
113112
114113test tanh64 {
114 const epsilon = math.floatEps(f64);
115115 const a = Complex(f64).init(5, 3);
116116 const c = tanh(a);
117117
118 try testing.expect(math.approxEqAbs(f64, c.re, 0.999913, epsilon));
119 try testing.expect(math.approxEqAbs(f64, c.im, -0.000025, epsilon));
118 try testing.expectApproxEqAbs(0.9999128201513536, c.re, epsilon);
119 try testing.expectApproxEqAbs(-0.00002536867620767604, c.im, epsilon);
120}
121
122test "tanh64 musl" {
123 const epsilon = math.floatEps(f64);
124 const a = Complex(f64).init(std.math.inf(f64), std.math.inf(f64));
125 const c = tanh(a);
126
127 try testing.expectApproxEqAbs(1, c.re, epsilon);
128 try testing.expectApproxEqAbs(0, c.im, epsilon);
120129}