libm/math/
rem_pio2_large.rs

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    #[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}
libm/math/rem_pio2_large.rs

libm/math/
rem_pio2_large.rs
