1 | #include "quadmath-imp.h" |
---|---|

2 | #include <math.h> |

3 | |

4 | |

5 | /* @(#)k_rem_pio2.c 5.1 93/09/24 */ |

6 | /* |

7 | * ==================================================== |

8 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |

9 | * |

10 | * Developed at SunPro, a Sun Microsystems, Inc. business. |

11 | * Permission to use, copy, modify, and distribute this |

12 | * software is freely granted, provided that this notice |

13 | * is preserved. |

14 | * ==================================================== |

15 | */ |

16 | |

17 | /* |

18 | * __quadmath_kernel_rem_pio2 (x,y,e0,nx,prec,ipio2) |

19 | * double x[],y[]; int e0,nx,prec; int ipio2[]; |

20 | * |

21 | * __quadmath_kernel_rem_pio2 return the last three digits of N with |

22 | * y = x - N*pi/2 |

23 | * so that |y| < pi/2. |

24 | * |

25 | * The method is to compute the integer (mod 8) and fraction parts of |

26 | * (2/pi)*x without doing the full multiplication. In general we |

27 | * skip the part of the product that are known to be a huge integer ( |

28 | * more accurately, = 0 mod 8 ). Thus the number of operations are |

29 | * independent of the exponent of the input. |

30 | * |

31 | * (2/pi) is represented by an array of 24-bit integers in ipio2[]. |

32 | * |

33 | * Input parameters: |

34 | * x[] The input value (must be positive) is broken into nx |

35 | * pieces of 24-bit integers in double precision format. |

36 | * x[i] will be the i-th 24 bit of x. The scaled exponent |

37 | * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 |

38 | * match x's up to 24 bits. |

39 | * |

40 | * Example of breaking a double positive z into x[0]+x[1]+x[2]: |

41 | * e0 = ilogb(z)-23 |

42 | * z = scalbn(z,-e0) |

43 | * for i = 0,1,2 |

44 | * x[i] = floor(z) |

45 | * z = (z-x[i])*2**24 |

46 | * |

47 | * |

48 | * y[] ouput result in an array of double precision numbers. |

49 | * The dimension of y[] is: |

50 | * 24-bit precision 1 |

51 | * 53-bit precision 2 |

52 | * 64-bit precision 2 |

53 | * 113-bit precision 3 |

54 | * The actual value is the sum of them. Thus for 113-bit |

55 | * precision, one may have to do something like: |

56 | * |

57 | * long double t,w,r_head, r_tail; |

58 | * t = (long double)y[2] + (long double)y[1]; |

59 | * w = (long double)y[0]; |

60 | * r_head = t+w; |

61 | * r_tail = w - (r_head - t); |

62 | * |

63 | * e0 The exponent of x[0] |

64 | * |

65 | * nx dimension of x[] |

66 | * |

67 | * prec an integer indicating the precision: |

68 | * 0 24 bits (single) |

69 | * 1 53 bits (double) |

70 | * 2 64 bits (extended) |

71 | * 3 113 bits (quad) |

72 | * |

73 | * ipio2[] |

74 | * integer array, contains the (24*i)-th to (24*i+23)-th |

75 | * bit of 2/pi after binary point. The corresponding |

76 | * floating value is |

77 | * |

78 | * ipio2[i] * 2^(-24(i+1)). |

79 | * |

80 | * External function: |

81 | * double scalbn(), floor(); |

82 | * |

83 | * |

84 | * Here is the description of some local variables: |

85 | * |

86 | * jk jk+1 is the initial number of terms of ipio2[] needed |

87 | * in the computation. The recommended value is 2,3,4, |

88 | * 6 for single, double, extended,and quad. |

89 | * |

90 | * jz local integer variable indicating the number of |

91 | * terms of ipio2[] used. |

92 | * |

93 | * jx nx - 1 |

94 | * |

95 | * jv index for pointing to the suitable ipio2[] for the |

96 | * computation. In general, we want |

97 | * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 |

98 | * is an integer. Thus |

99 | * e0-3-24*jv >= 0 or (e0-3)/24 >= jv |

100 | * Hence jv = max(0,(e0-3)/24). |

101 | * |

102 | * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. |

103 | * |

104 | * q[] double array with integral value, representing the |

105 | * 24-bits chunk of the product of x and 2/pi. |

106 | * |

107 | * q0 the corresponding exponent of q[0]. Note that the |

108 | * exponent for q[i] would be q0-24*i. |

109 | * |

110 | * PIo2[] double precision array, obtained by cutting pi/2 |

111 | * into 24 bits chunks. |

112 | * |

113 | * f[] ipio2[] in floating point |

114 | * |

115 | * iq[] integer array by breaking up q[] in 24-bits chunk. |

116 | * |

117 | * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] |

118 | * |

119 | * ih integer. If >0 it indicates q[] is >= 0.5, hence |

120 | * it also indicates the *sign* of the result. |

121 | * |

122 | */ |

