Return to sqrt.c CVS log

Up to [Research Unix] / researchv9 / libc / math

Annotation of researchv9/libc/math/sqrt.c, revision 1.1.1.1

1.1       root        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: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.