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}