lib/libm/arch/vax/n_cbrt.S

   1 /*      $NetBSD: n_cbrt.S,v 1.7 2013/11/24 15:12:10 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  *      @(#)cbrt.s      8.1 (Berkeley) 6/4/93
  31  */
  32
  33 #include <machine/asm.h>
  34
  35 /*
  36  * double cbrt(double arg)
  37  * W. Kahan, 10/13/80. revised 1/13/84 for keeping sign symmetry
  38  * error check by E LeBlanc, 8/18/82
  39  * Revised and tested by K.C. Ng, 5/2/85
  40  * Max error less than 0.667 ulps (unit in the last places)
  41  */
  42
  43 ENTRY(cbrtf, 0)
  44         cvtfd   4(%ap),-(%sp)
  45         calls   $2,_C_LABEL(cbrt)
  46         cvtdf   %r0,%r0
  47         ret
  48
  49 #ifdef WEAK_ALIAS
  50 WEAK_ALIAS(cbrtl, cbrt)
  51 #endif
  52
  53 ALTENTRY(cbrt)
  54 ENTRY(d_cbrt, 0x00c0)           # save %r6 & %r7
  55         movq    4(%ap),%r0      # %r0 = argument x
  56         jbr     dcbrt2
  57
  58 ENTRY(dcbrt_, 0x00c0)           # save %r6 & %r7
  59         movq    *4(%ap),%r0     # %r0 = argument x
  60
  61 dcbrt2: bicw3   $0x807f,%r0,%r2 # biased exponent of x
  62         jeql    return          # dcbrt(0)=0  dcbrt(res)=res. operand
  63         bicw3   $0x7fff,%r0,%ap # ap has sign(x)
  64         xorw2   %ap,%r0         # %r0 is abs(x)
  65         movl    %r0,%r2         # %r2 has abs(x)
  66         rotl    $16,%r2,%r2     # %r2 = |x| with bits unscrambled
  67         divl2   $3,%r2          # rough dcbrt with bias/3
  68         addl2   B,%r2           # restore bias, diminish fraction
  69         rotl    $16,%r2,%r2     # %r2=|q|=|dcbrt| to 5 bits
  70         mulf3   %r2,%r2,%r3     # %r3 =qq
  71         divf2   %r0,%r3         # %r3 = qq/x
  72         mulf2   %r2,%r3
  73         addf2   C,%r3           # %r3 = s = C + qqq/x
  74         divf3   %r3,D,%r4               # %r4 = D/s
  75         addf2   E,%r4
  76         addf2   %r4,%r3         # %r3 = s + E + D/s
  77         divf3   %r3,F,%r3               # %r3 = F / (s + E + D/s)
  78         addf2   G,%r3           # %r3 = G + F / (s + E + D/s)
  79         mulf2   %r3,%r2         # %r2 = q%r3 = new q to 23 bits
  80         clrl    %r3             # %r2:%r3 = q as double float
  81         muld3   %r2,%r2,%r4     # %r4:%r5 = qq exactly
  82         divd2   %r4,%r0         # %r0:%r1 = x/(q*q) rounded
  83         subd3   %r2,%r0,%r6     # %r6:%r7 = x/(q*q) - q exactly
  84         movq    %r2,%r4         # %r4:%r5 = q
  85         addw2   $0x80,%r4       # %r4:%r5 = 2 * q
  86         addd2   %r0,%r4         # %r4:%r5 = 2*q + x/(q*q)
  87         divd2   %r4,%r6         # %r6:%r7 = (x/(q*q)-q)/(2*q+x/(q*q))
  88         muld2   %r2,%r6         # %r6:%r7 = q*(x/(q*q)-q)/(2*q+x/(q*q))
  89         addd3   %r6,%r2,%r0     # %r0:%r1 = q + %r6:%r7
  90         bisw2   %ap,%r0         # restore the sign bit
  91 return:
  92         ret                     # error less than 0.667 ulps
  93
  94         _ALIGN_TEXT
  95 B :     .long            721142941              # (86-0.03306235651)*(2^23)
  96 C :     .float          0f0.5428571429          # 19/35
  97 D :     .float          0f-0.7053061224         # -864/1225
  98 E :     .float          0f1.414285714           # 99/70
  99 F :     .float          0f1.607142857           # 45/28
 100 G :     .float          0f0.3571428571          # 5/14