123 | |

124 | /* |

125 | * Constants: |

126 | * The hexadecimal values are the intended ones for the following |

127 | * constants. The decimal values may be used, provided that the |

128 | * compiler will convert from decimal to binary accurately enough |

129 | * to produce the hexadecimal values shown. |

130 | */ |

131 | |

132 | |

133 | static const int init_jk[] = {2,3,4,6}; /* initial value for jk */ |

134 | |

135 | static const double PIo2[] = { |

136 | 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ |

137 | 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ |

138 | 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ |

139 | 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ |

140 | 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ |

141 | 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */ |

142 | 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */ |

143 | 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ |

144 | }; |

145 | |

146 | static const double |

147 | zero = 0.0, |

148 | one = 1.0, |

149 | two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ |

150 | twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ |

151 | |

152 | |

153 | static int |

154 | __quadmath_kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2) |

155 | { |

156 | int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; |

157 | double z,fw,f[20],fq[20],q[20]; |

158 | |

159 | /* initialize jk*/ |

160 | jk = init_jk[prec]; |

161 | jp = jk; |

162 | |

163 | /* determine jx,jv,q0, note that 3>q0 */ |

164 | jx = nx-1; |

165 | jv = (e0-3)/24; if(jv<0) jv=0; |

166 | q0 = e0-24*(jv+1); |

167 | |

168 | /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ |

169 | j = jv-jx; m = jx+jk; |

170 | for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j]; |

171 | |

172 | /* compute q[0],q[1],...q[jk] */ |

173 | for (i=0;i<=jk;i++) { |

174 | for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw; |

175 | } |

176 | |

177 | jz = jk; |

178 | recompute: |

179 | /* distill q[] into iq[] reversingly */ |

180 | for(i=0,j=jz,z=q[jz];j>0;i++,j--) { |

181 | fw = (double)((int32_t)(twon24* z)); |

182 | iq[i] = (int32_t)(z-two24*fw); |

183 | z = q[j-1]+fw; |

184 | } |

185 | |

186 | /* compute n */ |

187 | z = scalbn(z,q0); /* actual value of z */ |

188 | z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */ |

189 | n = (int32_t) z; |

190 | z -= (double)n; |

191 | ih = 0; |

192 | if(q0>0) { /* need iq[jz-1] to determine n */ |

193 | i = (iq[jz-1]>>(24-q0)); n += i; |

194 | iq[jz-1] -= i<<(24-q0); |

195 | ih = iq[jz-1]>>(23-q0); |

196 | } |

197 | else if(q0==0) ih = iq[jz-1]>>23; |

198 | else if(z>=0.5) ih=2; |

199 | |

200 | if(ih>0) { /* q > 0.5 */ |

201 | n += 1; carry = 0; |

202 | for(i=0;i<jz ;i++) { /* compute 1-q */ |

203 | j = iq[i]; |

204 | if(carry==0) { |

205 | if(j!=0) { |

206 | carry = 1; iq[i] = 0x1000000- j; |

207 | } |

208 | } else iq[i] = 0xffffff - j; |

209 | } |

210 | if(q0>0) { /* rare case: chance is 1 in 12 */ |

211 | switch(q0) { |

212 | case 1: |

213 | iq[jz-1] &= 0x7fffff; break; |

214 | case 2: |

215 | iq[jz-1] &= 0x3fffff; break; |

216 | } |

217 | } |

218 | if(ih==2) { |

219 | z = one - z; |

220 | if(carry!=0) z -= scalbn(one,q0); |

221 | } |

222 | } |

