Return to sqrt.s CVS log

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

Annotation of 43BSD/usr.lib/libm/VAX/sqrt.s, revision 1.1.1.1

1.1       root        1: /* 
                      2:  * Copyright (c) 1985 Regents of the University of California.
                      3:  * 
                      4:  * Use and reproduction of this software are granted  in  accordance  with
                      5:  * the terms and conditions specified in  the  Berkeley  Software  License
                      6:  * Agreement (in particular, this entails acknowledgement of the programs'
                      7:  * source, and inclusion of this notice) with the additional understanding
                      8:  * that  all  recipients  should regard themselves as participants  in  an
                      9:  * ongoing  research  project and hence should  feel  obligated  to report
                     10:  * their  experiences (good or bad) with these elementary function  codes,
                     11:  * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
                     12:  *
                     13:  *
                     14:  * @(#)sqrt.s  1.1 (Berkeley) 8/21/85
                     15:  *
                     16:  * double sqrt(arg)   revised August 15,1982
                     17:  * double arg;
                     18:  * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); }
                     19:  * if arg is a reserved operand it is returned as it is
                     20:  * W. Kahan's magic square root
                     21:  * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82
                     22:  *
                     23:  * entry points:_d_sqrt                address of double arg is on the stack
                     24:  *             _sqrt           double arg is on the stack
                     25:  */
                     26:        .text
                     27:        .align  1
                     28:        .globl  _sqrt
                     29:        .globl  _d_sqrt
                     30:        .globl  libm$dsqrt_r5
                     31:        .set    EDOM,33
                     32: 
                     33: _d_sqrt:
                     34:        .word   0x003c          # save r5,r4,r3,r2
                     35:        movq    *4(ap),r0
                     36:        jmp     dsqrt2
                     37: _sqrt:
                     38:        .word   0x003c          # save r5,r4,r3,r2
                     39:        movq    4(ap),r0
                     40: dsqrt2:        bicw3   $0x807f,r0,r2   # check exponent of input
                     41:        jeql    noexp           # biased exponent is zero -> 0.0 or reserved
                     42:        bsbb    libm$dsqrt_r5
                     43: noexp: ret
                     44: 
                     45: /* **************************** internal procedure */
                     46: 
                     47: libm$dsqrt_r5:                 # ENTRY POINT FOR cdabs and cdsqrt
                     48:                                # returns double square root scaled by
                     49:                                # 2^r6
                     50: 
                     51:        movd    r0,r4
                     52:        jleq    nonpos          # argument is not positive
                     53:        movzwl  r4,r2
                     54:        ashl    $-1,r2,r0
                     55:        addw2   $0x203c,r0      # r0 has magic initial approximation
                     56: /*
                     57:  * Do two steps of Heron's rule
                     58:  * ((arg/guess) + guess) / 2 = better guess
                     59:  */
                     60:        divf3   r0,r4,r2
                     61:        addf2   r2,r0
                     62:        subw2   $0x80,r0        # divide by two
                     63: 
                     64:        divf3   r0,r4,r2
                     65:        addf2   r2,r0
                     66:        subw2   $0x80,r0        # divide by two
                     67: 
                     68: /* Scale argument and approximation to prevent over/underflow */
                     69: 
                     70:        bicw3   $0x807f,r4,r1
                     71:        subw2   $0x4080,r1              # r1 contains scaling factor
                     72:        subw2   r1,r4
                     73:        movl    r0,r2
                     74:        subw2   r1,r2
                     75: 
                     76: /* Cubic step
                     77:  *
                     78:  * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation,
                     79:  * a is approximation, and n is the original argument.
                     80:  * (let s be scale factor in the following comments)
                     81:  */
                     82:        clrl    r1
                     83:        clrl    r3
                     84:        muld2   r0,r2                   # r2:r3 = a*a/s
                     85:        subd2   r2,r4                   # r4:r5 = n/s - a*a/s
                     86:        addw2   $0x100,r2               # r2:r3 = 4*a*a/s
                     87:        addd2   r4,r2                   # r2:r3 = n/s + 3*a*a/s
                     88:        muld2   r0,r4                   # r4:r5 = a*n/s - a*a*a/s
                     89:        divd2   r2,r4                   # r4:r5 = a*(n-a*a)/(n+3*a*a)
                     90:        addw2   $0x80,r4                # r4:r5 = 2*a*(n-a*a)/(n+3*a*a)
                     91:        addd2   r4,r0                   # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a)
                     92:        rsb                             # DONE!
                     93: nonpos:
                     94:        jneq    negarg
                     95:        ret                     # argument and root are zero
                     96: negarg:
                     97:        pushl   $EDOM
                     98:        calls   $1,_infnan      # generate the reserved op fault
                     99:        ret

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