Subversion Repositories planix.SVN

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
#include <u.h>
2 
#include <libc.h>
3 
 
4 
double
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 
	flip = 0;
15 
	if(y < 0.){
16 
		y = -y;
17 
		flip = 1;
18 
	}
19 
	y1 = modf(y, &ye);
20 
	if(y1 != 0.0){
21 
		if(x <= 0.)
22 
			goto zreturn;
23 
		if(y1 > 0.5) {
24 
			y1 -= 1.;
25 
			ye += 1.;
26 
		}
27 
		xy = exp(y1 * log(x));
28 
	}else
29 
		xy = 1.0;
30 
	if(ye > 0x7FFFFFFF){	/* should be ~0UL but compiler can't convert double to ulong */
31 
		if(x <= 0.){
32 
 zreturn:
33 
			if(x==0. && !flip)
34 
				return 0.;
35 
			return NaN();
36 
		}
37 
		if(flip){
38 
			if(y == .5)
39 
				return 1./sqrt(x);
40 
			y = -y;
41 
		}else if(y == .5)
42 
			return sqrt(x);
43 
		return exp(y * log(x));
44 
	}
45 
	x = frexp(x, &ex);
46 
	ey = 0;
47 
	i = ye;
48 
	if(i)
49 
		for(;;){
50 
			if(i & 1){
51 
				xy *= x;
52 
				ey += ex;
53 
			}
54 
			i >>= 1;
55 
			if(i == 0)
56 
				break;
57 
			x *= x;
58 
			ex <<= 1;
59 
			if(x < .5){
60 
				x += x;
61 
				ex -= 1;
62 
			}
63 
		}
64 
	if(flip){
65 
		xy = 1. / xy;
66 
		ey = -ey;
67 
	}
68 
	return ldexp(xy, ey);
69 
}