zig

fork of https://codeberg.org/ziglang/zig
Log | Files | Refs | README | LICENSE

blob 5c14ec66 (3522B) - Raw

      1 // Ported from musl, which is licensed under the MIT license:
      2 // https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
      3 //
      4 // https://git.musl-libc.org/cgit/musl/tree/src/complex/ctanhf.c
      5 // https://git.musl-libc.org/cgit/musl/tree/src/complex/ctanh.c
      6 
      7 const std = @import("../../std.zig");
      8 const testing = std.testing;
      9 const math = std.math;
     10 const cmath = math.complex;
     11 const Complex = cmath.Complex;
     12 
     13 /// Returns the hyperbolic tangent of z.
     14 pub fn tanh(z: var) @typeOf(z) {
     15     const T = @typeOf(z.re);
     16     return switch (T) {
     17         f32 => tanh32(z),
     18         f64 => tanh64(z),
     19         else => @compileError("tan not implemented for " ++ @typeName(z)),
     20     };
     21 }
     22 
     23 fn tanh32(z: Complex(f32)) Complex(f32) {
     24     const x = z.re;
     25     const y = z.im;
     26 
     27     const hx = @bitCast(u32, x);
     28     const ix = hx & 0x7fffffff;
     29 
     30     if (ix >= 0x7f800000) {
     31         if (ix & 0x7fffff != 0) {
     32             const r = if (y == 0) y else x * y;
     33             return Complex(f32).new(x, r);
     34         }
     35         const xx = @bitCast(f32, hx - 0x40000000);
     36         const r = if (math.isInf(y)) y else math.sin(y) * math.cos(y);
     37         return Complex(f32).new(xx, math.copysign(f32, 0, r));
     38     }
     39 
     40     if (!math.isFinite(y)) {
     41         const r = if (ix != 0) y - y else x;
     42         return Complex(f32).new(r, y - y);
     43     }
     44 
     45     // x >= 11
     46     if (ix >= 0x41300000) {
     47         const exp_mx = math.exp(-math.fabs(x));
     48         return Complex(f32).new(math.copysign(f32, 1, x), 4 * math.sin(y) * math.cos(y) * exp_mx * exp_mx);
     49     }
     50 
     51     // Kahan's algorithm
     52     const t = math.tan(y);
     53     const beta = 1.0 + t * t;
     54     const s = math.sinh(x);
     55     const rho = math.sqrt(1 + s * s);
     56     const den = 1 + beta * s * s;
     57 
     58     return Complex(f32).new((beta * rho * s) / den, t / den);
     59 }
     60 
     61 fn tanh64(z: Complex(f64)) Complex(f64) {
     62     const x = z.re;
     63     const y = z.im;
     64 
     65     const fx = @bitCast(u64, x);
     66     // TODO: zig should allow this conversion implicitly because it can notice that the value necessarily
     67     // fits in range.
     68     const hx = @intCast(u32, fx >> 32);
     69     const lx = @truncate(u32, fx);
     70     const ix = hx & 0x7fffffff;
     71 
     72     if (ix >= 0x7ff00000) {
     73         if ((ix & 0x7fffff) | lx != 0) {
     74             const r = if (y == 0) y else x * y;
     75             return Complex(f64).new(x, r);
     76         }
     77 
     78         const xx = @bitCast(f64, (u64(hx - 0x40000000) << 32) | lx);
     79         const r = if (math.isInf(y)) y else math.sin(y) * math.cos(y);
     80         return Complex(f64).new(xx, math.copysign(f64, 0, r));
     81     }
     82 
     83     if (!math.isFinite(y)) {
     84         const r = if (ix != 0) y - y else x;
     85         return Complex(f64).new(r, y - y);
     86     }
     87 
     88     // x >= 22
     89     if (ix >= 0x40360000) {
     90         const exp_mx = math.exp(-math.fabs(x));
     91         return Complex(f64).new(math.copysign(f64, 1, x), 4 * math.sin(y) * math.cos(y) * exp_mx * exp_mx);
     92     }
     93 
     94     // Kahan's algorithm
     95     const t = math.tan(y);
     96     const beta = 1.0 + t * t;
     97     const s = math.sinh(x);
     98     const rho = math.sqrt(1 + s * s);
     99     const den = 1 + beta * s * s;
    100 
    101     return Complex(f64).new((beta * rho * s) / den, t / den);
    102 }
    103 
    104 const epsilon = 0.0001;
    105 
    106 test "complex.ctanh32" {
    107     const a = Complex(f32).new(5, 3);
    108     const c = tanh(a);
    109 
    110     testing.expect(math.approxEq(f32, c.re, 0.999913, epsilon));
    111     testing.expect(math.approxEq(f32, c.im, -0.000025, epsilon));
    112 }
    113 
    114 test "complex.ctanh64" {
    115     const a = Complex(f64).new(5, 3);
    116     const c = tanh(a);
    117 
    118     testing.expect(math.approxEq(f64, c.re, 0.999913, epsilon));
    119     testing.expect(math.approxEq(f64, c.im, -0.000025, epsilon));
    120 }