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_pio2_large.c
5
6const std = @import("std");
7const math = std.math;
8
9const init_jk = [_]i32{ 3, 4, 4, 6 }; // initial value for jk
10
11///
12/// Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
13///
14/// integer array, contains the (24*i)-th to (24*i+23)-th
15/// bit of 2/pi after binary point. The corresponding
16/// floating value is
17///
18/// ipio2[i] * 2^(-24(i+1)).
19///
20/// NB: This table must have at least (e0-3)/24 + jk terms.
21/// For quad precision (e0 <= 16360, jk = 6), this is 686.
22const ipio2 = [_]i32{
23 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
24 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
25 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
26 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
27 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
28 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
29 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
30 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
31 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
32 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
33 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
34
35 0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
36 0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
37 0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
38 0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
39 0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
40 0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
41 0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
42 0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
43 0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
44 0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
45 0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
46 0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
47 0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
48 0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
49 0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
50 0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
51 0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
52 0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
53 0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
54 0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
55 0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
56 0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
57 0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
58 0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
59 0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
60 0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
61 0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
62 0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
63 0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
64 0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
65 0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
66 0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
67 0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
68 0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
69 0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
70 0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
71 0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
72 0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
73 0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
74 0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
75 0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
76 0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
77 0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
78 0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
79 0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
80 0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
81 0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
82 0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
83 0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
84 0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
85 0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
86 0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
87 0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
88 0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
89 0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
90 0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
91 0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
92 0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
93 0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
94 0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
95 0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
96 0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
97 0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
98 0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
99 0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
100 0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
101 0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
102 0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
103 0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
104 0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
105 0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
106 0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
107 0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
108 0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
109 0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
110 0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
111 0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
112 0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
113 0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
114 0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
115 0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
116 0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
117 0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
118 0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
119 0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
120 0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
121 0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
122 0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
123 0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
124 0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
125 0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
126 0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
127 0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
128 0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
129 0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
130 0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
131 0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
132 0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
133 0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
134 0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
135 0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
136 0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
137 0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
138 0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
139};
140
141const PIo2 = [_]f64{
142 1.57079625129699707031e+00, // 0x3FF921FB, 0x40000000
143 7.54978941586159635335e-08, // 0x3E74442D, 0x00000000
144 5.39030252995776476554e-15, // 0x3CF84698, 0x80000000
145 3.28200341580791294123e-22, // 0x3B78CC51, 0x60000000
146 1.27065575308067607349e-29, // 0x39F01B83, 0x80000000
147 1.22933308981111328932e-36, // 0x387A2520, 0x40000000
148 2.73370053816464559624e-44, // 0x36E38222, 0x80000000
149 2.16741683877804819444e-51, // 0x3569F31D, 0x00000000
150};
151
152/// Returns the last three digits of N with y = x - N*pi/2 so that |y| < pi/2.
153///
154/// The method is to compute the integer (mod 8) and fraction parts of
155/// (2/pi)*x without doing the full multiplication. In general we
156/// skip the part of the product that are known to be a huge integer (
157/// more accurately, = 0 mod 8 ). Thus the number of operations are
158/// independent of the exponent of the input.
159///
160/// (2/pi) is represented by an array of 24-bit integers in ipio2[].
161///
162/// Input parameters:
163/// x[] The input value (must be positive) is broken into nx
164/// pieces of 24-bit integers in double precision format.
165/// x[i] will be the i-th 24 bit of x. The scaled exponent
166/// of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
167/// match x's up to 24 bits.
168///
169/// Example of breaking a double positive z into x[0]+x[1]+x[2]:
170/// e0 = ilogb(z)-23
171/// z = scalbn(z,-e0)
172/// for i = 0,1,2
173/// x[i] = floor(z)
174/// z = (z-x[i])*2**24
175///
176///
177/// y[] output result in an array of double precision numbers.
178/// The dimension of y[] is:
179/// 24-bit precision 1
180/// 53-bit precision 2
181/// 64-bit precision 2
182/// 113-bit precision 3
183/// The actual value is the sum of them. Thus for 113-bit
184/// precision, one may have to do something like:
185///
186/// long double t,w,r_head, r_tail;
187/// t = (long double)y[2] + (long double)y[1];
188/// w = (long double)y[0];
189/// r_head = t+w;
190/// r_tail = w - (r_head - t);
191///
192/// e0 The exponent of x[0]. Must be <= 16360 or you need to
193/// expand the ipio2 table.
194///
195/// nx dimension of x[]
196///
197/// prec an integer indicating the precision:
198/// 0 24 bits (single)
199/// 1 53 bits (double)
200/// 2 64 bits (extended)
201/// 3 113 bits (quad)
202///
203/// Here is the description of some local variables:
204///
205/// jk jk+1 is the initial number of terms of ipio2[] needed
206/// in the computation. The minimum and recommended value
207/// for jk is 3,4,4,6 for single, double, extended, and quad.
208/// jk+1 must be 2 larger than you might expect so that our
209/// recomputation test works. (Up to 24 bits in the integer
210/// part (the 24 bits of it that we compute) and 23 bits in
211/// the fraction part may be lost to cancelation before we
212/// recompute.)
213///
214/// jz local integer variable indicating the number of
215/// terms of ipio2[] used.
216///
217/// jx nx - 1
218///
219/// jv index for pointing to the suitable ipio2[] for the
220/// computation. In general, we want
221/// ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
222/// is an integer. Thus
223/// e0-3-24*jv >= 0 or (e0-3)/24 >= jv
224/// Hence jv = max(0,(e0-3)/24).
225///
226/// jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
227///
228/// q[] double array with integral value, representing the
229/// 24-bits chunk of the product of x and 2/pi.
230///
231/// q0 the corresponding exponent of q[0]. Note that the
232/// exponent for q[i] would be q0-24*i.
233///
234/// PIo2[] double precision array, obtained by cutting pi/2
235/// into 24 bits chunks.
236///
237/// f[] ipio2[] in floating point
238///
239/// iq[] integer array by breaking up q[] in 24-bits chunk.
240///
241/// fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
242///
243/// ih integer. If >0 it indicates q[] is >= 0.5, hence
244/// it also indicates the *sign* of the result.
245///
246///
247///
248/// Constants:
249/// The hexadecimal values are the intended ones for the following
250/// constants. The decimal values may be used, provided that the
251/// compiler will convert from decimal to binary accurately enough
252/// to produce the hexadecimal values shown.
253///
254pub fn rem_pio2_large(x: []const f64, y: []f64, e0: i32, nx: i32, prec: usize) i32 {
255 var jz: i32 = undefined;
256 var jx: i32 = undefined;
257 var jv: i32 = undefined;
258 var jp: i32 = undefined;
259 var jk: i32 = undefined;
260 var carry: i32 = undefined;
261 var n: i32 = undefined;
262 var iq: [20]i32 = undefined;
263 var i: i32 = undefined;
264 var j: i32 = undefined;
265 var k: i32 = undefined;
266 var m: i32 = undefined;
267 var q0: i32 = undefined;
268 var ih: i32 = undefined;
269
270 var z: f64 = undefined;
271 var fw: f64 = undefined;
272 var f: [20]f64 = undefined;
273 var fq: [20]f64 = undefined;
274 var q: [20]f64 = undefined;
275
276 // initialize jk
277 jk = init_jk[prec];
278 jp = jk;
279
280 // determine jx,jv,q0, note that 3>q0
281 jx = nx - 1;
282 jv = @divFloor(e0 - 3, 24);
283 if (jv < 0) jv = 0;
284 q0 = e0 - 24 * (jv + 1);
285
286 // set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk]
287 j = jv - jx;
288 m = jx + jk;
289 i = 0;
290 while (i <= m) : ({
291 i += 1;
292 j += 1;
293 }) {
294 f[@intCast(i)] = if (j < 0) 0.0 else @floatFromInt(ipio2[@intCast(j)]);
295 }
296
297 // compute q[0],q[1],...q[jk]
298 i = 0;
299 while (i <= jk) : (i += 1) {
300 j = 0;
301 fw = 0;
302 while (j <= jx) : (j += 1) {
303 fw += x[@intCast(j)] * f[@intCast(jx + i - j)];
304 }
305 q[@intCast(i)] = fw;
306 }
307
308 jz = jk;
309
310 // This is to handle a non-trivial goto translation from C.
311 // An unconditional return statement is found at the end of this loop.
312 recompute: while (true) {
313 // distill q[] into iq[] reversingly
314 i = 0;
315 j = jz;
316 z = q[@intCast(jz)];
317 while (j > 0) : ({
318 i += 1;
319 j -= 1;
320 }) {
321 fw = @floatFromInt(@as(i32, @intFromFloat(0x1p-24 * z)));
322 iq[@intCast(i)] = @intFromFloat(z - 0x1p24 * fw);
323 z = q[@intCast(j - 1)] + fw;
324 }
325
326 // compute n
327 z = math.scalbn(z, q0); // actual value of z
328 z -= 8.0 * @floor(z * 0.125); // trim off integer >= 8
329 n = @intFromFloat(z);
330 z -= @floatFromInt(n);
331 ih = 0;
332 if (q0 > 0) { // need iq[jz-1] to determine n
333 i = iq[@intCast(jz - 1)] >> @intCast(24 - q0);
334 n += i;
335 iq[@intCast(jz - 1)] -= i << @intCast(24 - q0);
336 ih = iq[@intCast(jz - 1)] >> @intCast(23 - q0);
337 } else if (q0 == 0) {
338 ih = iq[@intCast(jz - 1)] >> 23;
339 } else if (z >= 0.5) {
340 ih = 2;
341 }
342
343 if (ih > 0) { // q > 0.5
344 n += 1;
345 carry = 0;
346 i = 0;
347 while (i < jz) : (i += 1) { // compute 1-q
348 j = iq[@intCast(i)];
349 if (carry == 0) {
350 if (j != 0) {
351 carry = 1;
352 iq[@intCast(i)] = 0x1000000 - j;
353 }
354 } else {
355 iq[@intCast(i)] = 0xffffff - j;
356 }
357 }
358 if (q0 > 0) { // rare case: chance is 1 in 12
359 @branchHint(.unlikely);
360 switch (q0) {
361 1 => iq[@intCast(jz - 1)] &= 0x7fffff,
362 2 => iq[@intCast(jz - 1)] &= 0x3fffff,
363 else => unreachable,
364 }
365 }
366 if (ih == 2) {
367 z = 1.0 - z;
368 if (carry != 0) {
369 z -= math.scalbn(@as(f64, 1.0), q0);
370 }
371 }
372 }
373
374 // check if recomputation is needed
375 if (z == 0.0) {
376 j = 0;
377 i = jz - 1;
378 while (i >= jk) : (i -= 1) {
379 j |= iq[@intCast(i)];
380 }
381
382 if (j == 0) { // need recomputation
383 k = 1;
384 while (iq[@intCast(jk - k)] == 0) : (k += 1) {
385 // k = no. of terms needed
386 }
387
388 i = jz + 1;
389 while (i <= jz + k) : (i += 1) { // add q[jz+1] to q[jz+k]
390 f[@intCast(jx + i)] = @floatFromInt(ipio2[@intCast(jv + i)]);
391 j = 0;
392 fw = 0;
393 while (j <= jx) : (j += 1) {
394 fw += x[@intCast(j)] * f[@intCast(jx + i - j)];
395 }
396 q[@intCast(i)] = fw;
397 }
398 jz += k;
399 continue :recompute; // mimic goto recompute
400 }
401 }
402
403 // chop off zero terms
404 if (z == 0.0) {
405 jz -= 1;
406 q0 -= 24;
407 while (iq[@intCast(jz)] == 0) {
408 jz -= 1;
409 q0 -= 24;
410 }
411 } else { // break z into 24-bit if necessary
412 z = math.scalbn(z, -q0);
413 if (z >= 0x1p24) {
414 fw = @floatFromInt(@as(i32, @intFromFloat(0x1p-24 * z)));
415 iq[@intCast(jz)] = @intFromFloat(z - 0x1p24 * fw);
416 jz += 1;
417 q0 += 24;
418 iq[@intCast(jz)] = @intFromFloat(fw);
419 } else {
420 iq[@intCast(jz)] = @intFromFloat(z);
421 }
422 }
423
424 // convert integer "bit" chunk to floating-point value
425 fw = math.scalbn(@as(f64, 1.0), q0);
426 i = jz;
427 while (i >= 0) : (i -= 1) {
428 q[@intCast(i)] = fw * @as(f64, @floatFromInt(iq[@intCast(i)]));
429 fw *= 0x1p-24;
430 }
431
432 // compute PIo2[0,...,jp]*q[jz,...,0]
433 i = jz;
434 while (i >= 0) : (i -= 1) {
435 fw = 0;
436 k = 0;
437 while (k <= jp and k <= jz - i) : (k += 1) {
438 fw += PIo2[@intCast(k)] * q[@intCast(i + k)];
439 }
440 fq[@intCast(jz - i)] = fw;
441 }
442
443 // compress fq[] into y[]
444 switch (prec) {
445 0 => {
446 fw = 0.0;
447 i = jz;
448 while (i >= 0) : (i -= 1) {
449 fw += fq[@intCast(i)];
450 }
451 y[0] = if (ih == 0) fw else -fw;
452 },
453
454 1, 2 => {
455 fw = 0.0;
456 i = jz;
457 while (i >= 0) : (i -= 1) {
458 fw += fq[@intCast(i)];
459 }
460 // TODO: drop excess precision here once double_t is used
461 fw = fw;
462 y[0] = if (ih == 0) fw else -fw;
463 fw = fq[0] - fw;
464 i = 1;
465 while (i <= jz) : (i += 1) {
466 fw += fq[@intCast(i)];
467 }
468 y[1] = if (ih == 0) fw else -fw;
469 },
470 3 => { // painful
471 i = jz;
472 while (i > 0) : (i -= 1) {
473 fw = fq[@intCast(i - 1)] + fq[@intCast(i)];
474 fq[@intCast(i)] += fq[@intCast(i - 1)] - fw;
475 fq[@intCast(i - 1)] = fw;
476 }
477 i = jz;
478 while (i > 1) : (i -= 1) {
479 fw = fq[@intCast(i - 1)] + fq[@intCast(i)];
480 fq[@intCast(i)] += fq[@intCast(i - 1)] - fw;
481 fq[@intCast(i - 1)] = fw;
482 }
483 fw = 0;
484 i = jz;
485 while (i >= 2) : (i -= 1) {
486 fw += fq[@intCast(i)];
487 }
488 if (ih == 0) {
489 y[0] = fq[0];
490 y[1] = fq[1];
491 y[2] = fw;
492 } else {
493 y[0] = -fq[0];
494 y[1] = -fq[1];
495 y[2] = -fw;
496 }
497 },
498 else => unreachable,
499 }
500
501 return n & 7;
502 }
503}