1 #include <math.h>
2 #include <errno.h>
3 #include <limits.h>
4
5 double
pow(double x,double y)6 pow(double x, double y) /* return x ^ y (exponentiation) */
7 {
8 double xy, y1, ye;
9 long i;
10 int ex, ey, flip;
11
12 if(y == 0.0)
13 return 1.0;
14
15 /* prevent infinite loop */
16 if(isNaN(x) || isNaN(y))
17 return NaN();
18 if(isInf(x, 0))
19 return x;
20 if(isInf(y, 0))
21 return x == 0 || x == 1? x: y;
22
23 flip = 0;
24 if(y < 0.){
25 y = -y;
26 flip = 1;
27 }
28 y1 = modf(y, &ye);
29 if(y1 != 0.0){
30 if(x <= 0.)
31 goto zreturn;
32 if(y1 > 0.5) {
33 y1 -= 1.;
34 ye += 1.;
35 }
36 xy = exp(y1 * log(x));
37 }else
38 xy = 1.0;
39 if(ye > LONG_MAX){
40 if(x <= 0.){
41 zreturn:
42 if(x || flip)
43 errno = EDOM;
44 return 0.;
45 }
46 if(flip){
47 if(y == .5)
48 return 1./sqrt(x);
49 y = -y;
50 }else if(y == .5)
51 return sqrt(x);
52 return exp(y * log(x));
53 }
54 x = frexp(x, &ex);
55 ey = 0;
56 i = ye;
57 if(i)
58 for(;;){
59 if(i & 1){
60 xy *= x;
61 ey += ex;
62 }
63 i >>= 1;
64 if(i == 0)
65 break;
66 x *= x;
67 ex <<= 1;
68 if(x < .5){
69 x += x;
70 ex -= 1;
71 }
72 }
73 if(flip){
74 xy = 1. / xy;
75 ey = -ey;
76 }
77 errno = 0;
78 x = ldexp(xy, ey);
79 if(errno && ey < 0) {
80 errno = 0;
81 x = 0.;
82 }
83 return x;
84 }
85