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 <math.h>
9 
#include <errno.h>
10 
#define _RESEARCH_SOURCE
11 
#include <float.h>
12 
 
13 
double
14 
sqrt(double arg)
15 
{
16 
	double x, temp;
17 
	int exp, i;
18 
 
19 
	if(isInf(arg, 1))
20 
		return arg;
21 
	if(arg <= 0) {
22 
		if(arg < 0)
23 
			errno = EDOM;
24 
		return 0;
25 
	}
26 
	x = frexp(arg, &exp);
27 
	while(x < 0.5) {
28 
		x *= 2;
29 
		exp--;
30 
	}
31 
	/*
32 
	 * NOTE
33 
	 * this wont work on 1's comp
34 
	 */
35 
	if(exp & 1) {
36 
		x *= 2;
37 
		exp--;
38 
	}
39 
	temp = 0.5 * (1.0+x);
40 
 
41 
	while(exp > 60) {
42 
		temp *= (1L<<30);
43 
		exp -= 60;
44 
	}
45 
	while(exp < -60) {
46 
		temp /= (1L<<30);
47 
		exp += 60;
48 
	}
49 
	if(exp >= 0)
50 
		temp *= 1L << (exp/2);
51 
	else
52 
		temp /= 1L << (-exp/2);
53 
	for(i=0; i<=4; i++)
54 
		temp = 0.5*(temp + arg/temp);
55 
	return temp;
56 
}