| 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 | } |