libquadmath/math/rem_pio2q.c

   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  * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
  19  * double x[],y[]; int e0,nx,prec; int ipio2[];
  20  *
  21  * __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                     fw      = fq[i-1]+fq[i];
 286                     fq[i]  += fq[i-1]-fw;
 287                     fq[i-1] = fw;
 288                 }
 289                 for (i=jz;i>1;i--) {
 290                     fw      = fq[i-1]+fq[i];
 291                     fq[i]  += fq[i-1]-fw;
 292                     fq[i-1] = fw;
 293                 }
 294                 for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
 295                 if(ih==0) {
 296                     y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
 297                 } else {
 298                     y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
 299                 }
 300         }
 301         return n&7;
 302 }
 303
 304
 305 \f
 306
 307
 308 /* Quad-precision floating point argument reduction.
 309    Copyright (C) 1999 Free Software Foundation, Inc.
 310    This file is part of the GNU C Library.
 311    Contributed by Jakub Jelinek <jj@ultra.linux.cz>
 312
 313    The GNU C Library is free software; you can redistribute it and/or
 314    modify it under the terms of the GNU Lesser General Public
 315    License as published by the Free Software Foundation; either
 316    version 2.1 of the License, or (at your option) any later version.
 317
 318    The GNU C Library is distributed in the hope that it will be useful,
 319    but WITHOUT ANY WARRANTY; without even the implied warranty of
 320    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 321    Lesser General Public License for more details.
 322
 323    You should have received a copy of the GNU Lesser General Public
 324    License along with the GNU C Library; if not, write to the Free
 325    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 326    02111-1307 USA.  */
 327
 328 /*
 329  * Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi
 330  */
 331 static const int32_t two_over_pi[] = {
 332 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
 333 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
 334 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
 335 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
 336 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
 337 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
 338 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
 339 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
 340 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
 341 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
 342 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
 343 0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6,
 344 0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2,
 345 0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35,
 346 0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30,
 347 0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c,
 348 0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4,
 349 0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770,
 350 0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7,
 351 0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19,
 352 0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522,
 353 0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16,
 354 0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6,
 355 0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e,
 356 0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48,
 357 0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3,
 358 0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf,
 359 0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55,
 360 0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612,
 361 0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929,
 362 0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec,
 363 0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b,
 364 0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c,
 365 0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4,
 366 0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb,
 367 0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc,
 368 0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c,
 369 0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f,
 370 0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5,
 371 0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437,
 372 0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b,
 373 0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea,
 374 0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad,
 375 0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3,
 376 0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3,
 377 0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717,
 378 0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f,
 379 0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61,
 380 0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db,
 381 0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51,
 382 0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0,
 383 0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c,
 384 0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6,
 385 0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc,
 386 0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed,
 387 0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328,
 388 0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d,
 389 0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0,
 390 0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b,
 391 0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4,
 392 0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3,
 393 0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f,
 394 0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad,
 395 0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b,
 396 0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4,
 397 0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761,
 398 0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31,
 399 0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30,
 400 0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262,
 401 0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e,
 402 0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1,
 403 0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c,
 404 0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4,
 405 0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08,
 406 0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196,
 407 0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9,
 408 0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4,
 409 0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc,
 410 0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c,
 411 0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0,
 412 0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c,
 413 0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0,
 414 0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac,
 415 0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22,
 416 0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893,
 417 0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7,
 418 0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5,
 419 0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f,
 420 0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4,
 421 0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf,
 422 0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b,
 423 0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2,
 424 0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138,
 425 0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e,
 426 0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569,
 427 0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34,
 428 0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9,
 429 0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d,
 430 0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f,
 431 0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855,
 432 0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569,
 433 0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b,
 434 0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe,
 435 0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41,
 436 0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49,
 437 0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f,
 438 0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110,
 439 0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8,
 440 0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365,
 441 0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a,
 442 0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270,
 443 0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5,
 444 0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616,
 445 0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b,
 446 0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0,
 447 0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb,
 448 0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a,
 449 0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e,
 450 0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa,
 451 0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5,
 452 0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0,
 453 0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2,
 454 0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886,
 455 0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142,
 456 0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba,
 457 0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4,
 458 0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708,
 459 0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555,
 460 0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3,
 461 0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55,
 462 0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58,
 463 0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5,
 464 0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c,
 465 0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe,
 466 0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b,
 467 0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8,
 468 0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005,
 469 0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7,
 470 0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50,
 471 0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604,
 472 0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643,
 473 0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485,
 474 0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d,
 475 0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6,
 476 0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2,
 477 0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02,
 478 0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3,
 479 0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412,
 480 0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274,
 481 0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755,
 482 0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849,
 483 0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce,
 484 0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5,
 485 0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba,
 486 0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6,
 487 0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d,
 488 0x7b7b89, 0x483d38,
 489 };
 490
 491 static const __float128 c[] = {
 492 /* 93 bits of pi/2 */
 493 #define PI_2_1 c[0]
 494  1.57079632679489661923132169155131424e+00Q, /* 3fff921fb54442d18469898cc5100000 */
 495
 496 /* pi/2 - PI_2_1 */
 497 #define PI_2_1t c[1]
 498  8.84372056613570112025531863263659260e-29Q, /* 3fa1c06e0e68948127044533e63a0106 */
 499 };
 500
 501
 502 int32_t
 503 __quadmath_rem_pio2q (__float128 x, __float128 *y)
 504 {
 505   __float128 z, w, t;
 506   double tx[8];
 507   int64_t exp, n, ix, hx;
 508   uint64_t lx;
 509
 510   GET_FLT128_WORDS64 (hx, lx, x);
 511   ix = hx & 0x7fffffffffffffffLL;
 512   if (ix <= 0x3ffe921fb54442d1LL)       /* x in <-pi/4, pi/4> */
 513     {
 514       y[0] = x;
 515       y[1] = 0;
 516       return 0;
 517     }
 518
 519   if (ix < 0x40002d97c7f3321dLL)        /* |x| in <pi/4, 3pi/4) */
 520     {
 521       if (hx > 0)
 522         {
 523           /* 113 + 93 bit PI is ok */
 524           z = x - PI_2_1;
 525           y[0] = z - PI_2_1t;
 526           y[1] = (z - y[0]) - PI_2_1t;
 527           return 1;
 528         }
 529       else
 530         {
 531           /* 113 + 93 bit PI is ok */
 532           z = x + PI_2_1;
 533           y[0] = z + PI_2_1t;
 534           y[1] = (z - y[0]) + PI_2_1t;
 535           return -1;
 536         }
 537     }
 538
 539   if (ix >= 0x7fff000000000000LL)       /* x is +=oo or NaN */
 540     {
 541       y[0] = x - x;
 542       y[1] = y[0];
 543       return 0;
 544     }
 545
 546   /* Handle large arguments.
 547      We split the 113 bits of the mantissa into 5 24bit integers
 548      stored in a double array.  */
 549   exp = (ix >> 48) - 16383 - 23;
 550
 551   /* This is faster than doing this in floating point, because we
 552      have to convert it to integers anyway and like this we can keep
 553      both integer and floating point units busy.  */
 554   tx [0] = (double)(((ix >> 25) & 0x7fffff) | 0x800000);
 555   tx [1] = (double)((ix >> 1) & 0xffffff);
 556   tx [2] = (double)(((ix << 23) | (lx >> 41)) & 0xffffff);
 557   tx [3] = (double)((lx >> 17) & 0xffffff);
 558   tx [4] = (double)((lx << 7) & 0xffffff);
 559
 560   n = __quadmath_kernel_rem_pio2 (tx, tx + 5, exp,
 561                                   ((lx << 7) & 0xffffff) ? 5 : 4,
 562                                   3, two_over_pi);
 563
 564   /* The result is now stored in 3 double values, we need to convert it into
 565      two __float128 values.  */
 566   t = (__float128) tx [6] + (__float128) tx [7];
 567   w = (__float128) tx [5];
 568
 569   if (hx >= 0)
 570     {
 571       y[0] = w + t;
 572       y[1] = t - (y[0] - w);
 573       return n;
 574     }
 575   else
 576     {
 577       y[0] = -(w + t);
 578       y[1] = -t - (y[0] + w);
 579       return -n;
 580     }
 581 }