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