Back to index

lightning-sunbird  0.9+nobinonly
e_rem_pio2.c
Go to the documentation of this file.
00001 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
00002  *
00003  * ***** BEGIN LICENSE BLOCK *****
00004  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
00005  *
00006  * The contents of this file are subject to the Mozilla Public License Version
00007  * 1.1 (the "License"); you may not use this file except in compliance with
00008  * the License. You may obtain a copy of the License at
00009  * http://www.mozilla.org/MPL/
00010  *
00011  * Software distributed under the License is distributed on an "AS IS" basis,
00012  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
00013  * for the specific language governing rights and limitations under the
00014  * License.
00015  *
00016  * The Original Code is Mozilla Communicator client code, released
00017  * March 31, 1998.
00018  *
00019  * The Initial Developer of the Original Code is
00020  * Sun Microsystems, Inc.
00021  * Portions created by the Initial Developer are Copyright (C) 1998
00022  * the Initial Developer. All Rights Reserved.
00023  *
00024  * Contributor(s):
00025  *
00026  * Alternatively, the contents of this file may be used under the terms of
00027  * either of the GNU General Public License Version 2 or later (the "GPL"),
00028  * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00029  * in which case the provisions of the GPL or the LGPL are applicable instead
00030  * of those above. If you wish to allow use of your version of this file only
00031  * under the terms of either the GPL or the LGPL, and not to allow others to
00032  * use your version of this file under the terms of the MPL, indicate your
00033  * decision by deleting the provisions above and replace them with the notice
00034  * and other provisions required by the GPL or the LGPL. If you do not delete
00035  * the provisions above, a recipient may use your version of this file under
00036  * the terms of any one of the MPL, the GPL or the LGPL.
00037  *
00038  * ***** END LICENSE BLOCK ***** */
00039 
00040 /* @(#)e_rem_pio2.c 1.4 95/01/18 */
00041 /*
00042  * ====================================================
00043  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00044  *
00045  * Developed at SunSoft, a Sun Microsystems, Inc. business.
00046  * Permission to use, copy, modify, and distribute this
00047  * software is freely granted, provided that this notice 
00048  * is preserved.
00049  * ====================================================
00050  *
00051  */
00052 
00053 /* __ieee754_rem_pio2(x,y)
00054  * 
00055  * return the remainder of x rem pi/2 in y[0]+y[1] 
00056  * use __kernel_rem_pio2()
00057  */
00058 
00059 #include "fdlibm.h"
00060 
00061 /*
00062  * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 
00063  */
00064 #ifdef __STDC__
00065 static const int two_over_pi[] = {
00066 #else
00067 static int two_over_pi[] = {
00068 #endif
00069 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 
00070 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 
00071 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 
00072 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 
00073 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 
00074 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 
00075 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 
00076 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 
00077 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 
00078 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 
00079 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 
00080 };
00081 
00082 #ifdef __STDC__
00083 static const int npio2_hw[] = {
00084 #else
00085 static int npio2_hw[] = {
00086 #endif
00087 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
00088 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
00089 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
00090 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
00091 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
00092 0x404858EB, 0x404921FB,
00093 };
00094 
00095 /*
00096  * invpio2:  53 bits of 2/pi
00097  * pio2_1:   first  33 bit of pi/2
00098  * pio2_1t:  pi/2 - pio2_1
00099  * pio2_2:   second 33 bit of pi/2
00100  * pio2_2t:  pi/2 - (pio2_1+pio2_2)
00101  * pio2_3:   third  33 bit of pi/2
00102  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
00103  */
00104 
00105 #ifdef __STDC__
00106 static const double 
00107 #else
00108 static double 
00109 #endif
00110 zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
00111 half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
00112 two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
00113 invpio2 =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
00114 pio2_1  =  1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
00115 pio2_1t =  6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
00116 pio2_2  =  6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
00117 pio2_2t =  2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
00118 pio2_3  =  2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
00119 pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
00120 
00121 #ifdef __STDC__
00122        int __ieee754_rem_pio2(double x, double *y)
00123 #else
00124        int __ieee754_rem_pio2(x,y)
00125        double x,y[];
00126 #endif
00127 {
00128         fd_twoints u, ux, uz;
00129         double z = 0;
00130        double w,t,r,fn;
00131        double tx[3];
00132        int e0,i,j,nx,n,ix,hx;
00133 
00134         u.d = x;
00135        hx = __HI(u);        /* high word of x */
00136        ix = hx&0x7fffffff;
00137        if(ix<=0x3fe921fb)   /* |x| ~<= pi/4 , no need for reduction */
00138            {y[0] = x; y[1] = 0; return 0;}
00139        if(ix<0x4002d97c) {  /* |x| < 3pi/4, special case with n=+-1 */
00140            if(hx>0) { 
00141               z = x - pio2_1;
00142               if(ix!=0x3ff921fb) {        /* 33+53 bit pi is good enough */
00143                   y[0] = z - pio2_1t;
00144                   y[1] = (z-y[0])-pio2_1t;
00145               } else {             /* near pi/2, use 33+33+53 bit pi */
00146                   z -= pio2_2;
00147                   y[0] = z - pio2_2t;
00148                   y[1] = (z-y[0])-pio2_2t;
00149               }
00150               return 1;
00151            } else {  /* negative x */
00152               z = x + pio2_1;
00153               if(ix!=0x3ff921fb) {        /* 33+53 bit pi is good enough */
00154                   y[0] = z + pio2_1t;
00155                   y[1] = (z-y[0])+pio2_1t;
00156               } else {             /* near pi/2, use 33+33+53 bit pi */
00157                   z += pio2_2;
00158                   y[0] = z + pio2_2t;
00159                   y[1] = (z-y[0])+pio2_2t;
00160               }
00161               return -1;
00162            }
00163        }
00164        if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
00165            t  = fd_fabs(x);
00166            n  = (int) (t*invpio2+half);
00167            fn = (double)n;
00168            r  = t-fn*pio2_1;
00169            w  = fn*pio2_1t; /* 1st round good to 85 bit */
00170            if(n<32&&ix!=npio2_hw[n-1]) {  
00171               y[0] = r-w;   /* quick check no cancellation */
00172            } else {
00173                j  = ix>>20;
00174                y[0] = r-w;
00175                 u.d = y[0];
00176                i = j-(((__HI(u))>>20)&0x7ff);
00177                if(i>16) {  /* 2nd iteration needed, good to 118 */
00178                   t  = r;
00179                   w  = fn*pio2_2;  
00180                   r  = t-w;
00181                   w  = fn*pio2_2t-((t-r)-w);     
00182                   y[0] = r-w;
00183                     u.d = y[0];
00184                   i = j-(((__HI(u))>>20)&0x7ff);
00185                   if(i>49)  {      /* 3rd iteration need, 151 bits acc */
00186                      t  = r;       /* will cover all possible cases */
00187                      w  = fn*pio2_3;      
00188                      r  = t-w;
00189                      w  = fn*pio2_3t-((t-r)-w);  
00190                      y[0] = r-w;
00191                   }
00192               }
00193            }
00194            y[1] = (r-y[0])-w;
00195            if(hx<0)  {y[0] = -y[0]; y[1] = -y[1]; return -n;}
00196            else       return n;
00197        }
00198     /* 
00199      * all other (large) arguments
00200      */
00201        if(ix>=0x7ff00000) {        /* x is inf or NaN */
00202            y[0]=y[1]=x-x; return 0;
00203        }
00204     /* set z = scalbn(|x|,ilogb(x)-23) */
00205         ux.d = x; uz.d = z;
00206        __LO(uz) = __LO(ux);
00207         z = uz.d;
00208        e0     = (ix>>20)-1046;     /* e0 = ilogb(z)-23; */
00209         uz.d = z;
00210        __HI(uz) = ix - (e0<<20);
00211         z = uz.d;
00212        for(i=0;i<2;i++) {
00213               tx[i] = (double)((int)(z));
00214               z     = (z-tx[i])*two24;
00215        }
00216        tx[2] = z;
00217        nx = 3;
00218        while(tx[nx-1]==zero) nx--; /* skip zero term */
00219        n  =  __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
00220        if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
00221        return n;
00222 }