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