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