Subversion Repositories planix.SVN

Rev

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

Rev Author Line No. Line
2 - 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 NaN();
20 
		return 0;
21 
	}
22 
	if(isInf(arg, 1))
23 
		return arg;
24 
	x = frexp(arg, &exp);
25 
	while(x < 0.5) {
26 
		x *= 2;
27 
		exp--;
28 
	}
29 
	/*
30 
	 * NOTE
31 
	 * this wont work on 1's comp
32 
	 */
33 
	if(exp & 1) {
34 
		x *= 2;
35 
		exp--;
36 
	}
37 
	temp = 0.5 * (1.0+x);
38 
 
39 
	while(exp > 60) {
40 
		temp *= (1L<<30);
41 
		exp -= 60;
42 
	}
43 
	while(exp < -60) {
44 
		temp /= (1L<<30);
45 
		exp += 60;
46 
	}
47 
	if(exp >= 0)
48 
		temp *= 1L << (exp/2);
49 
	else
50 
		temp /= 1L << (-exp/2);
51 
	for(i=0; i<=4; i++)
52 
		temp = 0.5*(temp + arg/temp);
53 
	return temp;
54 
}