rem_pio2_large.zig (20581B) - Raw
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 6 const std = @import("std"); 7 const math = std.math; 8 9 const 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. 22 const 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 141 const 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 /// 254 pub 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 }