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 (any(target_pointer_width = "32" , target_pointer_width = "16" ))] |
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)] |
226 | 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 = 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 | } |
471 | |