223 | |

224 | /* check if recomputation is needed */ |

225 | if(z==zero) { |

226 | j = 0; |

227 | for (i=jz-1;i>=jk;i--) j |= iq[i]; |

228 | if(j==0) { /* need recomputation */ |

229 | for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */ |

230 | |

231 | for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */ |

232 | f[jx+i] = (double) ipio2[jv+i]; |

233 | for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; |

234 | q[i] = fw; |

235 | } |

236 | jz += k; |

237 | goto recompute; |

238 | } |

239 | } |

240 | |

241 | /* chop off zero terms */ |

242 | if(z==0.0) { |

243 | jz -= 1; q0 -= 24; |

244 | while(iq[jz]==0) { jz--; q0-=24;} |

245 | } else { /* break z into 24-bit if necessary */ |

246 | z = scalbn(z,-q0); |

247 | if(z>=two24) { |

248 | fw = (double)((int32_t)(twon24*z)); |

249 | iq[jz] = (int32_t)(z-two24*fw); |

250 | jz += 1; q0 += 24; |

251 | iq[jz] = (int32_t) fw; |

252 | } else iq[jz] = (int32_t) z ; |

253 | } |

254 | |

255 | /* convert integer "bit" chunk to floating-point value */ |

256 | fw = scalbn(one,q0); |

257 | for(i=jz;i>=0;i--) { |

258 | q[i] = fw*(double)iq[i]; fw*=twon24; |

259 | } |

260 | |

261 | /* compute PIo2[0,...,jp]*q[jz,...,0] */ |

262 | for(i=jz;i>=0;i--) { |

263 | for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k]; |

264 | fq[jz-i] = fw; |

265 | } |

266 | |

267 | /* compress fq[] into y[] */ |

268 | switch(prec) { |

269 | case 0: |

270 | fw = 0.0; |

271 | for (i=jz;i>=0;i--) fw += fq[i]; |

272 | y[0] = (ih==0)? fw: -fw; |

273 | break; |

274 | case 1: |

275 | case 2: |

276 | fw = 0.0; |

277 | for (i=jz;i>=0;i--) fw += fq[i]; |

278 | y[0] = (ih==0)? fw: -fw; |

279 | fw = fq[0]-fw; |

280 | for (i=1;i<=jz;i++) fw += fq[i]; |

281 | y[1] = (ih==0)? fw: -fw; |

282 | break; |

283 | case 3: /* painful */ |

284 | for (i=jz;i>0;i--) { |

285 | #if __FLT_EVAL_METHOD__ != 0 |

286 | volatile |

287 | #endif |

288 | double fv = (double)(fq[i-1]+fq[i]); |

289 | fq[i] += fq[i-1]-fv; |

290 | fq[i-1] = fv; |

291 | } |

292 | for (i=jz;i>1;i--) { |

293 | #if __FLT_EVAL_METHOD__ != 0 |

294 | volatile |

295 | #endif |

296 | double fv = (double)(fq[i-1]+fq[i]); |

297 | fq[i] += fq[i-1]-fv; |

298 | fq[i-1] = fv; |

299 | } |

300 | for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; |

301 | if(ih==0) { |

302 | y[0] = fq[0]; y[1] = fq[1]; y[2] = fw; |

303 | } else { |

304 | y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw; |

305 | } |

306 | } |

307 | return n&7; |

308 | } |

309 | |

310 | |

311 | |

312 | |

313 | |

314 | /* Quad-precision floating point argument reduction. |

315 | Copyright (C) 1999-2017 Free Software Foundation, Inc. |

316 | This file is part of the GNU C Library. |

317 | Contributed by Jakub Jelinek <jj@ultra.linux.cz> |

318 | |

319 | The GNU C Library is free software; you can redistribute it and/or |

320 | modify it under the terms of the GNU Lesser General Public |

321 | License as published by the Free Software Foundation; either |

322 | version 2.1 of the License, or (at your option) any later version. |

323 | |

324 | The GNU C Library is distributed in the hope that it will be useful, |

325 | but WITHOUT ANY WARRANTY; without even the implied warranty of |

326 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

327 | Lesser General Public License for more details. |

328 | |

329 | You should have received a copy of the GNU Lesser General Public |

330 | License along with the GNU C Library; if not, write to the Free |

331 | Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA |

332 | 02111-1307 USA. */ |

