Return to sqrtd.c CVS log

Up to [CSRG BSD Unix] / 43BSDTahoe / cci / libM

Annotation of 43BSDTahoe/cci/libM/sqrtd.c, revision 1.1

1.1     ! root        1: /*
        !             2:  * Square Root.
        !             3:  * New version by Les Powers (3/26/85).
        !             4:  * If the argument is less than zero, an error is reported.
        !             5:  * The single precision square root instruction "sqrtf" is
        !             6:  * used to generate the initial approximation.  Then two
        !             7:  * iterations of the Newton-Raphson method are used to improve
        !             8:  * the result to 56 bits.  Each iteration of the Newton-Raphson
        !             9:  * method uses the following replacement:
        !            10:  *     B = (B + A/B)/2
        !            11:  * where A is the argument and B is the approximation.
        !            12:  * In the following assembly code, the division by two is
        !            13:  * accomplished by subtracting one from the exponent field.
        !            14:  * To subtract one from the exponent field, the hex integer value
        !            15:  * 0x00800000 is subtracted from the first longword of the double
        !            16:  * precision number in registers r0 and r1.  The hex value 0x00800000
        !            17:  * is equal to decimal value 8838608.
        !            18:  */
        !            19: 
        !            20: #include <errno.h>
        !            21: 
        !            22: int    errno;
        !            23: 
        !            24: double
        !            25: sqrt(arg)
        !            26: double arg;
        !            27: {
        !            28:        if (arg > 0.0) {
        !            29:                asm("   ldf     4(fp)")
        !            30:                asm("   sqrtf")
        !            31:                asm("   clrl    r1")
        !            32:                asm("   stf     r0")
        !            33:                asm("   ldd     4(fp)")
        !            34:                asm("   divd    r0")
        !            35:                asm("   addd    r0")
        !            36:                asm("   std     r0")
        !            37:                asm("   subl2   $8388608,r0")
        !            38:                asm("   ldd     4(fp)")
        !            39:                asm("   divd    r0")
        !            40:                asm("   addd    r0")
        !            41:                asm("   std     r0")
        !            42:                asm("   subl2   $8388608,r0")
        !            43:        } else {
        !            44:                if (arg < 0.0)
        !            45:                        errno = EDOM;
        !            46:                return (0.0);
        !            47:        }
        !            48: }

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