1 /* Quad-precision floating point argument reduction.
2 Copyright (C) 1999, 2007, 2009-2011 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/>. */
27 /* Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi */
28 static const int two_over_pi[] = {
29 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
30 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
31 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
32 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
33 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
34 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
35 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
36 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
37 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
38 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
39 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
40 0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6,
41 0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2,
42 0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35,
43 0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30,
44 0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c,
45 0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4,
46 0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770,
47 0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7,
48 0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19,
49 0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522,
50 0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16,
51 0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6,
52 0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e,
53 0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48,
54 0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3,
55 0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf,
56 0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55,
57 0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612,
58 0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929,
59 0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec,
60 0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b,
61 0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c,
62 0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4,
63 0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb,
64 0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc,
65 0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c,
66 0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f,
67 0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5,
68 0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437,
69 0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b,
70 0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea,
71 0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad,
72 0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3,
73 0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3,
74 0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717,
75 0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f,
76 0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61,
77 0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db,
78 0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51,
79 0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0,
80 0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c,
81 0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6,
82 0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc,
83 0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed,
84 0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328,
85 0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d,
86 0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0,
87 0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b,
88 0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4,
89 0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3,
90 0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f,
91 0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad,
92 0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b,
93 0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4,
94 0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761,
95 0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31,
96 0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30,
97 0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262,
98 0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e,
99 0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1,
100 0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c,
101 0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4,
102 0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08,
103 0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196,
104 0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9,
105 0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4,
106 0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc,
107 0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c,
108 0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0,
109 0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c,
110 0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0,
111 0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac,
112 0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22,
113 0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893,
114 0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7,
115 0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5,
116 0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f,
117 0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4,
118 0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf,
119 0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b,
120 0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2,
121 0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138,
122 0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e,
123 0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569,
124 0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34,
125 0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9,
126 0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d,
127 0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f,
128 0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855,
129 0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569,
130 0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b,
131 0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe,
132 0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41,
133 0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49,
134 0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f,
135 0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110,
136 0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8,
137 0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365,
138 0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a,
139 0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270,
140 0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5,
141 0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616,
142 0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b,
143 0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0,
144 0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb,
145 0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a,
146 0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e,
147 0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa,
148 0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5,
149 0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0,
150 0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2,
151 0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886,
152 0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142,
153 0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba,
154 0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4,
155 0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708,
156 0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555,
157 0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3,
158 0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55,
159 0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58,
160 0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5,
161 0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c,
162 0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe,
163 0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b,
164 0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8,
165 0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005,
166 0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7,
167 0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50,
168 0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604,
169 0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643,
170 0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485,
171 0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d,
172 0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6,
173 0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2,
174 0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02,
175 0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3,
176 0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412,
177 0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274,
178 0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755,
179 0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849,
180 0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce,
181 0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5,
182 0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba,
183 0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6,
184 0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d,
188 static const long double c[] = {
189 /* 93 bits of pi/2 */
191 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */
195 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */
198 static int kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
202 ieee754_rem_pio2l (long double x, long double *y)
208 if (x >= -0.78539816339744830961566084581987572104929234984377
209 && x <= 0.78539816339744830961566084581987572104929234984377)
210 /* x in <-pi/4, pi/4> */
217 if (x > 0 && x < 2.35619449019234492884698253745962716314787704953131)
219 /* 113 + 93 bit PI is ok */
222 y[1] = (z - y[0]) - PI_2_1t;
226 if (x < 0 && x > -2.35619449019234492884698253745962716314787704953131)
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--)