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

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.