Annotation of 43BSD/usr.lib/libm/VAX/cabs.s, revision 1.1.1.1

1.1       root        1: # 
                      2: # Copyright (c) 1985 Regents of the University of California.
                      3: # 
                      4: # Use and reproduction of this software are granted  in  accordance  with
                      5: # the terms and conditions specified in  the  Berkeley  Software  License
                      6: # Agreement (in particular, this entails acknowledgement of the programs'
                      7: # source, and inclusion of this notice) with the additional understanding
                      8: # that  all  recipients  should regard themselves as participants  in  an
                      9: # ongoing  research  project and hence should  feel  obligated  to report
                     10: # their  experiences (good or bad) with these elementary function  codes,
                     11: # using "sendbug 4bsd-bugs@BERKELEY", to the authors.
                     12: #
                     13: 
                     14: # @(#)cabs.s   1.2 (Berkeley) 8/21/85
                     15: 
                     16: # double precision complex absolute value
                     17: # CABS by W. Kahan, 9/7/80.
                     18: # Revised for reserved operands by E. LeBlanc, 8/18/82
                     19: # argument for complex absolute value by reference, *4(ap)
                     20: # argument for cabs and hypot (C fcns) by value, 4(ap)
                     21: # output is in r0:r1 (error less than 0.86 ulps)
                     22: 
                     23:        .text
                     24:        .align  1
                     25:        .globl  _cabs
                     26:        .globl  _hypot
                     27:        .globl  _z_abs
                     28:        .globl  libm$cdabs_r6
                     29:        .globl  libm$dsqrt_r5
                     30: 
                     31: #      entry for c functions cabs and hypot
                     32: _cabs:
                     33: _hypot:
                     34:        .word   0x807c          # save r2-r6, enable floating overflow
                     35:        movq    4(ap),r0        # r0:1 = x
                     36:        movq    12(ap),r2       # r2:3 = y
                     37:        jmp     cabs2
                     38: #      entry for Fortran use, call by:   d = abs(z)
                     39: _z_abs:
                     40:        .word   0x807c          # save r2-r6, enable floating overflow
                     41:        movl    4(ap),r2        # indirect addressing is necessary here
                     42:        movq    (r2)+,r0        # r0:1 = x
                     43:        movq    (r2),r2         # r2:3 = y
                     44: 
                     45: cabs2:
                     46:        bicw3   $0x7f,r0,r4     # r4 has signed biased exp of x
                     47:        cmpw    $0x8000,r4
                     48:        jeql    return          # x is a reserved operand, so return it
                     49:        bicw3   $0x7f,r2,r5     # r5 has signed biased exp of y
                     50:        cmpw    $0x8000,r5
                     51:        jneq    cont            # y isn't a reserved operand
                     52:        movq    r2,r0           # return y if it's reserved
                     53:        ret
                     54: 
                     55: cont:
                     56:        bsbb    regs_set        # r0:1 = dsqrt(x^2+y^2)/2^r6
                     57:        addw2   r6,r0           # unscaled cdabs in r0:1
                     58:        jvc     return          # unless it overflows
                     59:        subw2   $0x80,r0        # halve r0 to get meaningful overflow
                     60:        addd2   r0,r0           # overflow; r0 is half of true abs value
                     61: return:
                     62:        ret
                     63: 
                     64: libm$cdabs_r6:                 # ENTRY POINT for cdsqrt
                     65:                                # calculates a scaled (factor in r6)
                     66:                                # complex absolute value
                     67: 
                     68:        movq    (r4)+,r0        # r0:r1 = x via indirect addressing
                     69:        movq    (r4),r2         # r2:r3 = y via indirect addressing
                     70: 
                     71:        bicw3   $0x7f,r0,r5     # r5 has signed biased exp of x
                     72:        cmpw    $0x8000,r5
                     73:        jeql    cdreserved      # x is a reserved operand
                     74:        bicw3   $0x7f,r2,r5     # r5 has signed biased exp of y
                     75:        cmpw    $0x8000,r5
                     76:        jneq    regs_set        # y isn't a reserved operand either?
                     77: 
                     78: cdreserved:
                     79:        movl    *4(ap),r4       # r4 -> (u,v), if x or y is reserved
                     80:        movq    r0,(r4)+        # copy u and v as is and return
                     81:        movq    r2,(r4)         # (again addressing is indirect)
                     82:        ret
                     83: 
                     84: regs_set:
                     85:        bicw2   $0x8000,r0      # r0:r1 = dabs(x)
                     86:        bicw2   $0x8000,r2      # r2:r3 = dabs(y)
                     87:        cmpw    r0,r2
                     88:        jgeq    ordered
                     89:        movq    r0,r4
                     90:        movq    r2,r0
                     91:        movq    r4,r2           # force y's exp <= x's exp
                     92: ordered:
                     93:        bicw3   $0x7f,r0,r6     # r6 = exponent(x) + bias(129)
                     94:        jeql    retsb           # if x = y = 0 then cdabs(x,y) = 0
                     95:        subw2   $0x4780,r6      # r6 = exponent(x) - 14
                     96:        subw2   r6,r0           # 2^14 <= scaled x < 2^15
                     97:        bitw    $0xff80,r2
                     98:        jeql    retsb           # if y = 0 return dabs(x)
                     99:        subw2   r6,r2
                    100:        cmpw    $0x3780,r2      # if scaled y < 2^-18
                    101:        jgtr    retsb           #   return dabs(x)
                    102:        emodd   r0,$0,r0,r4,r0  # r4 + r0:1 = scaled x^2
                    103:        emodd   r2,$0,r2,r5,r2  # r5 + r2:3 = scaled y^2
                    104:        addd2   r2,r0
                    105:        addl2   r5,r4
                    106:        cvtld   r4,r2
                    107:        addd2   r2,r0           # r0:1 = scaled x^2 + y^2
                    108:        jmp     libm$dsqrt_r5   # r0:1 = dsqrt(x^2+y^2)/2^r6
                    109: retsb:
                    110:        rsb                     # error < 0.86 ulp

unix.superglobalmegacorp.com

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