333 | |

334 | /* |

335 | * Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi |

336 | */ |

337 | static const int32_t two_over_pi[] = { |

338 | 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62, |

339 | 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a, |

340 | 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129, |

341 | 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41, |

342 | 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8, |

343 | 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf, |

344 | 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5, |

345 | 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08, |

346 | 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3, |

347 | 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880, |

348 | 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b, |

349 | 0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6, |

350 | 0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2, |

351 | 0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35, |

352 | 0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30, |

353 | 0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c, |

354 | 0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4, |

355 | 0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770, |

356 | 0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7, |

357 | 0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19, |

358 | 0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522, |

359 | 0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16, |

360 | 0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6, |

361 | 0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e, |

362 | 0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48, |

363 | 0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3, |

364 | 0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf, |

365 | 0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55, |

366 | 0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612, |

367 | 0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929, |

368 | 0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec, |

369 | 0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b, |

370 | 0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c, |

371 | 0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4, |

372 | 0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb, |

373 | 0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc, |

374 | 0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c, |

375 | 0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f, |

376 | 0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5, |

377 | 0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437, |

378 | 0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b, |

379 | 0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea, |

380 | 0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad, |

381 | 0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3, |

382 | 0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3, |

383 | 0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717, |

384 | 0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f, |

385 | 0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61, |

386 | 0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db, |

387 | 0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51, |

388 | 0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0, |

389 | 0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c, |

390 | 0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6, |

391 | 0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc, |

392 | 0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed, |

393 | 0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328, |

394 | 0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d, |

395 | 0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0, |

396 | 0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b, |

397 | 0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4, |

398 | 0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3, |

399 | 0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f, |

400 | 0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad, |

401 | 0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b, |

402 | 0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4, |

403 | 0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761, |

404 | 0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31, |

405 | 0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30, |

406 | 0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262, |

407 | 0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e, |

408 | 0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1, |

409 | 0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c, |

410 | 0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4, |

411 | 0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08, |

412 | 0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196, |

413 | 0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9, |

414 | 0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4, |

415 | 0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc, |

416 | 0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c, |

417 | 0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0, |

418 | 0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c, |

419 | 0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0, |

420 | 0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac, |

421 | 0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22, |

422 | 0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893, |

423 | 0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7, |

424 | 0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5, |

425 | 0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f, |

426 | 0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4, |

427 | 0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf, |

428 | 0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b, |

429 | 0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2, |

430 | 0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138, |

431 | 0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e, |

432 | 0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569, |

433 | 0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34, |

434 | 0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9, |

435 | 0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d, |

436 | 0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f, |

437 | 0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855, |

438 | 0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569, |

439 | 0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b, |

440 | 0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe, |

441 | 0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41, |

442 | 0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49, |

443 | 0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f, |

444 | 0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110, |

445 | 0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8, |

446 | 0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365, |

447 | 0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a, |

448 | 0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270, |

449 | 0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5, |

450 | 0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616, |

451 | 0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b, |

452 | 0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0, |

453 | 0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb, |

454 | 0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a, |

455 | 0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e, |

456 | 0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa, |

457 | 0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5, |

458 | 0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0, |

459 | 0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2, |

460 | 0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886, |

461 | 0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142, |

462 | 0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba, |

463 | 0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4, |

464 | 0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708, |

465 | 0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555, |

466 | 0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3, |

467 | 0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55, |

468 | 0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58, |

469 | 0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5, |

470 | 0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c, |

471 | 0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe, |

472 | 0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b, |

473 | 0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8, |

474 | 0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005, |

475 | 0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7, |

476 | 0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50, |

477 | 0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604, |

478 | 0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643, |

479 | 0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485, |

480 | 0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d, |

481 | 0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6, |

482 | 0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2, |

483 | 0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02, |

484 | 0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3, |

485 | 0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412, |

486 | 0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274, |

487 | 0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755, |

488 | 0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849, |

489 | 0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce, |

490 | 0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5, |

491 | 0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba, |

492 | 0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6, |

493 | 0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d, |

494 | 0x7b7b89, 0x483d38, |

495 | }; |

