Blob
1 /*2 sqrt returns the square root of its floating3 point argument. Newton's method.5 calls frexp6 */8 #include <u.h>9 #include <libc.h>11 double12 sqrt(double arg)13 {14 double x, temp;15 int exp, i;17 if(arg <= 0) {18 if(arg < 0)19 return 0.;20 return 0;21 }22 x = frexp(arg, &exp);23 while(x < 0.5) {24 x *= 2;25 exp--;26 }27 /*28 * NOTE29 * this wont work on 1's comp30 */31 if(exp & 1) {32 x *= 2;33 exp--;34 }35 temp = 0.5 * (1.0+x);37 while(exp > 60) {38 temp *= (1L<<30);39 exp -= 60;40 }41 while(exp < -60) {42 temp /= (1L<<30);43 exp += 60;44 }45 if(exp >= 0)46 temp *= 1L << (exp/2);47 else48 temp /= 1L << (-exp/2);49 for(i=0; i<=4; i++)50 temp = 0.5*(temp + arg/temp);51 return temp;52 }