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 debug = std.debug;
 9const math = std.math;
10const testing = std.testing;
11const cmath = math.complex;
12const Complex = cmath.Complex;
13
14/// Returns exp(z) scaled to avoid overflow.
15pub fn ldexp_cexp(z: anytype, expt: i32) Complex(@TypeOf(z.re, z.im)) {
16    const T = @TypeOf(z.re, z.im);
17
18    return switch (T) {
19        f32 => ldexp_cexp32(z, expt),
20        f64 => ldexp_cexp64(z, expt),
21        else => unreachable,
22    };
23}
24
25fn frexp_exp32(x: f32, expt: *i32) f32 {
26    const k = 235; // reduction constant
27    const kln2 = 162.88958740; // k * ln2
28
29    const exp_x = @exp(x - kln2);
30    const hx = @as(u32, @bitCast(exp_x));
31    // TODO zig should allow this cast implicitly because it should know the value is in range
32    expt.* = @as(i32, @intCast(hx >> 23)) - (0x7f + 127) + k;
33    return @as(f32, @bitCast((hx & 0x7fffff) | ((0x7f + 127) << 23)));
34}
35
36fn ldexp_cexp32(z: Complex(f32), expt: i32) Complex(f32) {
37    var ex_expt: i32 = undefined;
38    const exp_x = frexp_exp32(z.re, &ex_expt);
39    const exptf = expt + ex_expt;
40
41    const half_expt1 = @divTrunc(exptf, 2);
42    const scale1 = @as(f32, @bitCast((0x7f + half_expt1) << 23));
43
44    const half_expt2 = exptf - half_expt1;
45    const scale2 = @as(f32, @bitCast((0x7f + half_expt2) << 23));
46
47    return Complex(f32).init(
48        @cos(z.im) * exp_x * scale1 * scale2,
49        @sin(z.im) * exp_x * scale1 * scale2,
50    );
51}
52
53fn frexp_exp64(x: f64, expt: *i32) f64 {
54    const k = 1799; // reduction constant
55    const kln2 = 1246.97177782734161156; // k * ln2
56
57    const exp_x = @exp(x - kln2);
58
59    const fx = @as(u64, @bitCast(exp_x));
60    const hx = @as(u32, @intCast(fx >> 32));
61    const lx = @as(u32, @truncate(fx));
62
63    expt.* = @as(i32, @intCast(hx >> 20)) - (0x3ff + 1023) + k;
64
65    const high_word = (hx & 0xfffff) | ((0x3ff + 1023) << 20);
66    return @as(f64, @bitCast((@as(u64, high_word) << 32) | lx));
67}
68
69fn ldexp_cexp64(z: Complex(f64), expt: i32) Complex(f64) {
70    var ex_expt: i32 = undefined;
71    const exp_x = frexp_exp64(z.re, &ex_expt);
72    const exptf = @as(i64, expt + ex_expt);
73
74    const half_expt1 = @divTrunc(exptf, 2);
75    const scale1 = @as(f64, @bitCast((0x3ff + half_expt1) << (20 + 32)));
76
77    const half_expt2 = exptf - half_expt1;
78    const scale2 = @as(f64, @bitCast((0x3ff + half_expt2) << (20 + 32)));
79
80    return Complex(f64).init(
81        @cos(z.im) * exp_x * scale1 * scale2,
82        @sin(z.im) * exp_x * scale1 * scale2,
83    );
84}