lib/libm/common_source/pow.c

   1 /*
   2  * Copyright (c) 1985, 1993
   3  *      The Regents of the University of California.  All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  * 1. Redistributions of source code must retain the above copyright
   9  *    notice, this list of conditions and the following disclaimer.
  10  * 2. Redistributions in binary form must reproduce the above copyright
  11  *    notice, this list of conditions and the following disclaimer in the
  12  *    documentation and/or other materials provided with the distribution.
  13  * 3. All advertising materials mentioning features or use of this software
  14  *    must display the following acknowledgement:
  15  *      This product includes software developed by the University of
  16  *      California, Berkeley and its contributors.
  17  * 4. Neither the name of the University nor the names of its contributors
  18  *    may be used to endorse or promote products derived from this software
  19  *    without specific prior written permission.
  20  *
  21  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31  * SUCH DAMAGE.
  32  */
  33
  34 #ifndef lint
  35 static char sccsid[] = "@(#)pow.c       8.1 (Berkeley) 6/4/93";
  36 #endif /* not lint */
  37
  38 /* POW(X,Y)
  39  * RETURN X**Y
  40  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  41  * CODED IN C BY K.C. NG, 1/8/85;
  42  * REVISED BY K.C. NG on 7/10/85.
  43  * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92.
  44  * Required system supported functions:
  45  *      scalb(x,n)
  46  *      logb(x)
  47  *      copysign(x,y)
  48  *      finite(x)
  49  *      drem(x,y)
  50  *
  51  * Required kernel functions:
  52  *      exp__D(a,c)                     exp(a + c) for |a| << |c|
  53  *      struct d_double dlog(x)         r.a + r.b, |r.b| < |r.a|
  54  *
  55  * Method
  56  *      1. Compute and return log(x) in three pieces:
  57  *              log(x) = n*ln2 + hi + lo,
  58  *         where n is an integer.
  59  *      2. Perform y*log(x) by simulating muti-precision arithmetic and
  60  *         return the answer in three pieces:
  61  *              y*log(x) = m*ln2 + hi + lo,
  62  *         where m is an integer.
  63  *      3. Return x**y = exp(y*log(x))
  64  *              = 2^m * ( exp(hi+lo) ).
  65  *
  66  * Special cases:
  67  *      (anything) ** 0  is 1 ;
  68  *      (anything) ** 1  is itself;
  69  *      (anything) ** NaN is NaN;
  70  *      NaN ** (anything except 0) is NaN;
  71  *      +(anything > 1) ** +INF is +INF;
  72  *      -(anything > 1) ** +INF is NaN;
  73  *      +-(anything > 1) ** -INF is +0;
  74  *      +-(anything < 1) ** +INF is +0;
  75  *      +(anything < 1) ** -INF is +INF;
  76  *      -(anything < 1) ** -INF is NaN;
  77  *      +-1 ** +-INF is NaN and signal INVALID;
  78  *      +0 ** +(anything except 0, NaN)  is +0;
  79  *      -0 ** +(anything except 0, NaN, odd integer)  is +0;
  80  *      +0 ** -(anything except 0, NaN)  is +INF and signal DIV-BY-ZERO;
  81  *      -0 ** -(anything except 0, NaN, odd integer)  is +INF with signal;
  82  *      -0 ** (odd integer) = -( +0 ** (odd integer) );
  83  *      +INF ** +(anything except 0,NaN) is +INF;
  84  *      +INF ** -(anything except 0,NaN) is +0;
  85  *      -INF ** (odd integer) = -( +INF ** (odd integer) );
  86  *      -INF ** (even integer) = ( +INF ** (even integer) );
  87  *      -INF ** -(anything except integer,NaN) is NaN with signal;
  88  *      -(x=anything) ** (k=integer) is (-1)**k * (x ** k);
  89  *      -(anything except 0) ** (non-integer) is NaN with signal;
  90  *
  91  * Accuracy:
  92  *      pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX,
  93  *      and a Zilog Z8000,
  94  *                      pow(integer,integer)
  95  *      always returns the correct integer provided it is representable.
  96  *      In a test run with 100,000 random arguments with 0 < x, y < 20.0
  97  *      on a VAX, the maximum observed error was 1.79 ulps (units in the
  98  *      last place).
  99  *
 100  * Constants :
 101  * The hexadecimal values are the intended ones for the following constants.
 102  * The decimal values may be used, provided that the compiler will convert
 103  * from decimal to binary accurately enough to produce the hexadecimal values
 104  * shown.
 105  */
 106
 107 #include <errno.h>
 108 #include <math.h>
 109
 110 #include "mathimpl.h"
 111
 112 #if (defined(vax) || defined(tahoe))
 113 #define TRUNC(x)        x = (double) (float) x
 114 #define _IEEE           0
 115 #else
 116 #define _IEEE           1
 117 #define endian          (((*(int *) &one)) ? 1 : 0)
 118 #define TRUNC(x)        *(((int *) &x)+endian) &= 0xf8000000
 119 #define infnan(x)       0.0
 120 #endif          /* vax or tahoe */
 121
 122 const static double zero=0.0, one=1.0, two=2.0, negone= -1.0;
 123
 124 static double pow_P __P((double, double));
 125
 126 double pow(x,y)
 127 double x,y;
 128 {
 129         double t;
 130         if (y==zero)
 131                 return (one);
 132         else if (y==one || (_IEEE && x != x))
 133                 return (x);             /* if x is NaN or y=1 */
 134         else if (_IEEE && y!=y)         /* if y is NaN */
 135                 return (y);
 136         else if (!finite(y))            /* if y is INF */
 137                 if ((t=fabs(x))==one)   /* +-1 ** +-INF is NaN */
 138                         return (y - y);
 139                 else if (t>one)
 140                         return ((y<0)? zero : ((x<zero)? y-y : y));
 141                 else
 142                         return ((y>0)? zero : ((x<0)? y-y : -y));
 143         else if (y==two)
 144                 return (x*x);
 145         else if (y==negone)
 146                 return (one/x);
 147     /* x > 0, x == +0 */
 148         else if (copysign(one, x) == one)
 149                 return (pow_P(x, y));
 150
 151     /* sign(x)= -1 */
 152         /* if y is an even integer */
 153         else if ( (t=drem(y,two)) == zero)
 154                 return (pow_P(-x, y));
 155
 156         /* if y is an odd integer */
 157         else if (copysign(t,one) == one)
 158                 return (-pow_P(-x, y));
 159
 160         /* Henceforth y is not an integer */
 161         else if (x==zero)       /* x is -0 */
 162                 return ((y>zero)? -x : one/(-x));
 163         else if (_IEEE)
 164                 return (zero/zero);
 165         else
 166                 return (infnan(EDOM));
 167 }
 168 /* kernel function for x >= 0 */
 169 static double
 170 #ifdef _ANSI_SOURCE
 171 pow_P(double x, double y)
 172 #else
 173 pow_P(x, y) double x, y;
 174 #endif
 175 {
 176         struct Double s, t, __log__D();
 177         double  __exp__D();
 178         volatile double huge = 1e300, tiny = 1e-300;
 179
 180         if (x == zero)
 181                 if (y > zero)
 182                         return (zero);
 183                 else if (_IEEE)
 184                         return (huge*huge);
 185                 else
 186                         return (infnan(ERANGE));
 187         if (x == one)
 188                 return (one);
 189         if (!finite(x))
 190                 if (y < zero)
 191                         return (zero);
 192                 else if (_IEEE)
 193                         return (huge*huge);
 194                 else
 195                         return (infnan(ERANGE));
 196         if (y >= 7e18)          /* infinity */
 197                 if (x < 1)
 198                         return(tiny*tiny);
 199                 else if (_IEEE)
 200                         return (huge*huge);
 201                 else
 202                         return (infnan(ERANGE));
 203
 204         /* Return exp(y*log(x)), using simulated extended */
 205         /* precision for the log and the multiply.        */
 206
 207         s = __log__D(x);
 208         t.a = y;
 209         TRUNC(t.a);
 210         t.b = y - t.a;
 211         t.b = s.b*y + t.b*s.a;
 212         t.a *= s.a;
 213         s.a = t.a + t.b;
 214         s.b = (t.a - s.a) + t.b;
 215         return (__exp__D(s.a, s.b));
 216 }