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