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