lib/libm/arch/vax/n_sqrt.S

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