/* $NetBSD: n_pow.c,v 1.11 2014/10/11 07:19:27 martin Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#ifndef lint
#if 0
static char sccsid[] = "@(#)pow.c 8.1 (Berkeley) 6/4/93";
#endif
#endif /* not lint */
/* POW(X,Y)
* RETURN X**Y
* DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
* CODED IN C BY K.C. NG, 1/8/85;
* REVISED BY K.C. NG on 7/10/85.
* KERNEL pow_P() REPLACED BY P. McILROY 7/22/92.
* Required system supported functions:
* scalb(x,n)
* logb(x)
* copysign(x,y)
* finite(x)
* drem(x,y)
*
* Required kernel functions:
* exp__D(a,c) exp(a + c) for |a| << |c|
* struct d_double dlog(x) r.a + r.b, |r.b| < |r.a|
*
* Method
* 1. Compute and return log(x) in three pieces:
* log(x) = n*ln2 + hi + lo,
* where n is an integer.
* 2. Perform y*log(x) by simulating muti-precision arithmetic and
* return the answer in three pieces:
* y*log(x) = m*ln2 + hi + lo,
* where m is an integer.
* 3. Return x**y = exp(y*log(x))
* = 2^m * ( exp(hi+lo) ).
*
* Special cases:
* (anything) ** 0 is 1 ;
* (anything) ** 1 is itself;
* (anything) ** NaN is NaN;
* NaN ** (anything except 0) is NaN;
* +(anything > 1) ** +INF is +INF;
* -(anything > 1) ** +INF is NaN;
* +-(anything > 1) ** -INF is +0;
* +-(anything < 1) ** +INF is +0;
* +(anything < 1) ** -INF is +INF;
* -(anything < 1) ** -INF is NaN;
* +-1 ** +-INF is NaN and signal INVALID;
* +0 ** +(anything except 0, NaN) is +0;
* -0 ** +(anything except 0, NaN, odd integer) is +0;
* +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO;
* -0 ** -(anything except 0, NaN, odd integer) is +INF with signal;
* -0 ** (odd integer) = -( +0 ** (odd integer) );
* +INF ** +(anything except 0,NaN) is +INF;
* +INF ** -(anything except 0,NaN) is +0;
* -INF ** (odd integer) = -( +INF ** (odd integer) );
* -INF ** (even integer) = ( +INF ** (even integer) );
* -INF ** -(anything except integer,NaN) is NaN with signal;
* -(x=anything) ** (k=integer) is (-1)**k * (x ** k);
* -(anything except 0) ** (non-integer) is NaN with signal;
*
* Accuracy:
* pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX,
* and a Zilog Z8000,
* pow(integer,integer)
* always returns the correct integer provided it is representable.
* In a test run with 100,000 random arguments with 0 < x, y < 20.0
* on a VAX, the maximum observed error was 1.79 ulps (units in the
* last place).
*
* Constants :
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
#include <errno.h>
#include <math.h>
#include "mathimpl.h"
#if (defined(__vax__) || defined(tahoe))
#define TRUNC(x) x = (double) (float) x
#define _IEEE 0
#else
#define _IEEE 1
#define endian (((*(int *) &one)) ? 1 : 0)
#define TRUNC(x) *(((int *) &x)+endian) &= 0xf8000000
#define infnan(x) 0.0
#endif /* __vax__ or tahoe */
static const double zero=0.0, one=1.0, two=2.0, negone= -1.0;
static double pow_P (double, double);
#ifdef __weak_alias
__weak_alias(_powf, powf);
__weak_alias(_pow, pow);
__weak_alias(_powl, pow);
__weak_alias(powl, pow);
#endif
float
powf(float x, float y)
{
return pow((double) x, (double) (y));
}
double
pow(double x, double y)
{
double t;
if (y==zero)
return (one);
else if (y==one || (_IEEE && x != x))
return (x); /* if x is NaN or y=1 */
else if (_IEEE && y!=y) /* if y is NaN */
return (y);
else if (!finite(y)) /* if y is INF */
if ((t=fabs(x))==one) /* +-1 ** +-INF is NaN */
return (y - y);
else if (t>one)
return ((y<0)? zero : ((x<zero)? y-y : y));
else
return ((y>0)? zero : ((x<0)? y-y : -y));
else if (y==two)
return (x*x);
else if (y==negone)
return (one/x);
/* x > 0, x == +0 */
else if (copysign(one, x) == one)
return (pow_P(x, y));
/* sign(x)= -1 */
/* if y is an even integer */
else if ( (t=drem(y,two)) == zero)
return (pow_P(-x, y));
/* if y is an odd integer */
else if (copysign(t,one) == one)
return (-pow_P(-x, y));
/* Henceforth y is not an integer */
else if (x==zero) /* x is -0 */
return ((y>zero)? -x : one/(-x));
else if (_IEEE)
return (zero/zero);
else
return (infnan(EDOM));
}
/* kernel function for x >= 0 */
static double
pow_P(double x, double y)
{
struct Double s, t;
double huge = _HUGE, tiny = _TINY;
if (x == zero) {
if (y > zero)
return (zero);
else if (_IEEE)
return (huge*huge);
else
return (infnan(ERANGE));
}
if (x == one)
return (one);
if (!finite(x)) {
if (y < zero)
return (zero);
else if (_IEEE)
return (huge*huge);
else
return (infnan(ERANGE));
}
if (y >= 7e18) { /* infinity */
if (x < 1)
return(tiny*tiny);
else if (_IEEE)
return (huge*huge);
else
return (infnan(ERANGE));
}
/* Return exp(y*log(x)), using simulated extended */
/* precision for the log and the multiply. */
s = __log__D(x);
t.a = y;
TRUNC(t.a);
t.b = y - t.a;
t.b = s.b*y + t.b*s.a;
t.a *= s.a;
s.a = t.a + t.b;
s.b = (t.a - s.a) + t.b;
return (__exp__D(s.a, s.b));
}