lib/libm/arch/vax/n_sqrt.S

   1 /*      $NetBSD: n_sqrt.S,v 1.11 2014/10/11 06:59:29 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  *      @(#)sqrt.s      8.1 (Berkeley) 6/4/93
  31  */
  32
  33 #include <machine/asm.h>
  34
  35 #ifdef WEAK_ALIAS
  36 WEAK_ALIAS(_sqrtl, sqrt)
  37 WEAK_ALIAS(sqrtl, sqrt)
  38 #endif
  39
  40 /*
  41  * double sqrt(arg)   revised August 15,1982
  42  * double arg;
  43  * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); }
  44  * if arg is a reserved operand it is returned as it is
  45  * W. Kahan's magic square root
  46  * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82
  47  *
  48  * entry points:_d_sqrt         address of double arg is on the stack
  49  *              _sqrt           double arg is on the stack
  50  */
  51         .set    EDOM,33
  52
  53 ENTRY(d_sqrt, 0x003c)           # save %r5,%r4,%r3,%r2
  54         movq    *4(%ap),%r0
  55         jbr     dsqrt2
  56
  57 ENTRY(sqrt, 0x003c)             # save %r5,%r4,%r3,%r2
  58         movq    4(%ap),%r0
  59
  60 dsqrt2: bicw3   $0x807f,%r0,%r2 # check exponent of input
  61         jeql    noexp           # biased exponent is zero -> 0.0 or reserved
  62         bsbb    __libm_dsqrt_r5_lcl
  63 noexp:  ret
  64
  65 /* **************************** internal procedure */
  66
  67         .hidden __libm_dsqrt_r5
  68 ALTENTRY(__libm_dsqrt_r5)
  69         halt
  70         halt
  71 __libm_dsqrt_r5_lcl:
  72                                 /* ENTRY POINT FOR cdabs and cdsqrt     */
  73                                 /* returns double square root scaled by */
  74                                 /* 2^%r6        */
  75
  76         movd    %r0,%r4
  77         jleq    nonpos          # argument is not positive
  78         movzwl  %r4,%r2
  79         ashl    $-1,%r2,%r0
  80         addw2   $0x203c,%r0     # %r0 has magic initial approximation
  81 /*
  82  * Do two steps of Heron's rule
  83  * ((arg/guess) + guess) / 2 = better guess
  84  */
  85         divf3   %r0,%r4,%r2
  86         addf2   %r2,%r0
  87         subw2   $0x80,%r0       # divide by two
  88
  89         divf3   %r0,%r4,%r2
  90         addf2   %r2,%r0
  91         subw2   $0x80,%r0       # divide by two
  92
  93 /* Scale argument and approximation to prevent over/underflow */
  94
  95         bicw3   $0x807f,%r4,%r1
  96         subw2   $0x4080,%r1             # %r1 contains scaling factor
  97         subw2   %r1,%r4
  98         movl    %r0,%r2
  99         subw2   %r1,%r2
 100
 101 /* Cubic step
 102  *
 103  * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation,
 104  * a is approximation, and n is the original argument.
 105  * (let s be scale factor in the following comments)
 106  */
 107         clrl    %r1
 108         clrl    %r3
 109         muld2   %r0,%r2                 # %r2:%r3 = a*a/s
 110         subd2   %r2,%r4                 # %r4:%r5 = n/s - a*a/s
 111         addw2   $0x100,%r2              # %r2:%r3 = 4*a*a/s
 112         addd2   %r4,%r2                 # %r2:%r3 = n/s + 3*a*a/s
 113         muld2   %r0,%r4                 # %r4:%r5 = a*n/s - a*a*a/s
 114         divd2   %r2,%r4                 # %r4:%r5 = a*(n-a*a)/(n+3*a*a)
 115         addw2   $0x80,%r4               # %r4:%r5 = 2*a*(n-a*a)/(n+3*a*a)
 116         addd2   %r4,%r0                 # %r0:%r1 = a + 2*a*(n-a*a)/(n+3*a*a)
 117         rsb                             # DONE!
 118 nonpos:
 119         jneq    negarg
 120         ret                             # argument and root are zero
 121 negarg:
 122         pushl   $EDOM
 123         calls   $1,_C_LABEL(infnan)     # generate the reserved op fault
 124         ret
 125
 126 ENTRY(sqrtf, 0)
 127         cvtfd   4(%ap),-(%sp)
 128         calls   $2,_C_LABEL(sqrt)
 129         cvtdf   %r0,%r0
 130         ret