Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
#include <math.h>
2 
#include <errno.h>
3 
#include <limits.h>
4 
 
5 
double
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 
	flip = 0;
16 
	if(y < 0.){
17 
		y = -y;
18 
		flip = 1;
19 
	}
20 
	y1 = modf(y, &ye);
21 
	if(y1 != 0.0){
22 
		if(x <= 0.)
23 
			goto zreturn;
24 
		if(y1 > 0.5) {
25 
			y1 -= 1.;
26 
			ye += 1.;
27 
		}
28 
		xy = exp(y1 * log(x));
29 
	}else
30 
		xy = 1.0;
31 
	if(ye > LONG_MAX){
32 
		if(x <= 0.){
33 
 zreturn:
34 
			if(x || flip)
35 
				errno = EDOM;
36 
			return 0.;
37 
		}
38 
		if(flip){
39 
			if(y == .5)
40 
				return 1./sqrt(x);
41 
			y = -y;
42 
		}else if(y == .5)
43 
			return sqrt(x);
44 
		return exp(y * log(x));
45 
	}
46 
	x = frexp(x, &ex);
47 
	ey = 0;
48 
	i = ye;
49 
	if(i)
50 
		for(;;){
51 
			if(i & 1){
52 
				xy *= x;
53 
				ey += ex;
54 
			}
55 
			i >>= 1;
56 
			if(i == 0)
57 
				break;
58 
			x *= x;
59 
			ex <<= 1;
60 
			if(x < .5){
61 
				x += x;
62 
				ex -= 1;
63 
			}
64 
		}
65 
	if(flip){
66 
		xy = 1. / xy;
67 
		ey = -ey;
68 
	}
69 
	errno = 0;
70 
	x = ldexp(xy, ey);
71 
	if(errno && ey < 0) {
72 
		errno = 0;
73 
		x = 0.;
74 
	}
75 
	return x;
76 
}