2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001, 2002, 2004 Free Software Foundation
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 /***************************************************************************/
21 /* MODULE_NAME: upow.c */
28 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h */
29 /* halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c */
31 /* root.tbl uexp.tbl upow.tbl */
32 /* An ultimate power routine. Given two IEEE double machine numbers y,x */
33 /* it computes the correctly rounded (to nearest) value of x^y. */
34 /* Assumption: Machine arithmetic operations are performed in */
35 /* round to nearest mode of IEEE 754 standard. */
37 /***************************************************************************/
44 #include "math_private.h"
47 double __exp1(double x
, double xx
, double error
);
48 static double log1(double x
, double *delta
, double *error
);
49 static double my_log2(double x
, double *delta
, double *error
);
50 double __slowpow(double x
, double y
,double z
);
51 static double power1(double x
, double y
);
52 static int checkint(double x
);
54 /***************************************************************************/
55 /* An ultimate power routine. Given two IEEE double machine numbers y,x */
56 /* it computes the correctly rounded (to nearest) value of X^y. */
57 /***************************************************************************/
58 double __ieee754_pow(double x
, double y
) {
59 double z
,a
,aa
,error
, t
,a1
,a2
,y1
,y2
;
68 if (v
.i
[LOW_HALF
] == 0) { /* of y */
69 qx
= u
.i
[HIGH_HALF
]&0x7fffffff;
70 /* Checking if x is not too small to compute */
71 if (((qx
==0x7ff00000)&&(u
.i
[LOW_HALF
]!=0))||(qx
>0x7ff00000)) return NaNQ
.x
;
72 if (y
== 1.0) return x
;
73 if (y
== 2.0) return x
*x
;
74 if (y
== -1.0) return 1.0/x
;
75 if (y
== 0) return 1.0;
78 if(((u
.i
[HIGH_HALF
]>0 && u
.i
[HIGH_HALF
]<0x7ff00000)|| /* x>0 and not x->0 */
79 (u
.i
[HIGH_HALF
]==0 && u
.i
[LOW_HALF
]!=0)) &&
80 /* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
81 (v
.i
[HIGH_HALF
]&0x7fffffff) < 0x4ff00000) { /* if y<-1 or y>1 */
82 z
= log1(x
,&aa
,&error
); /* x^y =e^(y log (X)) */
94 t
= __exp1(a1
,a2
,1.9e16
*error
); /* return -10 or 0 if wasn't computed exactly */
95 return (t
>0)?t
:power1(x
,y
);
99 if (((v
.i
[HIGH_HALF
] & 0x7fffffff) == 0x7ff00000 && v
.i
[LOW_HALF
] != 0)
100 || (v
.i
[HIGH_HALF
] & 0x7fffffff) > 0x7ff00000)
102 if (ABS(y
) > 1.0e20
) return (y
>0)?0:INF
.x
;
105 return y
< 0 ? 1.0/x
: x
;
107 return y
< 0 ? 1.0/ABS(x
) : 0.0; /* return 0 */
110 if (u
.i
[HIGH_HALF
] < 0) {
113 if ((v
.i
[HIGH_HALF
] & 0x7fffffff) == 0x7ff00000 && v
.i
[LOW_HALF
] == 0) {
114 if (x
== -1.0) return 1.0;
115 else if (x
> -1.0) return v
.i
[HIGH_HALF
] < 0 ? INF
.x
: 0.0;
116 else return v
.i
[HIGH_HALF
] < 0 ? 0.0 : INF
.x
;
118 else if (u
.i
[HIGH_HALF
] == 0xfff00000 && u
.i
[LOW_HALF
] == 0)
119 return y
< 0 ? 0.0 : INF
.x
;
120 return NaNQ
.x
; /* y not integer and x<0 */
122 else if (u
.i
[HIGH_HALF
] == 0xfff00000 && u
.i
[LOW_HALF
] == 0)
125 return y
< 0 ? nZERO
.x
: nINF
.x
;
127 return y
< 0 ? 0.0 : INF
.x
;
129 return (k
==1)?__ieee754_pow(-x
,y
):-__ieee754_pow(-x
,y
); /* if y even or odd */
132 qx
= u
.i
[HIGH_HALF
]&0x7fffffff; /* no sign */
133 qy
= v
.i
[HIGH_HALF
]&0x7fffffff; /* no sign */
135 if (qx
> 0x7ff00000 || (qx
== 0x7ff00000 && u
.i
[LOW_HALF
] != 0)) return NaNQ
.x
;
136 /* if 0<x<2^-0x7fe */
137 if (qy
> 0x7ff00000 || (qy
== 0x7ff00000 && v
.i
[LOW_HALF
] != 0))
138 return x
== 1.0 ? 1.0 : NaNQ
.x
;
141 if (qx
== 0x7ff00000) /* x= 2^-0x3ff */
142 {if (y
== 0) return NaNQ
.x
;
145 if (qy
> 0x45f00000 && qy
< 0x7ff00000) {
146 if (x
== 1.0) return 1.0;
147 if (y
>0) return (x
>1.0)?INF
.x
:0;
148 if (y
<0) return (x
<1.0)?INF
.x
:0;
151 if (x
== 1.0) return 1.0;
152 if (y
>0) return (x
>1.0)?INF
.x
:0;
153 if (y
<0) return (x
<1.0)?INF
.x
:0;
154 return 0; /* unreachable, to make the compiler happy */
157 /**************************************************************************/
158 /* Computing x^y using more accurate but more slow log routine */
159 /**************************************************************************/
160 static double power1(double x
, double y
) {
161 double z
,a
,aa
,error
, t
,a1
,a2
,y1
,y2
;
162 z
= my_log2(x
,&aa
,&error
);
170 aa
= ((y1
*a1
-a
)+y1
*a2
+y2
*a1
)+y2
*a2
+aa
*y
;
173 error
= error
*ABS(y
);
174 t
= __exp1(a1
,a2
,1.9e16
*error
);
175 return (t
>= 0)?t
:__slowpow(x
,y
,z
);
178 /****************************************************************************/
179 /* Computing log(x) (x is left argument). The result is the returned double */
180 /* + the parameter delta. */
181 /* The result is bounded by error (rightmost argument) */
182 /****************************************************************************/
183 static double log1(double x
, double *delta
, double *error
) {
188 double uu
,vv
,eps
,nx
,e
,e1
,e2
,t
,t1
,t2
,res
,add
=0;
195 /**/ two52
= {{0x43300000, 0x00000000}}; /* 2**52 */
199 /**/ two52
= {{0x00000000, 0x43300000}}; /* 2**52 */
207 if (m
< 0x00100000) /* 1<x<2^-1007 */
208 { x
= x
*t52
.x
; add
= -52.0; u
.x
= x
; m
= u
.i
[HIGH_HALF
];}
210 if ((m
&0x000fffff) < 0x0006a09e)
211 {u
.i
[HIGH_HALF
] = (m
&0x000fffff)|0x3ff00000; two52
.i
[LOW_HALF
]=(m
>>20); }
213 {u
.i
[HIGH_HALF
] = (m
&0x000fffff)|0x3fe00000; two52
.i
[LOW_HALF
]=(m
>>20)+1; }
217 i
= (v
.i
[LOW_HALF
]&0x000003ff)<<2;
218 if (two52
.i
[LOW_HALF
] == 1023) /* nx = 0 */
220 if (i
> 1192 && i
< 1208) /* |x-1| < 1.5*2**-10 */
223 t1
= (t
+5.0e6
)-5.0e6
;
226 e2
= t
*t
*t
*(r3
+t
*(r4
+t
*(r5
+t
*(r6
+t
*(r7
+t
*r8
)))))-0.5*t2
*(t
+t1
);
228 *error
= 1.0e-21*ABS(t
);
229 *delta
= (e1
-res
)+e2
;
231 } /* |x-1| < 1.5*2**-10 */
234 v
.x
= u
.x
*(ui
.x
[i
]+ui
.x
[i
+1])+bigv
.x
;
236 j
= v
.i
[LOW_HALF
]&0x0007ffff;
240 e2
= eps
*(ui
.x
[i
+1]+vj
.x
[j
]*(ui
.x
[i
]+ui
.x
[i
+1]));
243 t
=ui
.x
[i
+2]+vj
.x
[j
+1];
245 t2
= (((t
-t1
)+e
)+(ui
.x
[i
+3]+vj
.x
[j
+2]))+e2
+e
*e
*(p2
+e
*(p3
+e
*p4
));
248 *delta
= (t1
-res
)+t2
;
255 nx
= (two52
.x
- two52e
.x
)+add
;
260 t
=nx
*ln2a
.x
+ui
.x
[i
+2];
262 t2
=(((t
-t1
)+e
)+nx
*ln2b
.x
+ui
.x
[i
+3]+e2
)+e
*e
*(q2
+e
*(q3
+e
*(q4
+e
*(q5
+e
*q6
))));
265 *delta
= (t1
-res
)+t2
;
270 /****************************************************************************/
271 /* More slow but more accurate routine of log */
272 /* Computing log(x)(x is left argument).The result is return double + delta.*/
273 /* The result is bounded by error (right argument) */
274 /****************************************************************************/
275 static double my_log2(double x
, double *delta
, double *error
) {
280 double uu
,vv
,eps
,nx
,e
,e1
,e2
,t
,t1
,t2
,res
,add
=0;
284 double ou1
,ou2
,lu1
,lu2
,ov
,lv1
,lv2
,a
,a1
,a2
;
285 double y
,yy
,z
,zz
,j1
,j2
,j3
,j4
,j5
,j6
,j7
,j8
;
289 /**/ two52
= {{0x43300000, 0x00000000}}; /* 2**52 */
293 /**/ two52
= {{0x00000000, 0x43300000}}; /* 2**52 */
302 if (m
<0x00100000) { /* x < 2^-1022 */
303 x
= x
*t52
.x
; add
= -52.0; u
.x
= x
; m
= u
.i
[HIGH_HALF
]; }
305 if ((m
&0x000fffff) < 0x0006a09e)
306 {u
.i
[HIGH_HALF
] = (m
&0x000fffff)|0x3ff00000; two52
.i
[LOW_HALF
]=(m
>>20); }
308 {u
.i
[HIGH_HALF
] = (m
&0x000fffff)|0x3fe00000; two52
.i
[LOW_HALF
]=(m
>>20)+1; }
312 i
= (v
.i
[LOW_HALF
]&0x000003ff)<<2;
313 /*------------------------------------- |x-1| < 2**-11------------------------------- */
314 if ((two52
.i
[LOW_HALF
] == 1023) && (i
== 1200))
317 EMULV(t
,s3
,y
,yy
,j1
,j2
,j3
,j4
,j5
);
318 ADD2(-0.5,0,y
,yy
,z
,zz
,j1
,j2
);
319 MUL2(t
,0,z
,zz
,y
,yy
,j1
,j2
,j3
,j4
,j5
,j6
,j7
,j8
);
320 MUL2(t
,0,y
,yy
,z
,zz
,j1
,j2
,j3
,j4
,j5
,j6
,j7
,j8
);
323 e2
= (((t
-e1
)+z
)+zz
)+t
*t
*t
*(ss3
+t
*(s4
+t
*(s5
+t
*(s6
+t
*(s7
+t
*s8
)))));
325 *error
= 1.0e-25*ABS(t
);
326 *delta
= (e1
-res
)+e2
;
329 /*----------------------------- |x-1| > 2**-11 -------------------------- */
331 { /*Computing log(x) according to log table */
332 nx
= (two52
.x
- two52e
.x
)+add
;
337 v
.x
= u
.x
*(ou1
+ou2
)+bigv
.x
;
339 j
= v
.i
[LOW_HALF
]&0x0007ffff;
345 a
= (ou1
+ou2
)*(1.0+ov
);
346 a1
= (a
+1.0e10
)-1.0e10
;
347 a2
= a
*(1.0-a1
*uu
*vv
);
354 t2
= (((t
-t1
)+e
)+(lu2
+lv2
+nx
*ln2b
.x
+e2
))+e
*e
*(p2
+e
*(p3
+e
*p4
));
357 *delta
= (t1
-res
)+t2
;
362 /**********************************************************************/
363 /* Routine receives a double x and checks if it is an integer. If not */
364 /* it returns 0, else it returns 1 if even or -1 if odd. */
365 /**********************************************************************/
366 static int checkint(double x
) {
367 union {int4 i
[2]; double x
;} u
;
373 m
= u
.i
[HIGH_HALF
]&0x7fffffff; /* no sign */
374 if (m
>= 0x7ff00000) return 0; /* x is +/-inf or NaN */
375 if (m
>= 0x43400000) return 1; /* |x| >= 2**53 */
376 if (m
< 0x40000000) return 0; /* |x| < 2, can not be 0 or 1 */
378 k
= (m
>>20)-1023; /* 1 <= k <= 52 */
379 if (k
== 52) return (n
&1)? -1:1; /* odd or even*/
381 if (n
<<(k
-20)) return 0; /* if not integer */
382 return (n
<<(k
-21))?-1:1;
384 if (n
) return 0; /*if not integer*/
385 if (k
== 20) return (m
&1)? -1:1;
386 if (m
<<(k
+12)) return 0;
387 return (m
<<(k
+11))?-1:1;