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