Meg's Mirrors

zig

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
7	const std = @import("../../std.zig");
8	const debug = std.debug;
9	const math = std.math;
10	const testing = std.testing;
11	const cmath = math.complex;
12	const Complex = cmath.Complex;
13
14	/// Returns exp(z) scaled to avoid overflow.
15	pub 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
25	fn 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
36	fn 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
53	fn 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
69	fn 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	}