More namespace issues.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / dbl-64 / e_pow.c
1 /*
2  * IBM Accurate Mathematical Library
3  * Copyright (c) International Business Machines Corp., 2001
4  *
5  * This program is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU Lesser General Public License as published by
7  * the Free Software Foundation; either version 2 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public License
16  * along with this program; if not, write to the Free Software
17  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18  */
19 /***************************************************************************/
20 /*  MODULE_NAME: upow.c                                                    */
21 /*                                                                         */
22 /*  FUNCTIONS: upow                                                        */
23 /*             power1                                                      */
24 /*             log2                                                        */
25 /*             log1                                                        */
26 /*             checkint                                                    */
27 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h                             */
28 /*               halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c       */
29 /*                          uexp.c  upow.c                                 */
30 /*               root.tbl uexp.tbl upow.tbl                                */
31 /* An ultimate power routine. Given two IEEE double machine numbers y,x    */
32 /* it computes the correctly rounded (to nearest) value of x^y.            */
33 /* Assumption: Machine arithmetic operations are performed in              */
34 /* round to nearest mode of IEEE 754 standard.                             */
35 /*                                                                         */
36 /***************************************************************************/
37 #include "endian.h"
38 #include "upow.h"
39 #include "dla.h"
40 #include "mydefs.h"
41 #include "MathLib.h"
42 #include "upow.tbl"
43
44
45 double __exp1(double x, double xx, double error);
46 static double log1(double x, double *delta, double *error);
47 static double log2(double x, double *delta, double *error);
48 double slowpow(double x, double y,double z);
49 static double power1(double x, double y);
50 static int checkint(double x);
51
52 /***************************************************************************/
53 /* An ultimate power routine. Given two IEEE double machine numbers y,x    */
54 /* it computes the correctly rounded (to nearest) value of X^y.            */
55 /***************************************************************************/
56 double __ieee754_pow(double x, double y) {
57   double z,a,aa,error, t,a1,a2,y1,y2;
58 #if 0
59   double gor=1.0;
60 #endif
61   mynumber u,v;
62   int k;
63   int4 qx,qy;
64   v.x=y;
65   u.x=x;
66   if (v.i[LOW_HALF] == 0) { /* of y */
67     qx = u.i[HIGH_HALF]&0x7fffffff;
68     /* Checking  if x is not too small to compute */
69     if (((qx==0x7ff00000)&&(u.i[LOW_HALF]!=0))||(qx>0x7ff00000)) return NaNQ.x;
70     if (y == 1.0) return x;
71     if (y == 2.0) return x*x;
72     if (y == -1.0) return (x!=0)?1.0/x:NaNQ.x;
73     if (y == 0) return ((x>0)&&(qx<0x7ff00000))?1.0:NaNQ.x;
74   }
75   /* else */
76   if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)||        /* x>0 and not x->0 */
77        (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0))  &&
78                                       /*   2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
79       (v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) {              /* if y<-1 or y>1   */
80     z = log1(x,&aa,&error);                                 /* x^y  =e^(y log (X)) */
81     t = y*134217729.0;
82     y1 = t - (t-y);
83     y2 = y - y1;
84     t = z*134217729.0;
85     a1 = t - (t-z);
86     a2 = (z - a1)+aa;
87     a = y1*a1;
88     aa = y2*a1 + y*a2;
89     a1 = a+aa;
90     a2 = (a-a1)+aa;
91     error = error*ABS(y);
92     t = __exp1(a1,a2,1.9e16*error);     /* return -10 or 0 if wasn't computed exactly */
93     return (t>0)?t:power1(x,y);
94   }
95
96   if (x == 0) {
97     if (ABS(y) > 1.0e20) return (y>0)?0:NaNQ.x;
98     k = checkint(y);
99     if (k == 0 || y < 0) return NaNQ.x;
100     else return (k==1)?0:x;                                            /* return 0 */
101   }
102   /* if x<0 */
103   if (u.i[HIGH_HALF] < 0) {
104     k = checkint(y);
105     if (k==0) return NaNQ.x;                              /* y not integer and x<0 */
106     return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */
107   }
108   /* x>0 */
109   qx = u.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
110   qy = v.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
111
112   if (qx > 0x7ff00000 || (qx == 0x7ff00000 && u.i[LOW_HALF] != 0)) return NaNQ.x;
113                                                                  /*  if 0<x<2^-0x7fe */
114   if (qy > 0x7ff00000 || (qy == 0x7ff00000 && v.i[LOW_HALF] != 0)) return NaNQ.x;
115                                                                  /*  if y<2^-0x7fe   */
116
117   if (qx == 0x7ff00000)                              /* x= 2^-0x3ff */
118     {if (y == 0) return NaNQ.x;
119     return (y>0)?x:0; }
120
121   if (qy > 0x45f00000 && qy < 0x7ff00000) {
122     if (x == 1.0) return 1.0;
123     if (y>0) return (x>1.0)?INF.x:0;
124     if (y<0) return (x<1.0)?INF.x:0;
125   }
126
127   if (x == 1.0) return NaNQ.x;
128   if (y>0) return (x>1.0)?INF.x:0;
129   if (y<0) return (x<1.0)?INF.x:0;
130   return 0;     /* unreachable, to make the compiler happy */
131 }
132
133 /**************************************************************************/
134 /* Computing x^y using more accurate but more slow log routine            */
135 /**************************************************************************/
136 static double power1(double x, double y) {
137   double z,a,aa,error, t,a1,a2,y1,y2;
138   z = log2(x,&aa,&error);
139   t = y*134217729.0;
140   y1 = t - (t-y);
141   y2 = y - y1;
142   t = z*134217729.0;
143   a1 = t - (t-z);
144   a2 = z - a1;
145   a = y*z;
146   aa = ((y1*a1-a)+y1*a2+y2*a1)+y2*a2+aa*y;
147   a1 = a+aa;
148   a2 = (a-a1)+aa;
149   error = error*ABS(y);
150   t = __exp1(a1,a2,1.9e16*error);
151   return (t >= 0)?t:slowpow(x,y,z);
152 }
153
154 /****************************************************************************/
155 /* Computing log(x) (x is left argument). The result is the returned double */
156 /* + the parameter delta.                                                   */
157 /* The result is bounded by error (rightmost argument)                      */
158 /****************************************************************************/
159 static double log1(double x, double *delta, double *error) {
160   int i,j,m;
161 #if 0
162   int n;
163 #endif
164   double uu,vv,eps,nx,e,e1,e2,t,t1,t2,res,add=0;
165 #if 0
166   double cor;
167 #endif
168   mynumber u,v;
169
170   u.x = x;
171   m = u.i[HIGH_HALF];
172   *error = 0;
173   *delta = 0;
174   if (m < 0x00100000)             /*  1<x<2^-1007 */
175     { x = x*t52.x; add = -52.0; u.x = x; m = u.i[HIGH_HALF];}
176
177   if ((m&0x000fffff) < 0x0006a09e)
178     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3ff00000; two52.i[LOW_HALF]=(m>>20); }
179   else
180     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3fe00000; two52.i[LOW_HALF]=(m>>20)+1; }
181
182   v.x = u.x + bigu.x;
183   uu = v.x - bigu.x;
184   i = (v.i[LOW_HALF]&0x000003ff)<<2;
185   if (two52.i[LOW_HALF] == 1023)         /* nx = 0              */
186   {
187       if (i > 1192 && i < 1208)          /* |x-1| < 1.5*2**-10  */
188       {
189           t = x - 1.0;
190           t1 = (t+5.0e6)-5.0e6;
191           t2 = t-t1;
192           e1 = t - 0.5*t1*t1;
193           e2 = t*t*t*(r3+t*(r4+t*(r5+t*(r6+t*(r7+t*r8)))))-0.5*t2*(t+t1);
194           res = e1+e2;
195           *error = 1.0e-21*ABS(t);
196           *delta = (e1-res)+e2;
197           return res;
198       }                  /* |x-1| < 1.5*2**-10  */
199       else
200       {
201           v.x = u.x*(ui.x[i]+ui.x[i+1])+bigv.x;
202           vv = v.x-bigv.x;
203           j = v.i[LOW_HALF]&0x0007ffff;
204           j = j+j+j;
205           eps = u.x - uu*vv;
206           e1 = eps*ui.x[i];
207           e2 = eps*(ui.x[i+1]+vj.x[j]*(ui.x[i]+ui.x[i+1]));
208           e = e1+e2;
209           e2 =  ((e1-e)+e2);
210           t=ui.x[i+2]+vj.x[j+1];
211           t1 = t+e;
212           t2 = (((t-t1)+e)+(ui.x[i+3]+vj.x[j+2]))+e2+e*e*(p2+e*(p3+e*p4));
213           res=t1+t2;
214           *error = 1.0e-24;
215           *delta = (t1-res)+t2;
216           return res;
217       }
218   }   /* nx = 0 */
219   else                            /* nx != 0   */
220   {
221       eps = u.x - uu;
222       nx = (two52.x - two52e.x)+add;
223       e1 = eps*ui.x[i];
224       e2 = eps*ui.x[i+1];
225       e=e1+e2;
226       e2 = (e1-e)+e2;
227       t=nx*ln2a.x+ui.x[i+2];
228       t1=t+e;
229       t2=(((t-t1)+e)+nx*ln2b.x+ui.x[i+3]+e2)+e*e*(q2+e*(q3+e*(q4+e*(q5+e*q6))));
230       res = t1+t2;
231       *error = 1.0e-21;
232       *delta = (t1-res)+t2;
233       return res;
234   }                                /* nx != 0   */
235 }
236
237 /****************************************************************************/
238 /* More slow but more accurate routine of log                               */
239 /* Computing log(x)(x is left argument).The result is return double + delta.*/
240 /* The result is bounded by error (right argument)                           */
241 /****************************************************************************/
242 static double log2(double x, double *delta, double *error) {
243   int i,j,m;
244 #if 0
245   int n;
246 #endif
247   double uu,vv,eps,nx,e,e1,e2,t,t1,t2,res,add=0;
248 #if 0
249   double cor;
250 #endif
251   double ou1,ou2,lu1,lu2,ov,lv1,lv2,a,a1,a2;
252   double y,yy,z,zz,j1,j2,j3,j4,j5,j6,j7,j8;
253   mynumber u,v;
254
255   u.x = x;
256   m = u.i[HIGH_HALF];
257   *error = 0;
258   *delta = 0;
259   add=0;
260   if (m<0x00100000) {  /* x < 2^-1022 */
261     x = x*t52.x;  add = -52.0; u.x = x; m = u.i[HIGH_HALF]; }
262
263   if ((m&0x000fffff) < 0x0006a09e)
264     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3ff00000; two52.i[LOW_HALF]=(m>>20); }
265   else
266     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3fe00000; two52.i[LOW_HALF]=(m>>20)+1; }
267
268   v.x = u.x + bigu.x;
269   uu = v.x - bigu.x;
270   i = (v.i[LOW_HALF]&0x000003ff)<<2;
271   /*------------------------------------- |x-1| < 2**-11-------------------------------  */
272   if ((two52.i[LOW_HALF] == 1023)  && (i == 1200))
273   {
274       t = x - 1.0;
275       EMULV(t,s3,y,yy,j1,j2,j3,j4,j5);
276       ADD2(-0.5,0,y,yy,z,zz,j1,j2);
277       MUL2(t,0,z,zz,y,yy,j1,j2,j3,j4,j5,j6,j7,j8);
278       MUL2(t,0,y,yy,z,zz,j1,j2,j3,j4,j5,j6,j7,j8);
279
280       e1 = t+z;
281       e2 = (((t-e1)+z)+zz)+t*t*t*(ss3+t*(s4+t*(s5+t*(s6+t*(s7+t*s8)))));
282       res = e1+e2;
283       *error = 1.0e-25*ABS(t);
284       *delta = (e1-res)+e2;
285       return res;
286   }
287   /*----------------------------- |x-1| > 2**-11  --------------------------  */
288   else
289   {          /*Computing log(x) according to log table                        */
290       nx = (two52.x - two52e.x)+add;
291       ou1 = ui.x[i];
292       ou2 = ui.x[i+1];
293       lu1 = ui.x[i+2];
294       lu2 = ui.x[i+3];
295       v.x = u.x*(ou1+ou2)+bigv.x;
296       vv = v.x-bigv.x;
297       j = v.i[LOW_HALF]&0x0007ffff;
298       j = j+j+j;
299       eps = u.x - uu*vv;
300       ov  = vj.x[j];
301       lv1 = vj.x[j+1];
302       lv2 = vj.x[j+2];
303       a = (ou1+ou2)*(1.0+ov);
304       a1 = (a+1.0e10)-1.0e10;
305       a2 = a*(1.0-a1*uu*vv);
306       e1 = eps*a1;
307       e2 = eps*a2;
308       e = e1+e2;
309       e2 = (e1-e)+e2;
310       t=nx*ln2a.x+lu1+lv1;
311       t1 = t+e;
312       t2 = (((t-t1)+e)+(lu2+lv2+nx*ln2b.x+e2))+e*e*(p2+e*(p3+e*p4));
313       res=t1+t2;
314       *error = 1.0e-27;
315       *delta = (t1-res)+t2;
316       return res;
317   }
318 }
319
320 /**********************************************************************/
321 /* Routine receives a double x and checks if it is an integer. If not */
322 /* it returns 0, else it returns 1 if even or -1 if odd.              */
323 /**********************************************************************/
324 static int checkint(double x) {
325   union {int4 i[2]; double x;} u;
326   int k,m,n;
327 #if 0
328   int l;
329 #endif
330   u.x = x;
331   m = u.i[HIGH_HALF]&0x7fffffff;    /* no sign */
332   if (m >= 0x7ff00000) return 0;    /*  x is +/-inf or NaN  */
333   if (m >= 0x43400000) return 1;    /*  |x| >= 2**53   */
334   if (m < 0x40000000) return 0;     /* |x| < 2,  can not be 0 or 1  */
335   n = u.i[LOW_HALF];
336   k = (m>>20)-1023;                 /*  1 <= k <= 52   */
337   if (k == 52) return (n&1)? -1:1;  /* odd or even*/
338   if (k>20) {
339     if (n<<(k-20)) return 0;        /* if not integer */
340     return (n<<(k-21))?-1:1;
341   }
342   if (n) return 0;                  /*if  not integer*/
343   if (k == 20) return (m&1)? -1:1;
344   if (m<<(k+12)) return 0;
345   return (m<<(k+11))?-1:1;
346 }