Commit 2bd7af63d7

Marc Tiehuis <marc@tiehu.is>
2024-07-29 06:26:09
std.math.complex: fix acosh/atan/cosh/sqrt
Some of these are upstream changes since the original port, others are translation errors.
1 parent f219286
Changed files (4)
lib/std/math/complex/acosh.zig
@@ -8,7 +8,11 @@ const Complex = cmath.Complex;
 pub fn acosh(z: anytype) Complex(@TypeOf(z.re, z.im)) {
     const T = @TypeOf(z.re, z.im);
     const q = cmath.acos(z);
-    return Complex(T).init(-q.im, q.re);
+
+    return if (math.signbit(z.im))
+        Complex(T).init(q.im, -q.re)
+    else
+        Complex(T).init(-q.im, q.re);
 }
 
 const epsilon = 0.0001;
lib/std/math/complex/atan.zig
@@ -32,37 +32,22 @@ fn redupif32(x: f32) f32 {
         t -= 0.5;
     }
 
-    const u = @as(f32, @floatFromInt(@as(i32, @intFromFloat(t))));
-    return ((x - u * DP1) - u * DP2) - t * DP3;
+    const u: f32 = @trunc(t);
+    return ((x - u * DP1) - u * DP2) - u * DP3;
 }
 
 fn atan32(z: Complex(f32)) Complex(f32) {
-    const maxnum = 1.0e38;
-
     const x = z.re;
     const y = z.im;
 
-    if ((x == 0.0) and (y > 1.0)) {
-        // overflow
-        return Complex(f32).init(maxnum, maxnum);
-    }
-
     const x2 = x * x;
     var a = 1.0 - x2 - (y * y);
-    if (a == 0.0) {
-        // overflow
-        return Complex(f32).init(maxnum, maxnum);
-    }
 
     var t = 0.5 * math.atan2(2.0 * x, a);
     const w = redupif32(t);
 
     t = y - 1.0;
     a = x2 + t * t;
-    if (a == 0.0) {
-        // overflow
-        return Complex(f32).init(maxnum, maxnum);
-    }
 
     t = y + 1.0;
     a = (x2 + (t * t)) / a;
@@ -81,37 +66,22 @@ fn redupif64(x: f64) f64 {
         t -= 0.5;
     }
 
-    const u = @as(f64, @floatFromInt(@as(i64, @intFromFloat(t))));
-    return ((x - u * DP1) - u * DP2) - t * DP3;
+    const u: f64 = @trunc(t);
+    return ((x - u * DP1) - u * DP2) - u * DP3;
 }
 
 fn atan64(z: Complex(f64)) Complex(f64) {
-    const maxnum = 1.0e308;
-
     const x = z.re;
     const y = z.im;
 
-    if ((x == 0.0) and (y > 1.0)) {
-        // overflow
-        return Complex(f64).init(maxnum, maxnum);
-    }
-
     const x2 = x * x;
     var a = 1.0 - x2 - (y * y);
-    if (a == 0.0) {
-        // overflow
-        return Complex(f64).init(maxnum, maxnum);
-    }
 
     var t = 0.5 * math.atan2(2.0 * x, a);
     const w = redupif64(t);
 
     t = y - 1.0;
     a = x2 + t * t;
-    if (a == 0.0) {
-        // overflow
-        return Complex(f64).init(maxnum, maxnum);
-    }
 
     t = y + 1.0;
     a = (x2 + (t * t)) / a;
lib/std/math/complex/cosh.zig
@@ -34,7 +34,7 @@ fn cosh32(z: Complex(f32)) Complex(f32) {
 
     if (ix < 0x7f800000 and iy < 0x7f800000) {
         if (iy == 0) {
-            return Complex(f32).init(math.cosh(x), y);
+            return Complex(f32).init(math.cosh(x), x * y);
         }
         // small x: normal case
         if (ix < 0x41100000) {
@@ -45,7 +45,7 @@ fn cosh32(z: Complex(f32)) Complex(f32) {
         if (ix < 0x42b17218) {
             // x < 88.7: exp(|x|) won't overflow
             const h = @exp(@abs(x)) * 0.5;
-            return Complex(f32).init(math.copysign(h, x) * @cos(y), h * @sin(y));
+            return Complex(f32).init(h * @cos(y), math.copysign(h, x) * @sin(y));
         }
         // x < 192.7: scale to avoid overflow
         else if (ix < 0x4340b1e7) {
@@ -68,7 +68,7 @@ fn cosh32(z: Complex(f32)) Complex(f32) {
         if (hx & 0x7fffff == 0) {
             return Complex(f32).init(x * x, math.copysign(@as(f32, 0.0), x) * y);
         }
-        return Complex(f32).init(x, math.copysign(@as(f32, 0.0), (x + x) * y));
+        return Complex(f32).init(x * x, math.copysign(@as(f32, 0.0), (x + x) * y));
     }
 
     if (ix < 0x7f800000 and iy >= 0x7f800000) {
lib/std/math/complex/sqrt.zig
@@ -43,7 +43,7 @@ fn sqrt32(z: Complex(f32)) Complex(f32) {
         // sqrt(-inf + i nan)   = nan +- inf i
         // sqrt(-inf + iy)      = 0 + inf i
         if (math.signbit(x)) {
-            return Complex(f32).init(@abs(x - y), math.copysign(x, y));
+            return Complex(f32).init(@abs(y - y), math.copysign(x, y));
         } else {
             return Complex(f32).init(x, math.copysign(y - y, y));
         }
@@ -94,7 +94,7 @@ fn sqrt64(z: Complex(f64)) Complex(f64) {
         // sqrt(-inf + i nan)   = nan +- inf i
         // sqrt(-inf + iy)      = 0 + inf i
         if (math.signbit(x)) {
-            return Complex(f64).init(@abs(x - y), math.copysign(x, y));
+            return Complex(f64).init(@abs(y - y), math.copysign(x, y));
         } else {
             return Complex(f64).init(x, math.copysign(y - y, y));
         }