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 133static const int init_jk[] = {2,3,4,6}; /* initial value for jk */ 134 135static 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 146static 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 153static 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; 178recompute: 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 */ 337static const int32_t two_over_pi[] = { 3380xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62, 3390x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a, 3400x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129, 3410xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41, 3420x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8, 3430x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf, 3440x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5, 3450xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08, 3460x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3, 3470x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880, 3480x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b, 3490x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6, 3500xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2, 3510xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35, 3520xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30, 3530x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c, 3540x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4, 3550xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770, 3560xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7, 3570xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19, 3580xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522, 3590x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16, 3600xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6, 3610x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e, 3620x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48, 3630xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3, 3640xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf, 3650xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55, 3660x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612, 3670xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929, 3680x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec, 3690xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b, 3700x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c, 3710x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4, 3720x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb, 3730xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc, 3740x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c, 3750x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f, 3760xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5, 3770x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437, 3780x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b, 3790x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea, 3800x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad, 3810x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3, 3820x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3, 3830xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717, 3840x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f, 3850xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61, 3860xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db, 3870xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51, 3880xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0, 3890x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c, 3900x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6, 3910x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc, 3920xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed, 3930x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328, 3940x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d, 3950x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0, 3960xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b, 3970x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4, 3980xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3, 3990xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f, 4000x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad, 4010x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b, 4020x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4, 4030xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761, 4040x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31, 4050x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30, 4060xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262, 4070x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e, 4080xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1, 4090xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c, 4100x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4, 4110xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08, 4120x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196, 4130xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9, 4140x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4, 4150x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc, 4160x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c, 4170xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0, 4180x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c, 4190xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0, 4200x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac, 4210xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22, 4220x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893, 4230x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7, 4240xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5, 4250x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f, 4260xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4, 4270x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf, 4280x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b, 4290x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2, 4300x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138, 4310xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e, 4320xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569, 4330x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34, 4340xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9, 4350x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d, 4360xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f, 4370xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855, 4380x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569, 4390xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b, 4400x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe, 4410x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41, 4420x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49, 4430xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f, 4440xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110, 4450x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8, 4460x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365, 4470xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a, 4480x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270, 4490x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5, 4500x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616, 4510x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b, 4520x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0, 4530xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb, 4540x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a, 4550x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e, 4560xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa, 4570x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5, 4580xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0, 4590xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2, 4600xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886, 4610x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142, 4620xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba, 4630xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4, 4640xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708, 4650x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555, 4660xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3, 4670x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55, 4680xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58, 4690xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5, 4700x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c, 4710x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe, 4720x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b, 4730x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8, 4740x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005, 4750x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7, 4760xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50, 4770x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604, 4780x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643, 4790x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485, 4800xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d, 4810xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6, 4820xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2, 4830x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02, 4840xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3, 4850x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412, 4860xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274, 4870x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755, 4880xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849, 4890xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce, 4900xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5, 4910x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba, 4920xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6, 4930xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d, 4940x7b7b89, 0x483d38, 495}; 496 497static 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 508int32_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