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/cexpf.c
  5// https://git.musl-libc.org/cgit/musl/tree/src/complex/cexp.c
  6
  7const std = @import("../../std.zig");
  8const testing = std.testing;
  9const math = std.math;
 10const cmath = math.complex;
 11const Complex = cmath.Complex;
 12
 13const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
 14
 15/// Returns e raised to the power of z (e^z).
 16pub fn exp(z: anytype) Complex(@TypeOf(z.re, z.im)) {
 17    const T = @TypeOf(z.re, z.im);
 18
 19    return switch (T) {
 20        f32 => exp32(z),
 21        f64 => exp64(z),
 22        else => @compileError("exp not implemented for " ++ @typeName(z)),
 23    };
 24}
 25
 26fn exp32(z: Complex(f32)) Complex(f32) {
 27    const exp_overflow = 0x42b17218; // max_exp * ln2 ~= 88.72283955
 28    const cexp_overflow = 0x43400074; // (max_exp - min_denom_exp) * ln2
 29
 30    const x = z.re;
 31    const y = z.im;
 32
 33    const hy = @as(u32, @bitCast(y)) & 0x7fffffff;
 34    // cexp(x + i0) = exp(x) + i0
 35    if (hy == 0) {
 36        return Complex(f32).init(@exp(x), y);
 37    }
 38
 39    const hx = @as(u32, @bitCast(x));
 40    // cexp(0 + iy) = cos(y) + isin(y)
 41    if ((hx & 0x7fffffff) == 0) {
 42        return Complex(f32).init(@cos(y), @sin(y));
 43    }
 44
 45    if (hy >= 0x7f800000) {
 46        // cexp(finite|nan +- i inf|nan) = nan + i nan
 47        if ((hx & 0x7fffffff) != 0x7f800000) {
 48            return Complex(f32).init(y - y, y - y);
 49        } // cexp(-inf +- i inf|nan) = 0 + i0
 50        else if (hx & 0x80000000 != 0) {
 51            return Complex(f32).init(0, 0);
 52        } // cexp(+inf +- i inf|nan) = inf + i nan
 53        else {
 54            return Complex(f32).init(x, y - y);
 55        }
 56    }
 57
 58    // 88.7 <= x <= 192 so must scale
 59    if (hx >= exp_overflow and hx <= cexp_overflow) {
 60        return ldexp_cexp(z, 0);
 61    } // - x < exp_overflow => exp(x) won't overflow (common)
 62    // - x > cexp_overflow, so exp(x) * s overflows for s > 0
 63    // - x = +-inf
 64    // - x = nan
 65    else {
 66        const exp_x = @exp(x);
 67        return Complex(f32).init(exp_x * @cos(y), exp_x * @sin(y));
 68    }
 69}
 70
 71fn exp64(z: Complex(f64)) Complex(f64) {
 72    const exp_overflow = 0x40862e42; // high bits of max_exp * ln2 ~= 710
 73    const cexp_overflow = 0x4096b8e4; // (max_exp - min_denorm_exp) * ln2
 74
 75    const x = z.re;
 76    const y = z.im;
 77
 78    const fy: u64 = @bitCast(y);
 79    const hy: u32 = @intCast((fy >> 32) & 0x7fffffff);
 80    const ly: u32 = @truncate(fy);
 81
 82    // cexp(x + i0) = exp(x) + i0
 83    if (hy | ly == 0) {
 84        return Complex(f64).init(@exp(x), y);
 85    }
 86
 87    const fx: u64 = @bitCast(x);
 88    const hx: u32 = @intCast(fx >> 32);
 89    const lx: u32 = @truncate(fx);
 90
 91    // cexp(0 + iy) = cos(y) + isin(y)
 92    if ((hx & 0x7fffffff) | lx == 0) {
 93        return Complex(f64).init(@cos(y), @sin(y));
 94    }
 95
 96    if (hy >= 0x7ff00000) {
 97        // cexp(finite|nan +- i inf|nan) = nan + i nan
 98        if (lx != 0 or (hx & 0x7fffffff) != 0x7ff00000) {
 99            return Complex(f64).init(y - y, y - y);
100        } // cexp(-inf +- i inf|nan) = 0 + i0
101        else if (hx & 0x80000000 != 0) {
102            return Complex(f64).init(0, 0);
103        } // cexp(+inf +- i inf|nan) = inf + i nan
104        else {
105            return Complex(f64).init(x, y - y);
106        }
107    }
108
109    // 709.7 <= x <= 1454.3 so must scale
110    if (hx >= exp_overflow and hx <= cexp_overflow) {
111        return ldexp_cexp(z, 0);
112    } // - x < exp_overflow => exp(x) won't overflow (common)
113    // - x > cexp_overflow, so exp(x) * s overflows for s > 0
114    // - x = +-inf
115    // - x = nan
116    else {
117        const exp_x = @exp(x);
118        return Complex(f64).init(exp_x * @cos(y), exp_x * @sin(y));
119    }
120}
121
122test exp32 {
123    const tolerance_f32 = @sqrt(math.floatEps(f32));
124
125    {
126        const a = Complex(f32).init(5, 3);
127        const c = exp(a);
128
129        try testing.expectApproxEqRel(@as(f32, -1.46927917e+02), c.re, tolerance_f32);
130        try testing.expectApproxEqRel(@as(f32, 2.0944065e+01), c.im, tolerance_f32);
131    }
132
133    {
134        const a = Complex(f32).init(88.8, 0x1p-149);
135        const c = exp(a);
136
137        try testing.expectApproxEqAbs(math.inf(f32), c.re, tolerance_f32);
138        try testing.expectApproxEqAbs(@as(f32, 5.15088629e-07), c.im, tolerance_f32);
139    }
140}
141
142test exp64 {
143    const tolerance_f64 = @sqrt(math.floatEps(f64));
144
145    {
146        const a = Complex(f64).init(5, 3);
147        const c = exp(a);
148
149        try testing.expectApproxEqRel(@as(f64, -1.469279139083189e+02), c.re, tolerance_f64);
150        try testing.expectApproxEqRel(@as(f64, 2.094406620874596e+01), c.im, tolerance_f64);
151    }
152
153    {
154        const a = Complex(f64).init(709.8, 0x1p-1074);
155        const c = exp(a);
156
157        try testing.expectApproxEqAbs(math.inf(f64), c.re, tolerance_f64);
158        try testing.expectApproxEqAbs(@as(f64, 9.036659362159884e-16), c.im, tolerance_f64);
159    }
160}