fn exp64(z: Complex(f64)) Complex(f64)z: Complex(f64)71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120fn exp64(z: Complex(f64)) Complex(f64) {
const exp_overflow = 0x40862e42; // high bits of max_exp * ln2 ~= 710
const cexp_overflow = 0x4096b8e4; // (max_exp - min_denorm_exp) * ln2
const x = z.re;
const y = z.im;
const fy: u64 = @bitCast(y);
const hy: u32 = @intCast((fy >> 32) & 0x7fffffff);
const ly: u32 = @truncate(fy);
// cexp(x + i0) = exp(x) + i0
if (hy | ly == 0) {
return Complex(f64).init(@exp(x), y);
}
const fx: u64 = @bitCast(x);
const hx: u32 = @intCast(fx >> 32);
const lx: u32 = @truncate(fx);
// cexp(0 + iy) = cos(y) + isin(y)
if ((hx & 0x7fffffff) | lx == 0) {
return Complex(f64).init(@cos(y), @sin(y));
}
if (hy >= 0x7ff00000) {
// cexp(finite|nan +- i inf|nan) = nan + i nan
if (lx != 0 or (hx & 0x7fffffff) != 0x7ff00000) {
return Complex(f64).init(y - y, y - y);
} // cexp(-inf +- i inf|nan) = 0 + i0
else if (hx & 0x80000000 != 0) {
return Complex(f64).init(0, 0);
} // cexp(+inf +- i inf|nan) = inf + i nan
else {
return Complex(f64).init(x, y - y);
}
}
// 709.7 <= x <= 1454.3 so must scale
if (hx >= exp_overflow and hx <= cexp_overflow) {
return ldexp_cexp(z, 0);
} // - x < exp_overflow => exp(x) won't overflow (common)
// - x > cexp_overflow, so exp(x) * s overflows for s > 0
// - x = +-inf
// - x = nan
else {
const exp_x = @exp(x);
return Complex(f64).init(exp_x * @cos(y), exp_x * @sin(y));
}
}