Annotation of 42BSD/usr.lib/libnm/sqrt.s, revision 1.1.1.1

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

unix.superglobalmegacorp.com

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