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