home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 4 / FreshFish_May-June1994.bin / bbs / gnu / libm-src.lha / src / amiga / libm / tahoe / cabs.s next >
Encoding:
Text File  |  1993-09-23  |  3.9 KB  |  119 lines
  1. # Copyright (c) 1987 Regents of the University of California.
  2. # All rights reserved.
  3. #
  4. # Redistribution and use in source and binary forms, with or without
  5. # modification, are permitted provided that the following conditions
  6. # are met:
  7. # 1. Redistributions of source code must retain the above copyright
  8. #    notice, this list of conditions and the following disclaimer.
  9. # 2. Redistributions in binary form must reproduce the above copyright
  10. #    notice, this list of conditions and the following disclaimer in the
  11. #    documentation and/or other materials provided with the distribution.
  12. # 3. All advertising materials mentioning features or use of this software
  13. #    must display the following acknowledgement:
  14. #    This product includes software developed by the University of
  15. #    California, Berkeley and its contributors.
  16. # 4. Neither the name of the University nor the names of its contributors
  17. #    may be used to endorse or promote products derived from this software
  18. #    without specific prior written permission.
  19. #
  20. # THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  21. # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  22. # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  23. # ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  24. # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  25. # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  26. # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  27. # HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  28. # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  29. # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  30. # SUCH DAMAGE.
  31. #
  32. #    @(#)cabs.s    5.5 (Berkeley) 10/9/90
  33. #
  34.     .data
  35.     .align    2
  36. _sccsid:
  37. .asciz    "@(#)cabs.s    5.5    5.5    (ucb.elefunt)    10/9/90"
  38.  
  39. # double precision complex absolute value
  40. # CABS by W. Kahan, 9/7/80.
  41. # Revised for reserved operands by E. LeBlanc, 8/18/82
  42. # argument for complex absolute value by reference, *4(fp)
  43. # argument for cabs and hypot (C fcns) by value, 4(fp)
  44. # output is in r0:r1
  45.  
  46.     .text
  47.     .align    2
  48.     .globl  _cabs
  49.     .globl  _hypot
  50.     .globl    _z_abs
  51.  
  52. #    entry for c functions cabs and hypot
  53. _cabs:
  54. _hypot:
  55.     .word    0x807c        # save r2-r6, enable floating overflow
  56.     movl    16(fp),r3
  57.     movl    12(fp),r2    # r2:3 = y
  58.     movl    8(fp),r1
  59.     movl    4(fp),r0    # r0:1 = x
  60.     brb    1f
  61. #    entry for Fortran use, call by:   d = abs(z)
  62. _z_abs:
  63.     .word    0x807c        # save r2-r6, enable floating overflow
  64.     movl    4(fp),r4    # indirect addressing is necessary here
  65.     movl    12(r4),r3    #
  66.     movl    8(r4),r2    # r2:3 = y
  67.     movl    4(r4),r1    #
  68.     movl    (r4),r0        # r0:1 = x
  69. 1:    andl3    $0xff800000,r0,r4    # r4 has signed biased exp of x
  70.     cmpl    $0x80000000,r4
  71.     beql    2f        # x is a reserved operand, so return it
  72.     andl3    $0xff800000,r2,r5    # r5 has signed biased exp of y
  73.     cmpl    $0x80000000,r5
  74.     bneq    3f        # y isn't a reserved operand
  75.     movl    r3,r1
  76.     movl    r2,r0        # return y if it's reserved
  77. 2:    ret
  78.  
  79. 3:    callf    $4,regs_set    # r0:1 = dsqrt(x^2+y^2)/2^r6
  80.     addl2    r6,r0        # unscaled cdabs in r0:1
  81.     jvc    2b        # unless it overflows
  82.     subl2    $0x800000,r0    # halve r0 to get meaningful overflow
  83.     ldd    r0
  84.     addd    r0        # overflow; r0 is half of true abs value
  85.     ret
  86.  
  87. regs_set:
  88.     .word    0x0000
  89.     andl2    $0x7fffffff,r0    # r0:r1 = dabs(x)
  90.     andl2    $0x7fffffff,r2    # r2:r3 = dabs(y)
  91.     cmpl    r0,r2
  92.     bgeq    4f
  93.     movl    r1,r5
  94.     movl    r0,r4
  95.     movl    r3,r1
  96.     movl    r2,r0
  97.     movl    r5,r3
  98.     movl    r4,r2        # force y's exp <= x's exp
  99. 4:    andl3    $0xff800000,r0,r6    # r6 = exponent(x) + bias(129)
  100.     beql    5f        # if x = y = 0 then cdabs(x,y) = 0
  101.     subl2    $0x47800000,r6    # r6 = exponent(x) - 14
  102.     subl2    r6,r0        # 2^14 <= scaled x < 2^15
  103.     bitl    $0xff800000,r2
  104.     beql    5f        # if y = 0 return dabs(x)
  105.     subl2    r6,r2
  106.     cmpl    $0x37800000,r2    # if scaled y < 2^-18
  107.     bgtr    5f        #   return dabs(x)
  108.     ldd    r0
  109.     muld    r0
  110.     std    r0        # r0:1 = scaled x^2
  111.     ldd    r2
  112.     muld    r2        # acc = scaled y^2
  113.     addd    r0
  114.     std    r0
  115.     pushl    r1
  116.     pushl    r0
  117.     callf    $12,_sqrt    # r0:1 = dsqrt(x^2+y^2)/2^r6
  118. 5:    ret
  119.