4.3BSD - /usr/src/usr.lib/libm/VAX/cabs.s

   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

Defined functions