libm/math/
rem_pio2_large.rs

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

libm/math/
rem_pio2_large.rs
