Return to sqrt.s CVS log

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

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

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