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     ! 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.