Return to sqrt.s CVS log

Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / vax

Annotation of 43BSDReno/lib/libm/vax/sqrt.s, revision 1.1

1.1     ! root        1: # Copyright (c) 1985 Regents of the University of California.
        !             2: # All rights reserved.
        !             3: #
        !             4: # Redistribution and use in source and binary forms are permitted
        !             5: # provided that the above copyright notice and this paragraph are
        !             6: # duplicated in all such forms and that any documentation,
        !             7: # advertising materials, and other materials related to such
        !             8: # distribution and use acknowledge that the software was developed
        !             9: # by the University of California, Berkeley.  The name of the
        !            10: # University may not be used to endorse or promote products derived
        !            11: # from this software without specific prior written permission.
        !            12: # THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
        !            13: # IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
        !            14: # WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
        !            15: #
        !            16: # All recipients should regard themselves as participants in an ongoing
        !            17: # research project and hence should feel obligated to report their
        !            18: # experiences (good or bad) with these elementary function codes, using
        !            19: # the sendbug(8) program, to the authors.
        !            20: #
        !            21: #      @(#)sqrt.s      5.3 (Berkeley) 6/30/88
        !            22: #
        !            23:        .data
        !            24:        .align  2
        !            25: _sccsid:
        !            26: .asciz "@(#)sqrt.s     1.1 (Berkeley) 8/21/85; 5.3 (ucb.elefunt) 6/30/88"
        !            27: 
        !            28: /*
        !            29:  * double sqrt(arg)   revised August 15,1982
        !            30:  * double arg;
        !            31:  * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); }
        !            32:  * if arg is a reserved operand it is returned as it is
        !            33:  * W. Kahan's magic square root
        !            34:  * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82
        !            35:  *
        !            36:  * entry points:_d_sqrt                address of double arg is on the stack
        !            37:  *             _sqrt           double arg is on the stack
        !            38:  */
        !            39:        .text
        !            40:        .align  1
        !            41:        .globl  _sqrt
        !            42:        .globl  _d_sqrt
        !            43:        .globl  libm$dsqrt_r5
        !            44:        .set    EDOM,33
        !            45: 
        !            46: _d_sqrt:
        !            47:        .word   0x003c          # save r5,r4,r3,r2
        !            48:        movq    *4(ap),r0
        !            49:        jmp     dsqrt2
        !            50: _sqrt:
        !            51:        .word   0x003c          # save r5,r4,r3,r2
        !            52:        movq    4(ap),r0
        !            53: dsqrt2:        bicw3   $0x807f,r0,r2   # check exponent of input
        !            54:        jeql    noexp           # biased exponent is zero -> 0.0 or reserved
        !            55:        bsbb    libm$dsqrt_r5
        !            56: noexp: ret
        !            57: 
        !            58: /* **************************** internal procedure */
        !            59: 
        !            60: libm$dsqrt_r5:                 # ENTRY POINT FOR cdabs and cdsqrt
        !            61:                                # returns double square root scaled by
        !            62:                                # 2^r6
        !            63: 
        !            64:        movd    r0,r4
        !            65:        jleq    nonpos          # argument is not positive
        !            66:        movzwl  r4,r2
        !            67:        ashl    $-1,r2,r0
        !            68:        addw2   $0x203c,r0      # r0 has magic initial approximation
        !            69: /*
        !            70:  * Do two steps of Heron's rule
        !            71:  * ((arg/guess) + guess) / 2 = better guess
        !            72:  */
        !            73:        divf3   r0,r4,r2
        !            74:        addf2   r2,r0
        !            75:        subw2   $0x80,r0        # divide by two
        !            76: 
        !            77:        divf3   r0,r4,r2
        !            78:        addf2   r2,r0
        !            79:        subw2   $0x80,r0        # divide by two
        !            80: 
        !            81: /* Scale argument and approximation to prevent over/underflow */
        !            82: 
        !            83:        bicw3   $0x807f,r4,r1
        !            84:        subw2   $0x4080,r1              # r1 contains scaling factor
        !            85:        subw2   r1,r4
        !            86:        movl    r0,r2
        !            87:        subw2   r1,r2
        !            88: 
        !            89: /* Cubic step
        !            90:  *
        !            91:  * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation,
        !            92:  * a is approximation, and n is the original argument.
        !            93:  * (let s be scale factor in the following comments)
        !            94:  */
        !            95:        clrl    r1
        !            96:        clrl    r3
        !            97:        muld2   r0,r2                   # r2:r3 = a*a/s
        !            98:        subd2   r2,r4                   # r4:r5 = n/s - a*a/s
        !            99:        addw2   $0x100,r2               # r2:r3 = 4*a*a/s
        !           100:        addd2   r4,r2                   # r2:r3 = n/s + 3*a*a/s
        !           101:        muld2   r0,r4                   # r4:r5 = a*n/s - a*a*a/s
        !           102:        divd2   r2,r4                   # r4:r5 = a*(n-a*a)/(n+3*a*a)
        !           103:        addw2   $0x80,r4                # r4:r5 = 2*a*(n-a*a)/(n+3*a*a)
        !           104:        addd2   r4,r0                   # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a)
        !           105:        rsb                             # DONE!
        !           106: nonpos:
        !           107:        jneq    negarg
        !           108:        ret                     # argument and root are zero
        !           109: negarg:
        !           110:        pushl   $EDOM
        !           111:        calls   $1,_infnan      # generate the reserved op fault
        !           112:        ret

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