Commit | Line | Data |
---|---|---|
94759b7b TL |
1 | /* |
2 | sqrt returns the square root of its floating | |
3 | point argument. Newton's method. | |
4 | ||
5 | calls frexp | |
6 | */ | |
7 | ||
8 | #include <errno.h> | |
9 | ||
10 | int errno; | |
11 | double frexp(); | |
12 | ||
13 | double | |
14 | sqrt(arg) | |
15 | double arg; | |
16 | { | |
17 | double x, temp; | |
18 | int exp; | |
19 | int i; | |
20 | ||
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 | } |