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}
477