Return to sqrt.s CVS log

Up to [CSRG BSD Unix] / 43BSDReno / old / libm / liboldnm

Annotation of 43BSDReno/old/libm/liboldnm/sqrt.s, revision 1.1.1.1

1.1       root        1: #
                      2: # Copyright (c) 1980 Regents of the University of California.
                      3: # All rights reserved.  The Berkeley software License Agreement
                      4: # specifies the terms and conditions for redistribution.
                      5: #
                      6: #      @(#)sqrt.s      5.1 (Berkeley) 5/8/85
                      7: # 
                      8: #
                      9: # double sqrt(arg):revised July 18,1980
                     10: # double arg
                     11: # if(arg<0.0) { _errno=EDOM; return(-sqrt(-arg)) }
                     12: # W. Kahan's magic sqrt
                     13: # coded by Heidi Stettner
                     14: 
                     15:        .set    EDOM,98
                     16:        .text
                     17:        .align  1
                     18:        .globl  _sqrt
                     19:        .globl  dsqrt_r5
                     20:        .globl  _errno
                     21: _sqrt:
                     22:        .word   0x003c          # save r5,r4,r3,r2
                     23:        bispsw  $0xe0
                     24:        movd    4(ap),r0
                     25:        bsbb    dsqrt_r5
                     26:        ret
                     27: dsqrt_r5:
                     28:        movd    r0,r4
                     29:        jleq    nonpos          #argument is not positive
                     30:        movzwl  r4,r2
                     31:        ashl    $-1,r2,r0
                     32:        addw2   $0x203c,r0      #r0 has magic initial appx
                     33: 
                     34: # Do two steps of Heron's rule
                     35: 
                     36:        divf3   r0,r4,r2        
                     37:        addf2   r2,r0
                     38:        subw2   $0x80,r0
                     39: 
                     40:        divf3   r0,r4,r2
                     41:        addf2   r2,r0
                     42:        subw2   $0x80,r0
                     43: 
                     44: 
                     45: # Scale argument and approximation to prevent over/underflow
                     46: # NOTE: The following four steps would not be necessary if underflow
                     47: #       were gentle.
                     48: 
                     49:        bicw3   $0xffff807f,r4,r1
                     50:        subw2   $0x4080,r1              # r1 contains scaling factor
                     51:        subw2   r1,r4
                     52:        movl    r0,r2
                     53:        subw2   r1,r2
                     54: 
                     55: # Cubic step
                     56: 
                     57:        clrl    r1
                     58:        clrl    r3
                     59:        muld2   r0,r2
                     60:        subd2   r2,r4
                     61:        addw2   $0x100,r2
                     62:        addd2   r4,r2
                     63:        muld2   r0,r4
                     64:        divd2   r2,r4
                     65:        addw2   $0x80,r4
                     66:        addd2   r4,r0
                     67:        rsb
                     68: nonpos:
                     69:        jneq    negarg
                     70:        clrd    r0              #argument is zero
                     71:        rsb
                     72: negarg:
                     73:        movl    $EDOM,_errno
                     74:        mnegd   r4,-(sp)
                     75:        calls   $2,_sqrt
                     76:        mnegd   r0,r0           # returns -sqrt(-arg)
                     77:        ret

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