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
14 use super::floor;
15 use super::scalbn;
16
17 // initial value for jk
18 const 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(target_pointer_width = "32")]
31 const 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")]
43 const 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
123 const 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)]
rem_pio2_large(x: &[f64], y: &mut [f64], e0: i32, prec: usize) -> i32226 pub(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 = 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 = (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 }
471