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/math/__rem_pio2l.c
5
6const std = @import("std");
7const math = std.math;
8const ld = math.long_double;
9
10const rem_pio2_large = @import("rem_pio2_large.zig").rem_pio2_large;
11
12pub fn rem_pio2l(comptime T: type, x: T, y: *[2]T) i32 {
13 const impl = switch (T) {
14 f80 => struct {
15 const round1: i8 = 22;
16 const round2: i8 = 61;
17 const nx: i8 = 3;
18 const ny: i8 = 2;
19
20 const pio4: T = 0x1.921fb54442d1846ap-1;
21 // 64 bits of 2/pi
22 const invpio2: T = 6.36619772367581343076e-01; // 0xa2f9836e4e44152a.0p-64
23 // first 39 bits of pi/2
24 const pio2_1: f64 = 1.57079632679597125389e+00; // 0x3FF921FB, 0x54444000
25 // pi/2 - pio2_1
26 const pio2_1t: T = -1.07463465549719416346e-12; // -0x973dcb3b399d747f.0p-103
27 // second 39 bits of pi/2
28 const pio2_2: f64 = -1.07463465549783099519e-12; // -0x12e7b967674000.0p-92
29 // pi/2 - (pio2_1+pio2_2)
30 const pio2_2t: T = 6.36831716351095013979e-25; // 0xc51701b839a25205.0p-144
31 // pi/2 - (pio2_1+pio2_2+pio2_3)
32 const pio2_3t: T = -2.75299651904407171810e-37; // -0xbb5bf6c7ddd660ce.0p-185
33 // third 39 bits of pi/2
34 const pio2_3: f64 = 6.36831716351370313614e-25; // 0x18a2e037074000.0p-133
35
36 fn small(x_val: T) bool {
37 const se = ld.signExponent(x_val);
38 const top = ld.mantissaTop(x_val);
39 const lhs = (@as(u32, se & 0x7fff) << 16) | top;
40 const rhs: u32 = ((0x3fff + 25) << 16) | 0x921f >> 1 | 0x8000;
41 return lhs < rhs;
42 }
43
44 fn quobits(v: T) i32 {
45 const q: i32 = @intFromFloat(v);
46 return @intCast(@as(u32, @bitCast(q)) & 0x7fffffff);
47 }
48 },
49 f128 => struct {
50 const round1: i8 = 51;
51 const round2: i8 = 119;
52 const nx: i8 = 5;
53 const ny: i8 = 3;
54
55 const pio4: T = 0x1.921fb54442d18469898cc51701b8p-1;
56 const invpio2: T = 6.3661977236758134307553505349005747e-01;
57 const pio2_1: T = 1.5707963267948966192292994253909555e+00;
58 const pio2_1t: T = 2.0222662487959507323996846200947577e-21;
59 const pio2_2: T = 2.0222662487959507323994779168837751e-21;
60 const pio2_2t: T = 2.0670321098263988236496903051604844e-43;
61 const pio2_3: T = 2.0670321098263988236499468110329591e-43;
62 const pio2_3t: T = -2.5650587247459238361625433492959285e-65;
63
64 fn small(x_val: T) bool {
65 const se = ld.signExponent(x_val);
66 const top = ld.mantissaTop(x_val);
67 const lhs = (@as(u32, se & 0x7fff) << 16) | top;
68 const rhs: u32 = ((0x3fff + 45) << 16) | 0x921f;
69 return lhs < rhs;
70 }
71
72 fn quobits(fn_val: T) i32 {
73 const q: i64 = @intFromFloat(fn_val);
74 return @intCast(@as(u64, @bitCast(q)) & 0x7fffffff);
75 }
76 },
77 else => @compileError("rem_pio2l supports only f80 and f128, got: " ++ @typeName(T)),
78 };
79
80 const x_se = ld.signExponent(x);
81 const ex: i32 = @intCast(x_se & 0x7fff);
82
83 if (impl.small(x)) {
84 // rint(x/(pi/2))
85 const toint: T = 1.5 / math.floatEps(T);
86 var fn_ = x * impl.invpio2 + toint - toint;
87 var n = impl.quobits(fn_);
88 var r = x - fn_ * @as(T, impl.pio2_1);
89 var w = fn_ * impl.pio2_1t; // 1st round good to 102/180 bits
90
91 // Matters with directed rounding.
92 if (r - w < -impl.pio4) {
93 @branchHint(.unlikely);
94 n -= 1;
95 fn_ -= 1;
96 r = x - fn_ * @as(T, impl.pio2_1);
97 w = fn_ * impl.pio2_1t;
98 } else if (r - w > impl.pio4) {
99 @branchHint(.unlikely);
100 n += 1;
101 fn_ += 1;
102 r = x - fn_ * @as(T, impl.pio2_1);
103 w = fn_ * impl.pio2_1t;
104 }
105
106 y[0] = r - w;
107
108 const ey: i32 = @intCast(ld.signExponent(y[0]) & 0x7fff);
109 if (ex - ey > impl.round1) {
110 var t = r;
111 w = fn_ * impl.pio2_2;
112 r = t - w;
113 w = fn_ * impl.pio2_2t - ((t - r) - w);
114 y[0] = r - w;
115 const ey2: i32 = @intCast(ld.signExponent(y[0]) & 0x7fff);
116 if (ex - ey2 > impl.round2) {
117 t = r;
118 w = fn_ * impl.pio2_3;
119 r = t - w;
120 w = fn_ * impl.pio2_3t - ((t - r) - w);
121 y[0] = r - w;
122 }
123 }
124 y[1] = (r - y[0]) - w;
125 return n;
126 }
127
128 // all other (large) arguments
129 if (ex == 0x7fff) { // x is inf or NaN
130 y[0] = x - x;
131 y[1] = y[0];
132 return 0;
133 }
134
135 var z: T = math.scalbn(@abs(x), -math.ilogb(x) + 23);
136 var tx: [impl.nx]f64 = undefined;
137 var ty: [impl.ny]f64 = undefined;
138 var i: usize = 0;
139
140 while (i < impl.nx - 1) : (i += 1) {
141 tx[i] = @floatFromInt(@as(i32, @intFromFloat(z)));
142 z = (z - @as(T, tx[i])) * 0x1p24;
143 }
144
145 tx[i] = @floatCast(z);
146 while (tx[i] == 0.0) {
147 i -= 1;
148 }
149
150 const n = rem_pio2_large(
151 tx[0..(i + 1)],
152 ty[0..impl.ny],
153 ex - 0x3fff - 23,
154 @intCast(i + 1),
155 impl.ny,
156 );
157 var w: f64 = ty[1];
158 if (impl.ny == 3) {
159 w += ty[2];
160 }
161 const r = ty[0] + w;
162 w -= r - ty[0];
163
164 if (x_se >> 15 != 0) {
165 y[0] = -@as(T, r);
166 y[1] = -@as(T, w);
167 return -n;
168 }
169
170 y[0] = @as(T, r);
171 y[1] = @as(T, w);
172 return n;
173}