Return to sqrt.c CVS log

Up to [CSRG BSD Unix] / 43BSD / usr.lib / libom

Annotation of 43BSD/usr.lib/libom/sqrt.c, revision 1.1.1.1

1.1       root        1: /*     @(#)sqrt.c      4.1     12/25/82        */
                      2: 
                      3: /*
                      4:        sqrt returns the square root of its floating
                      5:        point argument. Newton's method.
                      6: 
                      7:        calls frexp
                      8: */
                      9: 
                     10: #include <errno.h>
                     11: 
                     12: int errno;
                     13: double frexp();
                     14: 
                     15: double
                     16: sqrt(arg)
                     17: double arg;
                     18: {
                     19:        double x, temp;
                     20:        int exp;
                     21:        int i;
                     22: 
                     23:        if(arg <= 0.) {
                     24:                if(arg < 0.)
                     25:                        errno = EDOM;
                     26:                return(0.);
                     27:        }
                     28:        x = frexp(arg,&exp);
                     29:        while(x < 0.5) {
                     30:                x *= 2;
                     31:                exp--;
                     32:        }
                     33:        /*
                     34:         * NOTE
                     35:         * this wont work on 1's comp
                     36:         */
                     37:        if(exp & 1) {
                     38:                x *= 2;
                     39:                exp--;
                     40:        }
                     41:        temp = 0.5*(1.0+x);
                     42: 
                     43:        while(exp > 60) {
                     44:                temp *= (1L<<30);
                     45:                exp -= 60;
                     46:        }
                     47:        while(exp < -60) {
                     48:                temp /= (1L<<30);
                     49:                exp += 60;
                     50:        }
                     51:        if(exp >= 0)
                     52:                temp *= 1L << (exp/2);
                     53:        else
                     54:                temp /= 1L << (-exp/2);
                     55:        for(i=0; i<=4; i++)
                     56:                temp = 0.5*(temp + arg/temp);
                     57:        return(temp);
                     58: }

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