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::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[]     output 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(assert_no_panic, no_panic::no_panic)]
225pub(crate) fn rem_pio2_large(x: &[f64], y: &mut [f64], e0: i32, prec: usize) -> i32 {
226    // FIXME(rust-lang/rust#144518): Inline assembly would cause `no_panic` to fail
227    // on the callers of this function. As a workaround, avoid inlining `floor` here
228    // when implemented with assembly.
229    #[cfg_attr(x86_no_sse, inline(never))]
230    extern "C" fn floor(x: f64) -> f64 {
231        super::floor(x)
232    }
233
234    let x1p24 = f64::from_bits(0x4170000000000000); // 0x1p24 === 2 ^ 24
235    let x1p_24 = f64::from_bits(0x3e70000000000000); // 0x1p_24 === 2 ^ (-24)
236
237    if cfg!(target_pointer_width = "64") {
238        debug_assert!(e0 <= 16360);
239    }
240
241    let nx = x.len();
242
243    let mut fw: f64;
244    let mut n: i32;
245    let mut ih: i32;
246    let mut z: f64;
247    let mut f: [f64; 20] = [0.; 20];
248    let mut fq: [f64; 20] = [0.; 20];
249    let mut q: [f64; 20] = [0.; 20];
250    let mut iq: [i32; 20] = [0; 20];
251
252    /* initialize jk*/
253    let jk = i!(INIT_JK, prec);
254    let jp = jk;
255
256    /* determine jx,jv,q0, note that 3>q0 */
257    let jx = nx - 1;
258    let mut jv = div!(e0 - 3, 24);
259    if jv < 0 {
260        jv = 0;
261    }
262    let mut q0 = e0 - 24 * (jv + 1);
263    let jv = jv as usize;
264
265    /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
266    let mut j = (jv as i32) - (jx as i32);
267    let m = jx + jk;
268    for i in 0..=m {
269        i!(f, i, =, if j < 0 {
270            0.
271        } else {
272            i!(IPIO2, j as usize) as f64
273        });
274        j += 1;
275    }
276
277    /* compute q[0],q[1],...q[jk] */
278    for i in 0..=jk {
279        fw = 0f64;
280        for j in 0..=jx {
281            fw += i!(x, j) * i!(f, jx + i - j);
282        }
283        i!(q, i, =, fw);
284    }
285
286    let mut jz = jk;
287
288    'recompute: loop {
289        /* distill q[] into iq[] reversingly */
290        let mut i = 0i32;
291        z = i!(q, jz);
292        for j in (1..=jz).rev() {
293            fw = (x1p_24 * z) as i32 as f64;
294            i!(iq, i as usize, =, (z - x1p24 * fw) as i32);
295            z = i!(q, j - 1) + fw;
296            i += 1;
297        }
298
299        /* compute n */
300        z = scalbn(z, q0); /* actual value of z */
301        z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */
302        n = z as i32;
303        z -= n as f64;
304        ih = 0;
305        if q0 > 0 {
306            /* need iq[jz-1] to determine n */
307            i = i!(iq, jz - 1) >> (24 - q0);
308            n += i;
309            i!(iq, jz - 1, -=, i << (24 - q0));
310            ih = i!(iq, jz - 1) >> (23 - q0);
311        } else if q0 == 0 {
312            ih = i!(iq, jz - 1) >> 23;
313        } else if z >= 0.5 {
314            ih = 2;
315        }
316
317        if ih > 0 {
318            /* q > 0.5 */
319            n += 1;
320            let mut carry = 0i32;
321            for i in 0..jz {
322                /* compute 1-q */
323                let j = i!(iq, i);
324                if carry == 0 {
325                    if j != 0 {
326                        carry = 1;
327                        i!(iq, i, =, 0x1000000 - j);
328                    }
329                } else {
330                    i!(iq, i, =, 0xffffff - j);
331                }
332            }
333            if q0 > 0 {
334                /* rare case: chance is 1 in 12 */
335                match q0 {
336                    1 => {
337                        i!(iq, jz - 1, &=, 0x7fffff);
338                    }
339                    2 => {
340                        i!(iq, jz - 1, &=, 0x3fffff);
341                    }
342                    _ => {}
343                }
344            }
345            if ih == 2 {
346                z = 1. - z;
347                if carry != 0 {
348                    z -= scalbn(1., q0);
349                }
350            }
351        }
352
353        /* check if recomputation is needed */
354        if z == 0. {
355            let mut j = 0;
356            for i in (jk..=jz - 1).rev() {
357                j |= i!(iq, i);
358            }
359            if j == 0 {
360                /* need recomputation */
361                let mut k = 1;
362                while i!(iq, jk - k, ==, 0) {
363                    k += 1; /* k = no. of terms needed */
364                }
365
366                for i in (jz + 1)..=(jz + k) {
367                    /* add q[jz+1] to q[jz+k] */
368                    i!(f, jx + i, =, i!(IPIO2, jv + i) as f64);
369                    fw = 0f64;
370                    for j in 0..=jx {
371                        fw += i!(x, j) * i!(f, jx + i - j);
372                    }
373                    i!(q, i, =, fw);
374                }
375                jz += k;
376                continue 'recompute;
377            }
378        }
379
380        break;
381    }
382
383    /* chop off zero terms */
384    if z == 0. {
385        jz -= 1;
386        q0 -= 24;
387        while i!(iq, jz) == 0 {
388            jz -= 1;
389            q0 -= 24;
390        }
391    } else {
392        /* break z into 24-bit if necessary */
393        z = scalbn(z, -q0);
394        if z >= x1p24 {
395            fw = (x1p_24 * z) as i32 as f64;
396            i!(iq, jz, =, (z - x1p24 * fw) as i32);
397            jz += 1;
398            q0 += 24;
399            i!(iq, jz, =, fw as i32);
400        } else {
401            i!(iq, jz, =, z as i32);
402        }
403    }
404
405    /* convert integer "bit" chunk to floating-point value */
406    fw = scalbn(1., q0);
407    for i in (0..=jz).rev() {
408        i!(q, i, =, fw * (i!(iq, i) as f64));
409        fw *= x1p_24;
410    }
411
412    /* compute PIo2[0,...,jp]*q[jz,...,0] */
413    for i in (0..=jz).rev() {
414        fw = 0f64;
415        let mut k = 0;
416        while (k <= jp) && (k <= jz - i) {
417            fw += i!(PIO2, k) * i!(q, i + k);
418            k += 1;
419        }
420        i!(fq, jz - i, =, fw);
421    }
422
423    /* compress fq[] into y[] */
424    match prec {
425        0 => {
426            fw = 0f64;
427            for i in (0..=jz).rev() {
428                fw += i!(fq, i);
429            }
430            i!(y, 0, =, if ih == 0 { fw } else { -fw });
431        }
432        1 | 2 => {
433            fw = 0f64;
434            for i in (0..=jz).rev() {
435                fw += i!(fq, i);
436            }
437            i!(y, 0, =, if ih == 0 { fw } else { -fw });
438            fw = i!(fq, 0) - fw;
439            for i in 1..=jz {
440                fw += i!(fq, i);
441            }
442            i!(y, 1, =, if ih == 0 { fw } else { -fw });
443        }
444        3 => {
445            /* painful */
446            for i in (1..=jz).rev() {
447                fw = i!(fq, i - 1) + i!(fq, i);
448                i!(fq, i, +=, i!(fq, i - 1) - fw);
449                i!(fq, i - 1, =, fw);
450            }
451            for i in (2..=jz).rev() {
452                fw = i!(fq, i - 1) + i!(fq, i);
453                i!(fq, i, +=, i!(fq, i - 1) - fw);
454                i!(fq, i - 1, =, fw);
455            }
456            fw = 0f64;
457            for i in (2..=jz).rev() {
458                fw += i!(fq, i);
459            }
460            if ih == 0 {
461                i!(y, 0, =, i!(fq, 0));
462                i!(y, 1, =, i!(fq, 1));
463                i!(y, 2, =, fw);
464            } else {
465                i!(y, 0, =, -i!(fq, 0));
466                i!(y, 1, =, -i!(fq, 1));
467                i!(y, 2, =, -fw);
468            }
469        }
470        #[cfg(debug_assertions)]
471        _ => unreachable!(),
472        #[cfg(not(debug_assertions))]
473        _ => {}
474    }
475    n & 7
476}
libm/math/rem_pio2_large.rs

libm/math/
rem_pio2_large.rs
