This source file includes following definitions.
- ieee754_rem_pio2l
- kernel_rem_pio2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19 #include <config.h>
20
21
22 #include "trigl.h"
23
24 #include <float.h>
25 #include <math.h>
26
27
28
29
30
31 static const int two_over_pi[] = {
32 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
33 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
34 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
35 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
36 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
37 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
38 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
39 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
40 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
41 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
42 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
43 0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6,
44 0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2,
45 0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35,
46 0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30,
47 0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c,
48 0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4,
49 0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770,
50 0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7,
51 0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19,
52 0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522,
53 0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16,
54 0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6,
55 0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e,
56 0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48,
57 0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3,
58 0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf,
59 0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55,
60 0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612,
61 0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929,
62 0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec,
63 0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b,
64 0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c,
65 0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4,
66 0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb,
67 0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc,
68 0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c,
69 0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f,
70 0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5,
71 0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437,
72 0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b,
73 0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea,
74 0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad,
75 0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3,
76 0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3,
77 0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717,
78 0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f,
79 0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61,
80 0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db,
81 0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51,
82 0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0,
83 0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c,
84 0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6,
85 0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc,
86 0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed,
87 0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328,
88 0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d,
89 0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0,
90 0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b,
91 0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4,
92 0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3,
93 0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f,
94 0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad,
95 0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b,
96 0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4,
97 0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761,
98 0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31,
99 0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30,
100 0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262,
101 0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e,
102 0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1,
103 0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c,
104 0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4,
105 0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08,
106 0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196,
107 0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9,
108 0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4,
109 0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc,
110 0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c,
111 0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0,
112 0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c,
113 0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0,
114 0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac,
115 0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22,
116 0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893,
117 0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7,
118 0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5,
119 0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f,
120 0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4,
121 0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf,
122 0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b,
123 0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2,
124 0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138,
125 0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e,
126 0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569,
127 0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34,
128 0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9,
129 0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d,
130 0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f,
131 0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855,
132 0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569,
133 0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b,
134 0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe,
135 0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41,
136 0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49,
137 0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f,
138 0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110,
139 0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8,
140 0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365,
141 0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a,
142 0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270,
143 0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5,
144 0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616,
145 0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b,
146 0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0,
147 0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb,
148 0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a,
149 0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e,
150 0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa,
151 0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5,
152 0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0,
153 0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2,
154 0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886,
155 0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142,
156 0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba,
157 0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4,
158 0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708,
159 0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555,
160 0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3,
161 0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55,
162 0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58,
163 0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5,
164 0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c,
165 0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe,
166 0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b,
167 0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8,
168 0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005,
169 0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7,
170 0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50,
171 0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604,
172 0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643,
173 0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485,
174 0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d,
175 0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6,
176 0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2,
177 0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02,
178 0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3,
179 0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412,
180 0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274,
181 0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755,
182 0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849,
183 0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce,
184 0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5,
185 0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba,
186 0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6,
187 0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d,
188 0x7b7b89, 0x483d38,
189 };
190
191 static const long double c[] = {
192
193 #define PI_2_1 c[0]
194 1.57079632679489661923132169155131424e+00L,
195
196
197 #define PI_2_1t c[1]
198 8.84372056613570112025531863263659260e-29L,
199 };
200
201 static int kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
202 const int *ipio2);
203
204 int
205 ieee754_rem_pio2l (long double x, long double *y)
206 {
207 long double z, w, t;
208 double tx[8];
209 int exp, n;
210
211 if (x >= -0.78539816339744830961566084581987572104929234984377
212 && x <= 0.78539816339744830961566084581987572104929234984377)
213
214 {
215 y[0] = x;
216 y[1] = 0;
217 return 0;
218 }
219
220 if (x > 0 && x < 2.35619449019234492884698253745962716314787704953131)
221 {
222
223 z = x - PI_2_1;
224 y[0] = z - PI_2_1t;
225 y[1] = (z - y[0]) - PI_2_1t;
226 return 1;
227 }
228
229 if (x < 0 && x > -2.35619449019234492884698253745962716314787704953131)
230 {
231
232 z = x + PI_2_1;
233 y[0] = z + PI_2_1t;
234 y[1] = (z - y[0]) + PI_2_1t;
235 return -1;
236 }
237
238 if (x + x == x)
239 {
240 y[0] = x - x;
241 y[1] = y[0];
242 return 0;
243 }
244
245
246
247
248 z = frexp (x, &exp);
249 tx[0] = floorl (z *= 16777216.0);
250 z -= tx[0];
251 tx[1] = floorl (z *= 16777216.0);
252 z -= tx[1];
253 tx[2] = floorl (z *= 16777216.0);
254 z -= tx[2];
255 tx[3] = floorl (z *= 16777216.0);
256 z -= tx[3];
257 tx[4] = floorl (z *= 16777216.0);
258
259 n = kernel_rem_pio2 (tx, tx + 5, exp - 24, tx[4] ? 5 : 4, 3, two_over_pi);
260
261
262
263 t = (long double) tx[6] + (long double) tx[7];
264 w = (long double) tx[5];
265
266 if (x > 0)
267 {
268 y[0] = w + t;
269 y[1] = t - (y[0] - w);
270 return n;
271 }
272 else
273 {
274 y[0] = -(w + t);
275 y[1] = -t - (y[0] + w);
276 return -n;
277 }
278 }
279
280
281
282
283
284
285
286
287
288
289
290
291
292 #if defined LIBM_SCCS && !defined GCC_LINT && !defined lint
293 static char rcsid[] =
294 "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $";
295 #endif
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413 static const int init_jk[] = { 2, 3, 4, 6 };
414 static const double PIo2[] = {
415 1.57079625129699707031e+00,
416 7.54978941586159635335e-08,
417 5.39030252995776476554e-15,
418 3.28200341580791294123e-22,
419 1.27065575308067607349e-29,
420 1.22933308981111328932e-36,
421 2.73370053816464559624e-44,
422 2.16741683877804819444e-51,
423 };
424
425 static const double zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07,
426 twon24 = 5.96046447753906250000e-08;
427
428 static int
429 kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
430 const int *ipio2)
431 {
432 int jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
433 double z, fw, f[20], fq[20], q[20];
434
435
436 jk = init_jk[prec];
437 jp = jk;
438
439
440 jx = nx - 1;
441 jv = (e0 - 3) / 24;
442 if (jv < 0)
443 jv = 0;
444 q0 = e0 - 24 * (jv + 1);
445
446
447 j = jv - jx;
448 m = jx + jk;
449 for (i = 0; i <= m; i++, j++)
450 f[i] = (j < 0) ? zero : (double) ipio2[j];
451
452
453 for (i = 0; i <= jk; i++)
454 {
455 for (j = 0, fw = 0.0; j <= jx; j++)
456 fw += x[j] * f[jx + i - j];
457 q[i] = fw;
458 }
459
460 jz = jk;
461 recompute:
462
463 for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--)
464 {
465 fw = (double) ((int) (twon24 * z));
466 iq[i] = (int) (z - two24 * fw);
467 z = q[j - 1] + fw;
468 }
469
470
471 z = ldexp (z, q0);
472 z -= 8.0 * floor (z * 0.125);
473 n = (int) z;
474 z -= (double) n;
475 ih = 0;
476 if (q0 > 0)
477 {
478 i = (iq[jz - 1] >> (24 - q0));
479 n += i;
480 iq[jz - 1] -= i << (24 - q0);
481 ih = iq[jz - 1] >> (23 - q0);
482 }
483 else if (q0 == 0)
484 ih = iq[jz - 1] >> 23;
485 else if (z >= 0.5)
486 ih = 2;
487
488 if (ih > 0)
489 {
490 n += 1;
491 carry = 0;
492 for (i = 0; i < jz; i++)
493 {
494 j = iq[i];
495 if (carry == 0)
496 {
497 if (j != 0)
498 {
499 carry = 1;
500 iq[i] = 0x1000000 - j;
501 }
502 }
503 else
504 iq[i] = 0xffffff - j;
505 }
506 if (q0 > 0)
507 {
508 switch (q0)
509 {
510 case 1:
511 iq[jz - 1] &= 0x7fffff;
512 break;
513 case 2:
514 iq[jz - 1] &= 0x3fffff;
515 break;
516 }
517 }
518 if (ih == 2)
519 {
520 z = one - z;
521 if (carry != 0)
522 z -= ldexp (one, q0);
523 }
524 }
525
526
527 if (z == zero)
528 {
529 j = 0;
530 for (i = jz - 1; i >= jk; i--)
531 j |= iq[i];
532 if (j == 0)
533 {
534 for (k = 1; iq[jk - k] == 0; k++);
535
536 for (i = jz + 1; i <= jz + k; i++)
537 {
538 f[jx + i] = (double) ipio2[jv + i];
539 for (j = 0, fw = 0.0; j <= jx; j++)
540 fw += x[j] * f[jx + i - j];
541 q[i] = fw;
542 }
543 jz += k;
544 goto recompute;
545 }
546 }
547
548
549 if (z == 0.0)
550 {
551 jz -= 1;
552 q0 -= 24;
553 while (iq[jz] == 0)
554 {
555 jz--;
556 q0 -= 24;
557 }
558 }
559 else
560 {
561 z = ldexp (z, -q0);
562 if (z >= two24)
563 {
564 fw = (double) ((int) (twon24 * z));
565 iq[jz] = (int) (z - two24 * fw);
566 jz += 1;
567 q0 += 24;
568 iq[jz] = (int) fw;
569 }
570 else
571 iq[jz] = (int) z;
572 }
573
574
575 fw = ldexp (one, q0);
576 for (i = jz; i >= 0; i--)
577 {
578 q[i] = fw * (double) iq[i];
579 fw *= twon24;
580 }
581
582
583 for (i = jz; i >= 0; i--)
584 {
585 for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
586 fw += PIo2[k] * q[i + k];
587 fq[jz - i] = fw;
588 }
589
590
591 switch (prec)
592 {
593 case 0:
594 fw = 0.0;
595 for (i = jz; i >= 0; i--)
596 fw += fq[i];
597 y[0] = (ih == 0) ? fw : -fw;
598 break;
599 case 1:
600 case 2:
601 fw = 0.0;
602 for (i = jz; i >= 0; i--)
603 fw += fq[i];
604 y[0] = (ih == 0) ? fw : -fw;
605 fw = fq[0] - fw;
606 for (i = 1; i <= jz; i++)
607 fw += fq[i];
608 y[1] = (ih == 0) ? fw : -fw;
609 break;
610 case 3:
611 for (i = jz; i > 0; i--)
612 {
613 fw = fq[i - 1] + fq[i];
614 fq[i] += fq[i - 1] - fw;
615 fq[i - 1] = fw;
616 }
617 for (i = jz; i > 1; i--)
618 {
619 fw = fq[i - 1] + fq[i];
620 fq[i] += fq[i - 1] - fw;
621 fq[i - 1] = fw;
622 }
623 for (fw = 0.0, i = jz; i >= 2; i--)
624 fw += fq[i];
625 if (ih == 0)
626 {
627 y[0] = fq[0];
628 y[1] = fq[1];
629 y[2] = fw;
630 }
631 else
632 {
633 y[0] = -fq[0];
634 y[1] = -fq[1];
635 y[2] = -fw;
636 }
637 }
638 return n & 7;
639 }