lib/libm/arch/vax/n_cabs.S

   1 /*      $NetBSD: n_cabs.S,v 1.7 2014/10/10 20:58:09 martin 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 WEAK_ALIAS(hypotl, _hypot)
  59 WEAK_ALIAS(_hypotl, _hypot)
  60 #endif
  61
  62 ALTENTRY(cabs)
  63 ENTRY(_hypot, 0x8040)           # save %r6, enable floating overflow
  64         movq    4(%ap),%r0      # %r0:1 = x
  65         movq    12(%ap),%r2     # %r2:3 = y
  66         jbr     cabs2
  67
  68 /*      entry for Fortran use, call by:   d = abs(z) */
  69 ENTRY(z_abs, 0x8040)            # save %r6, enable floating overflow
  70         movl    4(%ap),%r2      # indirect addressing is necessary here
  71         movq    (%r2)+,%r0      # %r0:1 = x
  72         movq    (%r2),%r2               # %r2:3 = y
  73
  74 cabs2:
  75         bicw3   $0x7f,%r0,%r4   # %r4 has signed biased exp of x
  76         cmpw    $0x8000,%r4
  77         jeql    return          # x is a reserved operand, so return it
  78         bicw3   $0x7f,%r2,%r5   # %r5 has signed biased exp of y
  79         cmpw    $0x8000,%r5
  80         jneq    cont            /* y isn't a reserved operand */
  81         movq    %r2,%r0         /* return y if it's reserved */
  82         ret
  83
  84 cont:
  85         bsbb    regs_set        # %r0:1 = dsqrt(x^2+y^2)/2^%r6
  86         addw2   %r6,%r0         # unscaled cdabs in %r0:1
  87         jvc     return          # unless it overflows
  88         subw2   $0x80,%r0       # halve %r0 to get meaningful overflow
  89         addd2   %r0,%r0         # overflow; %r0 is half of true abs value
  90 return:
  91         ret
  92
  93 ENTRY(__libm_cdabs_r6,0)        # ENTRY POINT for cdsqrt
  94                                 # calculates a scaled (factor in %r6)
  95                                 # complex absolute value
  96
  97         movq    (%r4)+,%r0      # %r0:%r1 = x via indirect addressing
  98         movq    (%r4),%r2               # %r2:%r3 = y via indirect addressing
  99
 100         bicw3   $0x7f,%r0,%r5   # %r5 has signed biased exp of x
 101         cmpw    $0x8000,%r5
 102         jeql    cdreserved      # x is a reserved operand
 103         bicw3   $0x7f,%r2,%r5   # %r5 has signed biased exp of y
 104         cmpw    $0x8000,%r5
 105         jneq    regs_set        /* y isn't a reserved operand either? */
 106
 107 cdreserved:
 108         movl    *4(%ap),%r4     # %r4 -> (u,v), if x or y is reserved
 109         movq    %r0,(%r4)+      # copy u and v as is and return
 110         movq    %r2,(%r4)               # (again addressing is indirect)
 111         ret
 112
 113 regs_set:
 114         bicw2   $0x8000,%r0     # %r0:%r1 = dabs(x)
 115         bicw2   $0x8000,%r2     # %r2:%r3 = dabs(y)
 116         cmpw    %r0,%r2
 117         jgeq    ordered
 118         movq    %r0,%r4
 119         movq    %r2,%r0
 120         movq    %r4,%r2         # force y's exp <= x's exp
 121 ordered:
 122         bicw3   $0x7f,%r0,%r6   # %r6 = exponent(x) + bias(129)
 123         jeql    retsb           # if x = y = 0 then cdabs(x,y) = 0
 124         subw2   $0x4780,%r6     # %r6 = exponent(x) - 14
 125         subw2   %r6,%r0         # 2^14 <= scaled x < 2^15
 126         bitw    $0xff80,%r2
 127         jeql    retsb           # if y = 0 return dabs(x)
 128         subw2   %r6,%r2
 129         cmpw    $0x3780,%r2     # if scaled y < 2^-18
 130         jgtr    retsb           #   return dabs(x)
 131         emodd   %r0,$0,%r0,%r4,%r0      # %r4 + %r0:1 = scaled x^2
 132         emodd   %r2,$0,%r2,%r5,%r2      # %r5 + %r2:3 = scaled y^2
 133         addd2   %r2,%r0
 134         addl2   %r5,%r4
 135         cvtld   %r4,%r2
 136         addd2   %r2,%r0         # %r0:1 = scaled x^2 + y^2
 137         jmp     _C_LABEL(__libm_dsqrt_r5)+2
 138                                 # %r0:1 = dsqrt(x^2+y^2)/2^%r6
 139 retsb:
 140         rsb                     # error < 0.86 ulp