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 static double Inf(int); 14 static int isInf(double, int); 15 16 double 17 sqrt(double arg) 18 { 19 double x, temp; 20 int exp, i; 21 22 if(isInf(arg, 1)) 23 return arg; 24 if(arg <= 0) { 25 if(arg < 0) 26 errno = EDOM; 27 return 0; 28 } 29 x = frexp(arg, &exp); 30 while(x < 0.5) { 31 x *= 2; 32 exp--; 33 } 34 /* 35 * NOTE 36 * this wont work on 1's comp 37 */ 38 if(exp & 1) { 39 x *= 2; 40 exp--; 41 } 42 temp = 0.5 * (1.0+x); 43 44 while(exp > 60) { 45 temp *= (1L<<30); 46 exp -= 60; 47 } 48 while(exp < -60) { 49 temp /= (1L<<30); 50 exp += 60; 51 } 52 if(exp >= 0) 53 temp *= 1L << (exp/2); 54 else 55 temp /= 1L << (-exp/2); 56 for(i=0; i<=4; i++) 57 temp = 0.5*(temp + arg/temp); 58 return temp; 59 } 60 61 typedef union 62 { 63 double d; 64 struct 65 { 66 #ifdef IEEE_8087 67 long ls; 68 long ms; 69 #else 70 long ms; 71 long ls; 72 #endif 73 }; 74 } Cheat; 75 76 #define NANEXP (2047<<20) 77 #define NANMASK (2047<<20) 78 #define NANSIGN (1<<31) 79 80 static double 81 Inf(int sign) 82 { 83 Cheat a; 84 85 a.ms = NANEXP; 86 a.ls = 0; 87 if(sign < 0) 88 a.ms |= NANSIGN; 89 return a.d; 90 } 91 92 static int 93 isInf(double d, int sign) 94 { 95 Cheat a; 96 97 a.d = d; 98 if(a.ls != 0) 99 return 0; 100 if(a.ms == NANEXP) 101 return sign >= 0; 102 if(a.ms == (NANEXP|NANSIGN)) 103 return sign <= 0; 104 return 0; 105 } 106