Return to sqrt.s CVS log

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

Annotation of 43BSDTahoe/usr.lib/libm/vax/sqrt.s, revision 1.1.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.