2 28994509 2004-04-21 devnull sqrt returns the square root of its floating
3 28994509 2004-04-21 devnull point argument. Newton's method.
5 28994509 2004-04-21 devnull calls frexp
8 28994509 2004-04-21 devnull #include <u.h>
9 28994509 2004-04-21 devnull #include <libc.h>
12 28994509 2004-04-21 devnull sqrt(double arg)
14 28994509 2004-04-21 devnull double x, temp;
15 28994509 2004-04-21 devnull int exp, i;
17 28994509 2004-04-21 devnull if(arg <= 0) {
18 28994509 2004-04-21 devnull if(arg < 0)
19 28994509 2004-04-21 devnull return 0.;
20 28994509 2004-04-21 devnull return 0;
22 28994509 2004-04-21 devnull x = frexp(arg, &exp);
23 28994509 2004-04-21 devnull while(x < 0.5) {
29 28994509 2004-04-21 devnull * this wont work on 1's comp
31 28994509 2004-04-21 devnull if(exp & 1) {
35 28994509 2004-04-21 devnull temp = 0.5 * (1.0+x);
37 28994509 2004-04-21 devnull while(exp > 60) {
38 28994509 2004-04-21 devnull temp *= (1L<<30);
39 28994509 2004-04-21 devnull exp -= 60;
41 28994509 2004-04-21 devnull while(exp < -60) {
42 28994509 2004-04-21 devnull temp /= (1L<<30);
43 28994509 2004-04-21 devnull exp += 60;
45 28994509 2004-04-21 devnull if(exp >= 0)
46 28994509 2004-04-21 devnull temp *= 1L << (exp/2);
48 28994509 2004-04-21 devnull temp /= 1L << (-exp/2);
49 28994509 2004-04-21 devnull for(i=0; i<=4; i++)
50 28994509 2004-04-21 devnull temp = 0.5*(temp + arg/temp);
51 28994509 2004-04-21 devnull return temp;