1 /* Quad-precision floating point argument reduction.
2 Copyright (C) 1999, 2007 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Jakub Jelinek <jj@ultra.linux.cz>
6 This program is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 3 of the License, or
9 (at your option) any later version.
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with this program. If not, see <http://www.gnu.org/licenses/>. */
26 /* Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi */
27 static const int two_over_pi[] = {
28 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
29 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
30 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
31 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
32 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
33 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
34 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
35 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
36 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
37 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
38 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
39 0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6,
40 0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2,
41 0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35,
42 0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30,
43 0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c,
44 0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4,
45 0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770,
46 0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7,
47 0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19,
48 0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522,
49 0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16,
50 0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6,
51 0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e,
52 0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48,
53 0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3,
54 0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf,
55 0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55,
56 0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612,
57 0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929,
58 0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec,
59 0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b,
60 0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c,
61 0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4,
62 0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb,
63 0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc,
64 0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c,
65 0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f,
66 0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5,
67 0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437,
68 0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b,
69 0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea,
70 0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad,
71 0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3,
72 0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3,
73 0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717,
74 0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f,
75 0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61,
76 0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db,
77 0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51,
78 0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0,
79 0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c,
80 0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6,
81 0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc,
82 0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed,
83 0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328,
84 0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d,
85 0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0,
86 0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b,
87 0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4,
88 0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3,
89 0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f,
90 0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad,
91 0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b,
92 0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4,
93 0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761,
94 0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31,
95 0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30,
96 0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262,
97 0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e,
98 0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1,
99 0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c,
100 0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4,
101 0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08,
102 0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196,
103 0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9,
104 0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4,
105 0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc,
106 0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c,
107 0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0,
108 0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c,
109 0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0,
110 0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac,
111 0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22,
112 0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893,
113 0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7,
114 0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5,
115 0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f,
116 0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4,
117 0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf,
118 0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b,
119 0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2,
120 0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138,
121 0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e,
122 0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569,
123 0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34,
124 0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9,
125 0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d,
126 0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f,
127 0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855,
128 0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569,
129 0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b,
130 0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe,
131 0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41,
132 0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49,
133 0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f,
134 0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110,
135 0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8,
136 0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365,
137 0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a,
138 0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270,
139 0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5,
140 0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616,
141 0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b,
142 0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0,
143 0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb,
144 0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a,
145 0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e,
146 0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa,
147 0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5,
148 0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0,
149 0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2,
150 0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886,
151 0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142,
152 0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba,
153 0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4,
154 0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708,
155 0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555,
156 0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3,
157 0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55,
158 0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58,
159 0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5,
160 0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c,
161 0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe,
162 0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b,
163 0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8,
164 0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005,
165 0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7,
166 0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50,
167 0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604,
168 0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643,
169 0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485,
170 0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d,
171 0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6,
172 0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2,
173 0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02,
174 0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3,
175 0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412,
176 0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274,
177 0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755,
178 0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849,
179 0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce,
180 0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5,
181 0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba,
182 0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6,
183 0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d,
187 static const long double c[] = {
188 /* 93 bits of pi/2 */
190 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */
194 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */
197 static int kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
201 ieee754_rem_pio2l (long double x, long double *y)
207 if (x >= -0.78539816339744830961566084581987572104929234984377
208 && x < 0.78539816339744830961566084581987572104929234984377)
209 /* x in <-pi/4, pi/4> */
216 if (x >= 2.35619449019234492884698253745962716314787704953131
217 && x < 2.35619449019234492884698253745962716314787704953131)
220 /* 113 + 93 bit PI is ok */
223 y[1] = (z - y[0]) - PI_2_1t;
228 /* 113 + 93 bit PI is ok */
231 y[1] = (z - y[0]) + PI_2_1t;
235 if (x + x == x) /* x is ±oo */
242 /* Handle large arguments.
243 We split the 113 bits of the mantissa into 5 24bit integers
244 stored in a double array. */
246 tx[0] = floorl (z *= 16777216.0);
248 tx[1] = floorl (z *= 16777216.0);
250 tx[2] = floorl (z *= 16777216.0);
252 tx[3] = floorl (z *= 16777216.0);
254 tx[4] = floorl (z *= 16777216.0);
256 n = kernel_rem_pio2 (tx, tx + 5, exp - 24, tx[4] ? 5 : 4, 3, two_over_pi);
258 /* The result is now stored in 3 double values, we need to convert it into
259 two long double values. */
260 t = (long double) tx[6] + (long double) tx[7];
261 w = (long double) tx[5];
266 y[1] = t - (y[0] - w);
272 y[1] = -t - (y[0] + w);
277 /* @(#)k_rem_pio2.c 5.1 93/09/24 */
279 * ====================================================
280 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
282 * Developed at SunPro, a Sun Microsystems, Inc. business.
283 * Permission to use, copy, modify, and distribute this
284 * software is freely granted, provided that this notice
286 * ====================================================
289 #if defined(LIBM_SCCS) && !defined(lint)
290 static char rcsid[] =
291 "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $";
295 * kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
296 * double x[],y[]; int e0,nx,prec; int ipio2[];
298 * kernel_rem_pio2 return the last three digits of N with
300 * so that |y| < pi/2.
302 * The method is to compute the integer (mod 8) and fraction parts of
303 * (2/pi)*x without doing the full multiplication. In general we
304 * skip the part of the product that are known to be a huge integer (
305 * more accurately, = 0 mod 8 ). Thus the number of operations are
306 * independent of the exponent of the input.
308 * (2/pi) is represented by an array of 24-bit integers in ipio2[].
311 * x[] The input value (must be positive) is broken into nx
312 * pieces of 24-bit integers in double precision format.
313 * x[i] will be the i-th 24 bit of x. The scaled exponent
314 * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
315 * match x's up to 24 bits.
317 * Example of breaking a double positive z into x[0]+x[1]+x[2]:
325 * y[] ouput result in an array of double precision numbers.
326 * The dimension of y[] is:
330 * 113-bit precision 3
331 * The actual value is the sum of them. Thus for 113-bit
332 * precision, one may have to do something like:
334 * long double t,w,r_head, r_tail;
335 * t = (long double)y[2] + (long double)y[1];
336 * w = (long double)y[0];
338 * r_tail = w - (r_head - t);
340 * e0 The exponent of x[0]
342 * nx dimension of x[]
344 * prec an integer indicating the precision:
347 * 2 64 bits (extended)
351 * integer array, contains the (24*i)-th to (24*i+23)-th
352 * bit of 2/pi after binary point. The corresponding
355 * ipio2[i] * 2^(-24(i+1)).
358 * double scalbn(), floor();
361 * Here is the description of some local variables:
363 * jk jk+1 is the initial number of terms of ipio2[] needed
364 * in the computation. The recommended value is 2,3,4,
365 * 6 for single, double, extended,and quad.
367 * jz local integer variable indicating the number of
368 * terms of ipio2[] used.
372 * jv index for pointing to the suitable ipio2[] for the
373 * computation. In general, we want
374 * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
375 * is an integer. Thus
376 * e0-3-24*jv >= 0 or (e0-3)/24 >= jv
377 * Hence jv = max(0,(e0-3)/24).
379 * jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
381 * q[] double array with integral value, representing the
382 * 24-bits chunk of the product of x and 2/pi.
384 * q0 the corresponding exponent of q[0]. Note that the
385 * exponent for q[i] would be q0-24*i.
387 * PIo2[] double precision array, obtained by cutting pi/2
388 * into 24 bits chunks.
390 * f[] ipio2[] in floating point
392 * iq[] integer array by breaking up q[] in 24-bits chunk.
394 * fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
396 * ih integer. If >0 it indicates q[] is >= 0.5, hence
397 * it also indicates the *sign* of the result.
404 * The hexadecimal values are the intended ones for the following
405 * constants. The decimal values may be used, provided that the
406 * compiler will convert from decimal to binary accurately enough
407 * to produce the hexadecimal values shown.
410 static const int init_jk[] = { 2, 3, 4, 6 }; /* initial value for jk */
411 static const double PIo2[] = {
412 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
413 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
414 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
415 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
416 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
417 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
418 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
419 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
422 static const double zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
423 twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
426 kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
429 int jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
430 double z, fw, f[20], fq[20], q[20];
436 /* determine jx,jv,q0, note that 3>q0 */
441 q0 = e0 - 24 * (jv + 1);
443 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
446 for (i = 0; i <= m; i++, j++)
447 f[i] = (j < 0) ? zero : (double) ipio2[j];
449 /* compute q[0],q[1],...q[jk] */
450 for (i = 0; i <= jk; i++)
452 for (j = 0, fw = 0.0; j <= jx; j++)
453 fw += x[j] * f[jx + i - j];
459 /* distill q[] into iq[] reversingly */
460 for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--)
462 fw = (double) ((int) (twon24 * z));
463 iq[i] = (int) (z - two24 * fw);
468 z = ldexp (z, q0); /* actual value of z */
469 z -= 8.0 * floor (z * 0.125); /* trim off integer >= 8 */
474 { /* need iq[jz-1] to determine n */
475 i = (iq[jz - 1] >> (24 - q0));
477 iq[jz - 1] -= i << (24 - q0);
478 ih = iq[jz - 1] >> (23 - q0);
481 ih = iq[jz - 1] >> 23;
489 for (i = 0; i < jz; i++)
497 iq[i] = 0x1000000 - j;
501 iq[i] = 0xffffff - j;
504 { /* rare case: chance is 1 in 12 */
508 iq[jz - 1] &= 0x7fffff;
511 iq[jz - 1] &= 0x3fffff;
519 z -= ldexp (one, q0);
523 /* check if recomputation is needed */
527 for (i = jz - 1; i >= jk; i--)
530 { /* need recomputation */
531 for (k = 1; iq[jk - k] == 0; k++); /* k = no. of terms needed */
533 for (i = jz + 1; i <= jz + k; i++)
534 { /* add q[jz+1] to q[jz+k] */
535 f[jx + i] = (double) ipio2[jv + i];
536 for (j = 0, fw = 0.0; j <= jx; j++)
537 fw += x[j] * f[jx + i - j];
545 /* chop off zero terms */
557 { /* break z into 24-bit if necessary */
561 fw = (double) ((int) (twon24 * z));
562 iq[jz] = (int) (z - two24 * fw);
571 /* convert integer "bit" chunk to floating-point value */
572 fw = ldexp (one, q0);
573 for (i = jz; i >= 0; i--)
575 q[i] = fw * (double) iq[i];
579 /* compute PIo2[0,...,jp]*q[jz,...,0] */
580 for (i = jz; i >= 0; i--)
582 for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
583 fw += PIo2[k] * q[i + k];
587 /* compress fq[] into y[] */
592 for (i = jz; i >= 0; i--)
594 y[0] = (ih == 0) ? fw : -fw;
599 for (i = jz; i >= 0; i--)
601 y[0] = (ih == 0) ? fw : -fw;
603 for (i = 1; i <= jz; i++)
605 y[1] = (ih == 0) ? fw : -fw;
607 case 3: /* painful */
608 for (i = jz; i > 0; i--)
610 fw = fq[i - 1] + fq[i];
611 fq[i] += fq[i - 1] - fw;
614 for (i = jz; i > 1; i--)
616 fw = fq[i - 1] + fq[i];
617 fq[i] += fq[i - 1] - fw;
620 for (fw = 0.0, i = jz; i >= 2; i--)