lib/libm/arch/vax/n_cabs.S

   1 /*      $NetBSD: n_cabs.S,v 1.5 2003/08/07 16:44:45 agc Exp $   */
   2 /*
   3  * Copyright (c) 1985, 1993
   4  *      The Regents of the University of California.  All rights reserved.
   5  *
   6  * Redistribution and use in source and binary forms, with or without
   7  * modification, are permitted provided that the following conditions
   8  * are met:
   9  * 1. Redistributions of source code must retain the above copyright
  10  *    notice, this list of conditions and the following disclaimer.
  11  * 2. Redistributions in binary form must reproduce the above copyright
  12  *    notice, this list of conditions and the following disclaimer in the
  13  *    documentation and/or other materials provided with the distribution.
  14  * 3. Neither the name of the University nor the names of its contributors
  15  *    may be used to endorse or promote products derived from this software
  16  *    without specific prior written permission.
  17  *
  18  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  19  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  22  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  23  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  24  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  25  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  26  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  27  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  28  * SUCH DAMAGE.
  29  *
  30  *      @(#)cabs.s      8.1 (Berkeley) 6/4/93
  31  */
  32
  33 #include <machine/asm.h>
  34
  35         .globl  _C_LABEL(__libm_dsqrt_r5)
  36 /*
  37  * double precision complex absolute value
  38  * CABS by W. Kahan, 9/7/80.
  39  * Revised for reserved operands by E. LeBlanc, 8/18/82
  40  * argument for complex absolute value by reference, *4(%ap)
  41  * argument for cabs and hypot (C fcns) by value, 4(%ap)
  42  * output is in %r0:%r1 (error less than 0.86 ulps)
  43  */
  44
  45 /*      entry for c functions cabs and hypot */
  46 #ifdef WEAK_ALIAS
  47 WEAK_ALIAS(hypotf, _hypotf)
  48 #endif
  49
  50 ENTRY(_hypotf, 0)
  51         cvtfd   4(%ap),-(%sp)
  52         calls   $2,_C_LABEL(_hypot)
  53         cvtdf   %r0,%r0
  54         ret
  55
  56 #ifdef WEAK_ALIAS
  57 WEAK_ALIAS(hypot, _hypot)
  58 #endif
  59
  60 ALTENTRY(cabs)
  61 ENTRY(_hypot, 0x8040)           # save %r6, enable floating overflow
  62         movq    4(%ap),%r0      # %r0:1 = x
  63         movq    12(%ap),%r2     # %r2:3 = y
  64         jbr     cabs2
  65
  66 /*      entry for Fortran use, call by:   d = abs(z) */
  67 ENTRY(z_abs, 0x8040)            # save %r6, enable floating overflow
  68         movl    4(%ap),%r2      # indirect addressing is necessary here
  69         movq    (%r2)+,%r0      # %r0:1 = x
  70         movq    (%r2),%r2               # %r2:3 = y
  71
  72 cabs2:
  73         bicw3   $0x7f,%r0,%r4   # %r4 has signed biased exp of x
  74         cmpw    $0x8000,%r4
  75         jeql    return          # x is a reserved operand, so return it
  76         bicw3   $0x7f,%r2,%r5   # %r5 has signed biased exp of y
  77         cmpw    $0x8000,%r5
  78         jneq    cont            /* y isn't a reserved operand */
  79         movq    %r2,%r0         /* return y if it's reserved */
  80         ret
  81
  82 cont:
  83         bsbb    regs_set        # %r0:1 = dsqrt(x^2+y^2)/2^%r6
  84         addw2   %r6,%r0         # unscaled cdabs in %r0:1
  85         jvc     return          # unless it overflows
  86         subw2   $0x80,%r0       # halve %r0 to get meaningful overflow
  87         addd2   %r0,%r0         # overflow; %r0 is half of true abs value
  88 return:
  89         ret
  90
  91 ENTRY(__libm_cdabs_r6,0)        # ENTRY POINT for cdsqrt
  92                                 # calculates a scaled (factor in %r6)
  93                                 # complex absolute value
  94
  95         movq    (%r4)+,%r0      # %r0:%r1 = x via indirect addressing
  96         movq    (%r4),%r2               # %r2:%r3 = y via indirect addressing
  97
  98         bicw3   $0x7f,%r0,%r5   # %r5 has signed biased exp of x
  99         cmpw    $0x8000,%r5
 100         jeql    cdreserved      # x is a reserved operand
 101         bicw3   $0x7f,%r2,%r5   # %r5 has signed biased exp of y
 102         cmpw    $0x8000,%r5
 103         jneq    regs_set        /* y isn't a reserved operand either? */
 104
 105 cdreserved:
 106         movl    *4(%ap),%r4     # %r4 -> (u,v), if x or y is reserved
 107         movq    %r0,(%r4)+      # copy u and v as is and return
 108         movq    %r2,(%r4)               # (again addressing is indirect)
 109         ret
 110
 111 regs_set:
 112         bicw2   $0x8000,%r0     # %r0:%r1 = dabs(x)
 113         bicw2   $0x8000,%r2     # %r2:%r3 = dabs(y)
 114         cmpw    %r0,%r2
 115         jgeq    ordered
 116         movq    %r0,%r4
 117         movq    %r2,%r0
 118         movq    %r4,%r2         # force y's exp <= x's exp
 119 ordered:
 120         bicw3   $0x7f,%r0,%r6   # %r6 = exponent(x) + bias(129)
 121         jeql    retsb           # if x = y = 0 then cdabs(x,y) = 0
 122         subw2   $0x4780,%r6     # %r6 = exponent(x) - 14
 123         subw2   %r6,%r0         # 2^14 <= scaled x < 2^15
 124         bitw    $0xff80,%r2
 125         jeql    retsb           # if y = 0 return dabs(x)
 126         subw2   %r6,%r2
 127         cmpw    $0x3780,%r2     # if scaled y < 2^-18
 128         jgtr    retsb           #   return dabs(x)
 129         emodd   %r0,$0,%r0,%r4,%r0      # %r4 + %r0:1 = scaled x^2
 130         emodd   %r2,$0,%r2,%r5,%r2      # %r5 + %r2:3 = scaled y^2
 131         addd2   %r2,%r0
 132         addl2   %r5,%r4
 133         cvtld   %r4,%r2
 134         addd2   %r2,%r0         # %r0:1 = scaled x^2 + y^2
 135         jmp     _C_LABEL(__libm_dsqrt_r5)+2
 136                                 # %r0:1 = dsqrt(x^2+y^2)/2^%r6
 137 retsb:
 138         rsb                     # error < 0.86 ulp