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, scalbn};
15
16// initial value for jk
17const 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"))]
30const 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")]
42const 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
122const 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)]
225pub(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    if cfg!(target_pointer_width = "64") {
230        debug_assert!(e0 <= 16360);
231    }
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            i!(y, 0, =, if ih == 0 { fw } else { -fw });
430            fw = i!(fq, 0) - fw;
431            for i in 1..=jz {
432                fw += i!(fq, i);
433            }
434            i!(y, 1, =, if ih == 0 { fw } else { -fw });
435        }
436        3 => {
437            /* painful */
438            for i in (1..=jz).rev() {
439                fw = i!(fq, i - 1) + i!(fq, i);
440                i!(fq, i, +=, i!(fq, i - 1) - fw);
441                i!(fq, i - 1, =, fw);
442            }
443            for i in (2..=jz).rev() {
444                fw = i!(fq, i - 1) + i!(fq, i);
445                i!(fq, i, +=, i!(fq, i - 1) - fw);
446                i!(fq, i - 1, =, fw);
447            }
448            fw = 0f64;
449            for i in (2..=jz).rev() {
450                fw += i!(fq, i);
451            }
452            if ih == 0 {
453                i!(y, 0, =, i!(fq, 0));
454                i!(y, 1, =, i!(fq, 1));
455                i!(y, 2, =, fw);
456            } else {
457                i!(y, 0, =, -i!(fq, 0));
458                i!(y, 1, =, -i!(fq, 1));
459                i!(y, 2, =, -fw);
460            }
461        }
462        #[cfg(debug_assertions)]
463        _ => unreachable!(),
464        #[cfg(not(debug_assertions))]
465        _ => {}
466    }
467    n & 7
468}
libm/math/rem_pio2_large.rs

libm/math/
rem_pio2_large.rs
