lib/libm/noieee_src/n_pow.c

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