Return to sqrt.s CVS log

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

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

1.1     ! root        1: /*
        !             2:  * Copyright (c) 1987 Regents of the University of California.
        !             3:  * All rights reserved.
        !             4:  *
        !             5:  * Redistribution and use in source and binary forms are permitted
        !             6:  * provided that: (1) source distributions retain this entire copyright
        !             7:  * notice and comment, and (2) distributions including binaries display
        !             8:  * the following acknowledgement:  ``This product includes software
        !             9:  * developed by the University of California, Berkeley and its contributors''
        !            10:  * in the documentation or other materials provided with the distribution
        !            11:  * and in all advertising materials mentioning features or use of this
        !            12:  * software. Neither the name of the University nor the names of its
        !            13:  * contributors may be used to endorse or promote products derived
        !            14:  * from this software without specific prior written permission.
        !            15:  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
        !            16:  * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
        !            17:  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
        !            18:  *
        !            19:  * All recipients should regard themselves as participants in an ongoing
        !            20:  * research project and hence should feel obligated to report their
        !            21:  * experiences (good or bad) with these elementary function codes, using
        !            22:  * the sendbug(8) program, to the authors.
        !            23:  */
        !            24:        .data
        !            25:        .align  2
        !            26: _sccsid:
        !            27: .asciz "@(#)sqrt.s     5.5     (ucb.elefunt)   6/1/90"
        !            28: 
        !            29: /*
        !            30:  * double sqrt(arg)   revised August 15,1982
        !            31:  * double arg;
        !            32:  * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); }
        !            33:  * if arg is a reserved operand it is returned as it is
        !            34:  * W. Kahan's magic square root
        !            35:  * Coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82.
        !            36:  * Re-coded in tahoe assembly language by Z. Alex Liu 7/13/87.
        !            37:  *
        !            38:  * entry points:_d_sqrt                address of double arg is on the stack
        !            39:  *             _sqrt           double arg is on the stack
        !            40:  */
        !            41:        .text
        !            42:        .align  2
        !            43:        .globl  _sqrt
        !            44:        .globl  _d_sqrt
        !            45:        .globl  libm$dsqrt_r5
        !            46:        .set    EDOM,33
        !            47: 
        !            48: _d_sqrt:
        !            49:        .word   0x003c          # save r2-r5
        !            50:        movl    4(fp),r2
        !            51:        movl    (r2),r0
        !            52:        movl    4(r2),r1        # r0:r1 = x
        !            53:        brb     1f
        !            54: _sqrt:
        !            55:        .word   0x003c          # save r2-r5
        !            56:        movl    4(fp),r0
        !            57:        movl    8(fp),r1        # r0:r1 = x
        !            58: 1:     andl3   $0x7f800000,r0,r2       # r2 = biased exponent
        !            59:        bneq    2f
        !            60:        ret                     # biased exponent is zero -> 0 or reserved op.
        !            61: /*
        !            62:  *                             # internal procedure
        !            63:  *                             # ENTRY POINT FOR cdabs and cdsqrt
        !            64:  */
        !            65: libm$dsqrt_r5:                 # returns double square root scaled by 2^r6
        !            66:        .word   0x0000          # save nothing
        !            67: 2:     ldd     r0
        !            68:        std     r4
        !            69:        bleq    nonpos          # argument is not positive
        !            70:        andl3   $0xfffe0000,r4,r2
        !            71:        shar    $1,r2,r0
        !            72:        addl2   $0x203c0000,r0  # r0 has magic initial approximation
        !            73: /*
        !            74:  *                             # Do two steps of Heron's rule
        !            75:  *                             # ((arg/guess)+guess)/2 = better guess
        !            76:  */
        !            77:        ldf     r4
        !            78:        divf    r0
        !            79:        addf    r0
        !            80:        stf     r0
        !            81:        subl2   $0x800000,r0    # divide by two
        !            82:        ldf     r4
        !            83:        divf    r0
        !            84:        addf    r0
        !            85:        stf     r0
        !            86:        subl2   $0x800000,r0    # divide by two
        !            87: /*
        !            88:  *                             # Scale argument and approximation
        !            89:  *                             # to prevent over/underflow
        !            90:  */
        !            91:        andl3   $0x7f800000,r4,r1
        !            92:        subl2   $0x40800000,r1  # r1 contains scaling factor
        !            93:        subl2   r1,r4           # r4:r5 = n/s
        !            94:        movl    r0,r2
        !            95:        subl2   r1,r2           # r2 = a/s
        !            96: /*
        !            97:  *                             # Cubic step
        !            98:  *                             # b = a+2*a*(n-a*a)/(n+3*a*a) where
        !            99:  *                             # b is better approximation, a is approximation
        !           100:  *                             # and n is the original argument.
        !           101:  *                             # s := scale factor.
        !           102:  */
        !           103:        clrl    r1              # r0:r1 = a
        !           104:        clrl    r3              # r2:r3 = a/s
        !           105:        ldd     r0              # acc = a
        !           106:        muld    r2              # acc = a*a/s
        !           107:        std     r2              # r2:r3 = a*a/s
        !           108:        negd                    # acc = -a*a/s
        !           109:        addd    r4              # acc = n/s-a*a/s
        !           110:        std     r4              # r4:r5 = n/s-a*a/s
        !           111:        addl2   $0x1000000,r2   # r2:r3 = 4*a*a/s
        !           112:        ldd     r2              # acc = 4*a*a/s
        !           113:        addd    r4              # acc = n/s+3*a*a/s
        !           114:        std     r2              # r2:r3 = n/s+3*a*a/s
        !           115:        ldd     r0              # acc = a
        !           116:        muld    r4              # acc = a*n/s-a*a*a/s
        !           117:        divd    r2              # acc = a*(n-a*a)/(n+3*a*a)
        !           118:        std     r4              # r4:r5 = a*(n-a*a)/(n+3*a*a)
        !           119:        addl2   $0x800000,r4    # r4:r5 = 2*a*(n-a*a)/(n+3*a*a)
        !           120:        ldd     r4              # acc = 2*a*(n-a*a)/(n+3*a*a)
        !           121:        addd    r0              # acc = a+2*a*(n-a*a)/(n+3*a*a)
        !           122:        std     r0              # r0:r1 = a+2*a*(n-a*a)/(n+3*a*a)
        !           123:        ret                     # rsb
        !           124: nonpos:
        !           125:        bneq    negarg
        !           126:        ret                     # argument and root are zero
        !           127: negarg:
        !           128:        pushl   $EDOM
        !           129:        callf   $8,_infnan      # generate the reserved op fault
        !           130:        ret

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