src/dom/js/fdlibm/e_rem_pio2.c

   1 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
   2  *
   3  * ***** BEGIN LICENSE BLOCK *****
   4  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
   5  *
   6  * The contents of this file are subject to the Mozilla Public License Version
   7  * 1.1 (the "License"); you may not use this file except in compliance with
   8  * the License. You may obtain a copy of the License at
   9  * http://www.mozilla.org/MPL/
  10  *
  11  * Software distributed under the License is distributed on an "AS IS" basis,
  12  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
  13  * for the specific language governing rights and limitations under the
  14  * License.
  15  *
  16  * The Original Code is Mozilla Communicator client code, released
  17  * March 31, 1998.
  18  *
  19  * The Initial Developer of the Original Code is
  20  * Sun Microsystems, Inc.
  21  * Portions created by the Initial Developer are Copyright (C) 1998
  22  * the Initial Developer. All Rights Reserved.
  23  *
  24  * Contributor(s):
  25  *
  26  * Alternatively, the contents of this file may be used under the terms of
  27  * either of the GNU General Public License Version 2 or later (the "GPL"),
  28  * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
  29  * in which case the provisions of the GPL or the LGPL are applicable instead
  30  * of those above. If you wish to allow use of your version of this file only
  31  * under the terms of either the GPL or the LGPL, and not to allow others to
  32  * use your version of this file under the terms of the MPL, indicate your
  33  * decision by deleting the provisions above and replace them with the notice
  34  * and other provisions required by the GPL or the LGPL. If you do not delete
  35  * the provisions above, a recipient may use your version of this file under
  36  * the terms of any one of the MPL, the GPL or the LGPL.
  37  *
  38  * ***** END LICENSE BLOCK ***** */
  39
  40 /* @(#)e_rem_pio2.c 1.4 95/01/18 */
  41 /*
  42  * ====================================================
  43  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  44  *
  45  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  46  * Permission to use, copy, modify, and distribute this
  47  * software is freely granted, provided that this notice
  48  * is preserved.
  49  * ====================================================
  50  *
  51  */
  52
  53 /* __ieee754_rem_pio2(x,y)
  54  *
  55  * return the remainder of x rem pi/2 in y[0]+y[1]
  56  * use __kernel_rem_pio2()
  57  */
  58
  59 #include "fdlibm.h"
  60
  61 /*
  62  * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
  63  */
  64 #ifdef __STDC__
  65 static const int two_over_pi[] = {
  66 #else
  67 static int two_over_pi[] = {
  68 #endif
  69 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
  70 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
  71 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
  72 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
  73 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
  74 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
  75 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
  76 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
  77 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
  78 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
  79 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
  80 };
  81
  82 #ifdef __STDC__
  83 static const int npio2_hw[] = {
  84 #else
  85 static int npio2_hw[] = {
  86 #endif
  87 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
  88 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
  89 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
  90 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
  91 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
  92 0x404858EB, 0x404921FB,
  93 };
  94
  95 /*
  96  * invpio2:  53 bits of 2/pi
  97  * pio2_1:   first  33 bit of pi/2
  98  * pio2_1t:  pi/2 - pio2_1
  99  * pio2_2:   second 33 bit of pi/2
 100  * pio2_2t:  pi/2 - (pio2_1+pio2_2)
 101  * pio2_3:   third  33 bit of pi/2
 102  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
 103  */
 104
 105 #ifdef __STDC__
 106 static const double
 107 #else
 108 static double
 109 #endif
 110 zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
 111 half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 112 two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
 113 invpio2 =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
 114 pio2_1  =  1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
 115 pio2_1t =  6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
 116 pio2_2  =  6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
 117 pio2_2t =  2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
 118 pio2_3  =  2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
 119 pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 120
 121 #ifdef __STDC__
 122         int __ieee754_rem_pio2(double x, double *y)
 123 #else
 124         int __ieee754_rem_pio2(x,y)
 125         double x,y[];
 126 #endif
 127 {
 128         fd_twoints u, ux, uz;
 129         double z,w,t,r,fn;
 130         double tx[3];
 131         int e0,i,j,nx,n,ix,hx;
 132
 133         u.d = x;
 134         hx = __HI(u);           /* high word of x */
 135         ix = hx&0x7fffffff;
 136         if(ix<=0x3fe921fb)   /* |x| ~<= pi/4 , no need for reduction */
 137             {y[0] = x; y[1] = 0; return 0;}
 138         if(ix<0x4002d97c) {  /* |x| < 3pi/4, special case with n=+-1 */
 139             if(hx>0) {
 140                 z = x - pio2_1;
 141                 if(ix!=0x3ff921fb) {    /* 33+53 bit pi is good enough */
 142                     y[0] = z - pio2_1t;
 143                     y[1] = (z-y[0])-pio2_1t;
 144                 } else {                /* near pi/2, use 33+33+53 bit pi */
 145                     z -= pio2_2;
 146                     y[0] = z - pio2_2t;
 147                     y[1] = (z-y[0])-pio2_2t;
 148                 }
 149                 return 1;
 150             } else {    /* negative x */
 151                 z = x + pio2_1;
 152                 if(ix!=0x3ff921fb) {    /* 33+53 bit pi is good enough */
 153                     y[0] = z + pio2_1t;
 154                     y[1] = (z-y[0])+pio2_1t;
 155                 } else {                /* near pi/2, use 33+33+53 bit pi */
 156                     z += pio2_2;
 157                     y[0] = z + pio2_2t;
 158                     y[1] = (z-y[0])+pio2_2t;
 159                 }
 160                 return -1;
 161             }
 162         }
 163         if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
 164             t  = fd_fabs(x);
 165             n  = (int) (t*invpio2+half);
 166             fn = (double)n;
 167             r  = t-fn*pio2_1;
 168             w  = fn*pio2_1t;    /* 1st round good to 85 bit */
 169             if(n<32&&ix!=npio2_hw[n-1]) {
 170                 y[0] = r-w;     /* quick check no cancellation */
 171             } else {
 172                 j  = ix>>20;
 173                 y[0] = r-w;
 174                 u.d = y[0];
 175                 i = j-(((__HI(u))>>20)&0x7ff);
 176                 if(i>16) {  /* 2nd iteration needed, good to 118 */
 177                     t  = r;
 178                     w  = fn*pio2_2;
 179                     r  = t-w;
 180                     w  = fn*pio2_2t-((t-r)-w);
 181                     y[0] = r-w;
 182                     u.d = y[0];
 183                     i = j-(((__HI(u))>>20)&0x7ff);
 184                     if(i>49)  { /* 3rd iteration need, 151 bits acc */
 185                         t  = r; /* will cover all possible cases */
 186                         w  = fn*pio2_3;
 187                         r  = t-w;
 188                         w  = fn*pio2_3t-((t-r)-w);
 189                         y[0] = r-w;
 190                     }
 191                 }
 192             }
 193             y[1] = (r-y[0])-w;
 194             if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}
 195             else         return n;
 196         }
 197     /*
 198      * all other (large) arguments
 199      */
 200         if(ix>=0x7ff00000) {            /* x is inf or NaN */
 201             y[0]=y[1]=x-x; return 0;
 202         }
 203     /* set z = scalbn(|x|,ilogb(x)-23) */
 204         ux.d = x; uz.d = z;
 205         __LO(uz) = __LO(ux);
 206         z = uz.d;
 207         e0      = (ix>>20)-1046;        /* e0 = ilogb(z)-23; */
 208         uz.d = z;
 209         __HI(uz) = ix - (e0<<20);
 210         z = uz.d;
 211         for(i=0;i<2;i++) {
 212                 tx[i] = (double)((int)(z));
 213                 z     = (z-tx[i])*two24;
 214         }
 215         tx[2] = z;
 216         nx = 3;
 217         while(tx[nx-1]==zero) nx--;     /* skip zero term */
 218         n  =  __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
 219         if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
 220         return n;
 221 }