1 /*
2 	sqrt returns the square root of its floating
3 	point argument. Newton's method.
4 
5 	calls frexp
6 */
7 
8 #include <u.h>
9 #include <libc.h>
10 
11 double
12 sqrt(double arg)
13 {
14 	double x, temp;
15 	int exp, i;
16 
17 	if(arg <= 0) {
18 		if(arg < 0)
19 			return 0.;
20 		return 0;
21 	}
22 	x = frexp(arg, &exp);
23 	while(x < 0.5) {
24 		x *= 2;
25 		exp--;
26 	}
27 	/*
28 	 * NOTE
29 	 * this wont work on 1's comp
30 	 */
31 	if(exp & 1) {
32 		x *= 2;
33 		exp--;
34 	}
35 	temp = 0.5 * (1.0+x);
36 
37 	while(exp > 60) {
38 		temp *= (1L<<30);
39 		exp -= 60;
40 	}
41 	while(exp < -60) {
42 		temp /= (1L<<30);
43 		exp += 60;
44 	}
45 	if(exp >= 0)
46 		temp *= 1L << (exp/2);
47 	else
48 		temp /= 1L << (-exp/2);
49 	for(i=0; i<=4; i++)
50 		temp = 0.5*(temp + arg/temp);
51 	return temp;
52 }