Return to sqrt.c CVS log

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

Annotation of researchv10dc/libc/math/sqrt.c, revision 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.