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