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  * $DragonFly: src/lib/libm/common_source/Attic/pow.c,v 1.4 2003/11/12 20:21:31 eirikn Exp $
  35  */
  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 const static double zero=0.0, one=1.0, two=2.0, negone= -1.0;
 122
 123 static double pow_P (double, double);
 124
 125 double pow(x,y)
 126 double x,y;
 127 {
 128         double t;
 129         if (y==zero)
 130                 return (one);
 131         else if (y==one || (_IEEE && x != x))
 132                 return (x);             /* if x is NaN or y=1 */
 133         else if (_IEEE && y!=y)         /* if y is NaN */
 134                 return (y);
 135         else if (!finite(y))            /* if y is INF */
 136                 if ((t=fabs(x))==one)   /* +-1 ** +-INF is NaN */
 137                         return (y - y);
 138                 else if (t>one)
 139                         return ((y<0)? zero : ((x<zero)? y-y : y));
 140                 else
 141                         return ((y>0)? zero : ((x<0)? y-y : -y));
 142         else if (y==two)
 143                 return (x*x);
 144         else if (y==negone)
 145                 return (one/x);
 146     /* x > 0, x == +0 */
 147         else if (copysign(one, x) == one)
 148                 return (pow_P(x, y));
 149
 150     /* sign(x)= -1 */
 151         /* if y is an even integer */
 152         else if ( (t=drem(y,two)) == zero)
 153                 return (pow_P(-x, y));
 154
 155         /* if y is an odd integer */
 156         else if (copysign(t,one) == one)
 157                 return (-pow_P(-x, y));
 158
 159         /* Henceforth y is not an integer */
 160         else if (x==zero)       /* x is -0 */
 161                 return ((y>zero)? -x : one/(-x));
 162         else if (_IEEE)
 163                 return (zero/zero);
 164         else
 165                 return (infnan(EDOM));
 166 }
 167 /* kernel function for x >= 0 */
 168 static double
 169 #ifdef _ANSI_SOURCE
 170 pow_P(double x, double y)
 171 #else
 172 pow_P(x, y) double x, y;
 173 #endif
 174 {
 175         struct Double s, t, __log__D();
 176         double  __exp__D();
 177         volatile double huge = 1e300, tiny = 1e-300;
 178
 179         if (x == zero)
 180                 if (y > zero)
 181                         return (zero);
 182                 else if (_IEEE)
 183                         return (huge*huge);
 184                 else
 185                         return (infnan(ERANGE));
 186         if (x == one)
 187                 return (one);
 188         if (!finite(x))
 189                 if (y < zero)
 190                         return (zero);
 191                 else if (_IEEE)
 192                         return (huge*huge);
 193                 else
 194                         return (infnan(ERANGE));
 195         if (y >= 7e18)          /* infinity */
 196                 if (x < 1)
 197                         return(tiny*tiny);
 198                 else if (_IEEE)
 199                         return (huge*huge);
 200                 else
 201                         return (infnan(ERANGE));
 202
 203         /* Return exp(y*log(x)), using simulated extended */
 204         /* precision for the log and the multiply.        */
 205
 206         s = __log__D(x);
 207         t.a = y;
 208         TRUNC(t.a);
 209         t.b = y - t.a;
 210         t.b = s.b*y + t.b*s.a;
 211         t.a *= s.a;
 212         s.a = t.a + t.b;
 213         s.b = (t.a - s.a) + t.b;
 214         return (__exp__D(s.a, s.b));
 215 }