496 | |

497 | static const __float128 c[] = { |

498 | /* 113 bits of pi/2 */ |

499 | #define PI_2_1 c[0] |

500 | 0x1.921fb54442d18469898cc51701b8p+0Q, |

501 | |

502 | /* pi/2 - PI_2_1 */ |

503 | #define PI_2_1t c[1] |

504 | 0x3.9a252049c1114cf98e804177d4c8p-116Q, |

505 | }; |

506 | |

507 | |

508 | int32_t |

509 | __quadmath_rem_pio2q (__float128 x, __float128 *y) |

510 | { |

511 | __float128 z, w, t; |

512 | double tx[8]; |

513 | int64_t exp, n, ix, hx; |

514 | uint64_t lx; |

515 | |

516 | GET_FLT128_WORDS64 (hx, lx, x); |

517 | ix = hx & 0x7fffffffffffffffLL; |

518 | if (ix <= 0x3ffe921fb54442d1LL) /* x in <-pi/4, pi/4> */ |

519 | { |

520 | y[0] = x; |

521 | y[1] = 0; |

522 | return 0; |

523 | } |

524 | |

525 | if (ix < 0x40002d97c7f3321dLL) /* |x| in <pi/4, 3pi/4) */ |

526 | { |

527 | if (hx > 0) |

528 | { |

529 | /* 113 + 113 bit PI is ok */ |

530 | z = x - PI_2_1; |

531 | y[0] = z - PI_2_1t; |

532 | y[1] = (z - y[0]) - PI_2_1t; |

533 | return 1; |

534 | } |

535 | else |

536 | { |

537 | /* 113 + 113 bit PI is ok */ |

538 | z = x + PI_2_1; |

539 | y[0] = z + PI_2_1t; |

540 | y[1] = (z - y[0]) + PI_2_1t; |

541 | return -1; |

542 | } |

543 | } |

544 | |

545 | if (ix >= 0x7fff000000000000LL) /* x is +=oo or NaN */ |

546 | { |

547 | y[0] = x - x; |

548 | y[1] = y[0]; |

549 | return 0; |

550 | } |

551 | |

552 | /* Handle large arguments. |

553 | We split the 113 bits of the mantissa into 5 24bit integers |

554 | stored in a double array. */ |

555 | exp = (ix >> 48) - 16383 - 23; |

556 | |

557 | /* This is faster than doing this in floating point, because we |

558 | have to convert it to integers anyway and like this we can keep |

559 | both integer and floating point units busy. */ |

560 | tx [0] = (double)(((ix >> 25) & 0x7fffff) | 0x800000); |

561 | tx [1] = (double)((ix >> 1) & 0xffffff); |

562 | tx [2] = (double)(((ix << 23) | (lx >> 41)) & 0xffffff); |

563 | tx [3] = (double)((lx >> 17) & 0xffffff); |

564 | tx [4] = (double)((lx << 7) & 0xffffff); |

565 | |

566 | n = __quadmath_kernel_rem_pio2 (tx, tx + 5, exp, |

567 | ((lx << 7) & 0xffffff) ? 5 : 4, |

568 | 3, two_over_pi); |

569 | |

570 | /* The result is now stored in 3 double values, we need to convert it into |

571 | two __float128 values. */ |

572 | t = (__float128) tx [6] + (__float128) tx [7]; |

573 | w = (__float128) tx [5]; |

574 | |

575 | if (hx >= 0) |

576 | { |

577 | y[0] = w + t; |

578 | y[1] = t - (y[0] - w); |

579 | return n; |

580 | } |

581 | else |

582 | { |

583 | y[0] = -(w + t); |

584 | y[1] = -t - (y[0] + w); |

585 | return -n; |

586 | } |

587 | } |

588 |