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