206  *                      e0-3-24*jv >= 0 or (e0-3)/24 >= jv
 207  *              Hence jv = max(0,(e0-3)/24).
 208  *
 209  *      jp      jp+1 is the number of terms in PIo2[] needed, jp = jk.
 210  *
 211  *      q[]     double array with integral value, representing the
 212  *              24-bits chunk of the product of x and 2/pi.
 213  *
 214  *      q0      the corresponding exponent of q[0]. Note that the
 215  *              exponent for q[i] would be q0-24*i.
 216  *
 217  *      PIo2[]  double precision array, obtained by cutting pi/2
 218  *              into 24 bits chunks.
 219  *
 220  *      f[]     ipio2[] in floating point
 221  *
 222  *      iq[]    integer array by breaking up q[] in 24-bits chunk.
 223  *
 224  *      fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]
 225  *
 226  *      ih      integer. If >0 it indicats q[] is >= 0.5, hence
 227  *              it also indicates the *sign* of the result.
 228  *
 229  */
 230 
 231 
 232 /*
 233  * Constants:
 234  * The hexadecimal values are the intended ones for the following
 235  * constants. The decimal values may be used, provided that the
 236  * compiler will convert from decimal to binary accurately enough
 237  * to produce the hexadecimal values shown.
 238  */
 239 
 240 
 241 static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
 242 
 243 static const double PIo2[] = {
 244   1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
 245   7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
 246   5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */

 330   if(z==zeroB) {
 331     j = 0;
 332     for (i=jz-1;i>=jk;i--) j |= iq[i];
 333     if(j==0) { /* need recomputation */
 334       for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
 335 
 336       for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
 337         f[jx+i] = (double) ipio2[jv+i];
 338         for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
 339         q[i] = fw;
 340       }
 341       jz += k;
 342       goto recompute;
 343     }
 344   }
 345 
 346   /* chop off zero terms */
 347   if(z==0.0) {
 348     jz -= 1; q0 -= 24;
 349     while(iq[jz]==0) { jz--; q0-=24;}
 350   } else { /* break z into 24-bit if neccessary */
 351     z = scalbnA(z,-q0);
 352     if(z>=two24B) {
 353       fw = (double)((int)(twon24*z));
 354       iq[jz] = (int)(z-two24B*fw);
 355       jz += 1; q0 += 24;
 356       iq[jz] = (int) fw;
 357     } else iq[jz] = (int) z ;
 358   }
 359 
 360   /* convert integer "bit" chunk to floating-point value */
 361   fw = scalbnA(one,q0);
 362   for(i=jz;i>=0;i--) {
 363     q[i] = fw*(double)iq[i]; fw*=twon24;
 364   }
 365 
 366   /* compute PIo2[0,...,jp]*q[jz,...,0] */
 367   for(i=jz;i>=0;i--) {
 368     for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
 369     fq[jz-i] = fw;
 370   }

 392       fq[i-1] = fw;
 393     }
 394     for (i=jz;i>1;i--) {
 395       fw      = fq[i-1]+fq[i];
 396       fq[i]  += fq[i-1]-fw;
 397       fq[i-1] = fw;
 398     }
 399     for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
 400     if(ih==0) {
 401       y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
 402     } else {
 403       y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
 404     }
 405   }
 406   return n&7;
 407 }
 408 
 409 
 410 /*
 411  * ====================================================
 412  * Copyright (c) 1993 Oracle and/or its affilates. All rights reserved.
 413  *
 414  * Developed at SunPro, a Sun Microsystems, Inc. business.
 415  * Permission to use, copy, modify, and distribute this
 416  * software is freely granted, provided that this notice
 417  * is preserved.
 418  * ====================================================
 419  *
 420  */
 421 
 422 /* __ieee754_rem_pio2(x,y)
 423  *
 424  * return the remainder of x rem pi/2 in y[0]+y[1]
 425  * use __kernel_rem_pio2()
 426  */
 427 
 428 /*
 429  * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
 430  */
 431 static const int two_over_pi[] = {
 432   0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,

 206  *                      e0-3-24*jv >= 0 or (e0-3)/24 >= jv
 207  *              Hence jv = max(0,(e0-3)/24).
 208  *
 209  *      jp      jp+1 is the number of terms in PIo2[] needed, jp = jk.
 210  *
 211  *      q[]     double array with integral value, representing the
 212  *              24-bits chunk of the product of x and 2/pi.
 213  *
 214  *      q0      the corresponding exponent of q[0]. Note that the
 215  *              exponent for q[i] would be q0-24*i.
 216  *
 217  *      PIo2[]  double precision array, obtained by cutting pi/2
 218  *              into 24 bits chunks.
 219  *
 220  *      f[]     ipio2[] in floating point
 221  *
 222  *      iq[]    integer array by breaking up q[] in 24-bits chunk.
 223  *
 224  *      fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]
 225  *
 226  *      ih      integer. If >0 it indicates q[] is >= 0.5, hence
 227  *              it also indicates the *sign* of the result.
 228  *
 229  */
 230 
 231 
 232 /*
 233  * Constants:
 234  * The hexadecimal values are the intended ones for the following
 235  * constants. The decimal values may be used, provided that the
 236  * compiler will convert from decimal to binary accurately enough
 237  * to produce the hexadecimal values shown.
 238  */
 239 
 240 
 241 static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
 242 
 243 static const double PIo2[] = {
 244   1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
 245   7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
 246   5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */

 330   if(z==zeroB) {
 331     j = 0;
 332     for (i=jz-1;i>=jk;i--) j |= iq[i];
 333     if(j==0) { /* need recomputation */
 334       for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
 335 
 336       for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
 337         f[jx+i] = (double) ipio2[jv+i];
 338         for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
 339         q[i] = fw;
 340       }
 341       jz += k;
 342       goto recompute;
 343     }
 344   }
 345 
 346   /* chop off zero terms */
 347   if(z==0.0) {
 348     jz -= 1; q0 -= 24;
 349     while(iq[jz]==0) { jz--; q0-=24;}
 350   } else { /* break z into 24-bit if necessary */
 351     z = scalbnA(z,-q0);
 352     if(z>=two24B) {
 353       fw = (double)((int)(twon24*z));
 354       iq[jz] = (int)(z-two24B*fw);
 355       jz += 1; q0 += 24;
 356       iq[jz] = (int) fw;
 357     } else iq[jz] = (int) z ;
 358   }
 359 
 360   /* convert integer "bit" chunk to floating-point value */
 361   fw = scalbnA(one,q0);
 362   for(i=jz;i>=0;i--) {
 363     q[i] = fw*(double)iq[i]; fw*=twon24;
 364   }
 365 
 366   /* compute PIo2[0,...,jp]*q[jz,...,0] */
 367   for(i=jz;i>=0;i--) {
 368     for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
 369     fq[jz-i] = fw;
 370   }

 392       fq[i-1] = fw;
 393     }
 394     for (i=jz;i>1;i--) {
 395       fw      = fq[i-1]+fq[i];
 396       fq[i]  += fq[i-1]-fw;
 397       fq[i-1] = fw;
 398     }
 399     for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
 400     if(ih==0) {
 401       y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
 402     } else {
 403       y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
 404     }
 405   }
 406   return n&7;
 407 }
 408 
 409 
 410 /*
 411  * ====================================================
 412  * Copyright (c) 1993 Oracle and/or its affiliates. All rights reserved.
 413  *
 414  * Developed at SunPro, a Sun Microsystems, Inc. business.
 415  * Permission to use, copy, modify, and distribute this
 416  * software is freely granted, provided that this notice
 417  * is preserved.
 418  * ====================================================
 419  *
 420  */
 421 
 422 /* __ieee754_rem_pio2(x,y)
 423  *
 424  * return the remainder of x rem pi/2 in y[0]+y[1]
 425  * use __kernel_rem_pio2()
 426  */
 427 
 428 /*
 429  * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
 430  */
 431 static const int two_over_pi[] = {
 432   0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,