lib/libm/noieee_src/n_pow.c

   1 /*      $NetBSD: n_pow.c,v 1.11 2014/10/11 07:19:27 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
  31 #ifndef lint
  32 #if 0
  33 static char sccsid[] = "@(#)pow.c       8.1 (Berkeley) 6/4/93";
  34 #endif
  35 #endif /* not lint */
  36
  37 /* POW(X,Y)
  38  * RETURN X**Y
  39  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  40  * CODED IN C BY K.C. NG, 1/8/85;
  41  * REVISED BY K.C. NG on 7/10/85.
  42  * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92.
  43  * Required system supported functions:
  44  *      scalb(x,n)
  45  *      logb(x)
  46  *      copysign(x,y)
  47  *      finite(x)
  48  *      drem(x,y)
  49  *
  50  * Required kernel functions:
  51  *      exp__D(a,c)                     exp(a + c) for |a| << |c|
  52  *      struct d_double dlog(x)         r.a + r.b, |r.b| < |r.a|
  53  *
  54  * Method
  55  *      1. Compute and return log(x) in three pieces:
  56  *              log(x) = n*ln2 + hi + lo,
  57  *         where n is an integer.
  58  *      2. Perform y*log(x) by simulating muti-precision arithmetic and
  59  *         return the answer in three pieces:
  60  *              y*log(x) = m*ln2 + hi + lo,
  61  *         where m is an integer.
  62  *      3. Return x**y = exp(y*log(x))
  63  *              = 2^m * ( exp(hi+lo) ).
  64  *
  65  * Special cases:
  66  *      (anything) ** 0  is 1 ;
  67  *      (anything) ** 1  is itself;
  68  *      (anything) ** NaN is NaN;
  69  *      NaN ** (anything except 0) is NaN;
  70  *      +(anything > 1) ** +INF is +INF;
  71  *      -(anything > 1) ** +INF is NaN;
  72  *      +-(anything > 1) ** -INF is +0;
  73  *      +-(anything < 1) ** +INF is +0;
  74  *      +(anything < 1) ** -INF is +INF;
  75  *      -(anything < 1) ** -INF is NaN;
  76  *      +-1 ** +-INF is NaN and signal INVALID;
  77  *      +0 ** +(anything except 0, NaN)  is +0;
  78  *      -0 ** +(anything except 0, NaN, odd integer)  is +0;
  79  *      +0 ** -(anything except 0, NaN)  is +INF and signal DIV-BY-ZERO;
  80  *      -0 ** -(anything except 0, NaN, odd integer)  is +INF with signal;
  81  *      -0 ** (odd integer) = -( +0 ** (odd integer) );
  82  *      +INF ** +(anything except 0,NaN) is +INF;
  83  *      +INF ** -(anything except 0,NaN) is +0;
  84  *      -INF ** (odd integer) = -( +INF ** (odd integer) );
  85  *      -INF ** (even integer) = ( +INF ** (even integer) );
  86  *      -INF ** -(anything except integer,NaN) is NaN with signal;
  87  *      -(x=anything) ** (k=integer) is (-1)**k * (x ** k);
  88  *      -(anything except 0) ** (non-integer) is NaN with signal;
  89  *
  90  * Accuracy:
  91  *      pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX,
  92  *      and a Zilog Z8000,
  93  *                      pow(integer,integer)
  94  *      always returns the correct integer provided it is representable.
  95  *      In a test run with 100,000 random arguments with 0 < x, y < 20.0
  96  *      on a VAX, the maximum observed error was 1.79 ulps (units in the
  97  *      last place).
  98  *
  99  * Constants :
 100  * The hexadecimal values are the intended ones for the following constants.
 101  * The decimal values may be used, provided that the compiler will convert
 102  * from decimal to binary accurately enough to produce the hexadecimal values
 103  * shown.
 104  */
 105
 106 #include <errno.h>
 107 #include <math.h>
 108
 109 #include "mathimpl.h"
 110
 111 #if (defined(__vax__) || defined(tahoe))
 112 #define TRUNC(x)        x = (double) (float) x
 113 #define _IEEE           0
 114 #else
 115 #define _IEEE           1
 116 #define endian          (((*(int *) &one)) ? 1 : 0)
 117 #define TRUNC(x)        *(((int *) &x)+endian) &= 0xf8000000
 118 #define infnan(x)       0.0
 119 #endif          /* __vax__ or tahoe */
 120
 121 static const double zero=0.0, one=1.0, two=2.0, negone= -1.0;
 122
 123 static double pow_P (double, double);
 124
 125 #ifdef __weak_alias
 126 __weak_alias(_powf, powf);
 127 __weak_alias(_pow, pow);
 128 __weak_alias(_powl, pow);
 129 __weak_alias(powl, pow);
 130 #endif
 131
 132 float
 133 powf(float x, float y)
 134 {
 135    return pow((double) x, (double) (y));
 136 }
 137
 138 double
 139 pow(double x, double y)
 140 {
 141         double t;
 142         if (y==zero)
 143                 return (one);
 144         else if (y==one || (_IEEE && x != x))
 145                 return (x);             /* if x is NaN or y=1 */
 146         else if (_IEEE && y!=y)         /* if y is NaN */
 147                 return (y);
 148         else if (!finite(y))            /* if y is INF */
 149                 if ((t=fabs(x))==one)   /* +-1 ** +-INF is NaN */
 150                         return (y - y);
 151                 else if (t>one)
 152                         return ((y<0)? zero : ((x<zero)? y-y : y));
 153                 else
 154                         return ((y>0)? zero : ((x<0)? y-y : -y));
 155         else if (y==two)
 156                 return (x*x);
 157         else if (y==negone)
 158                 return (one/x);
 159     /* x > 0, x == +0 */
 160         else if (copysign(one, x) == one)
 161                 return (pow_P(x, y));
 162
 163     /* sign(x)= -1 */
 164         /* if y is an even integer */
 165         else if ( (t=drem(y,two)) == zero)
 166                 return (pow_P(-x, y));
 167
 168         /* if y is an odd integer */
 169         else if (copysign(t,one) == one)
 170                 return (-pow_P(-x, y));
 171
 172         /* Henceforth y is not an integer */
 173         else if (x==zero)       /* x is -0 */
 174                 return ((y>zero)? -x : one/(-x));
 175         else if (_IEEE)
 176                 return (zero/zero);
 177         else
 178                 return (infnan(EDOM));
 179 }
 180
 181 /* kernel function for x >= 0 */
 182 static double
 183 pow_P(double x, double y)
 184 {
 185         struct Double s, t;
 186         double  huge = _HUGE, tiny = _TINY;
 187
 188         if (x == zero) {
 189                 if (y > zero)
 190                         return (zero);
 191                 else if (_IEEE)
 192                         return (huge*huge);
 193                 else
 194                         return (infnan(ERANGE));
 195         }
 196         if (x == one)
 197                 return (one);
 198         if (!finite(x)) {
 199                 if (y < zero)
 200                         return (zero);
 201                 else if (_IEEE)
 202                         return (huge*huge);
 203                 else
 204                         return (infnan(ERANGE));
 205         }
 206         if (y >= 7e18) {        /* infinity */
 207                 if (x < 1)
 208                         return(tiny*tiny);
 209                 else if (_IEEE)
 210                         return (huge*huge);
 211                 else
 212                         return (infnan(ERANGE));
 213         }
 214
 215         /* Return exp(y*log(x)), using simulated extended */
 216         /* precision for the log and the multiply.        */
 217
 218         s = __log__D(x);
 219         t.a = y;
 220         TRUNC(t.a);
 221         t.b = y - t.a;
 222         t.b = s.b*y + t.b*s.a;
 223         t.a *= s.a;
 224         s.a = t.a + t.b;
 225         s.b = (t.a - s.a) + t.b;
 226         return (__exp__D(s.a, s.b));
 227 }