tsqrt.c - plan9port - [fork] Plan 9 from user space
 (HTM) git clone git://src.adamsgaard.dk/plan9port
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
       tsqrt.c (675B)
       ---
            1 /*
            2         sqrt returns the square root of its floating
            3         point argument. Newton's method.
            4 
            5         calls frexp
            6 */
            7 
            8 #include <u.h>
            9 #include <libc.h>
           10 
           11 double
           12 sqrt(double arg)
           13 {
           14         double x, temp;
           15         int exp, i;
           16 
           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          * NOTE
           29          * this wont work on 1's comp
           30          */
           31         if(exp & 1) {
           32                 x *= 2;
           33                 exp--;
           34         }
           35         temp = 0.5 * (1.0+x);
           36 
           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         else
           48                 temp /= 1L << (-exp/2);
           49         for(i=0; i<=4; i++)
           50                 temp = 0.5*(temp + arg/temp);
           51         return temp;
           52 }