2000-02-18 03:39:52 +08:00
|
|
|
/* wf_pow.c -- float version of w_pow.c.
|
|
|
|
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
|
|
|
*/
|
|
|
|
|
|
|
|
/*
|
|
|
|
* ====================================================
|
|
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
|
|
*
|
|
|
|
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
|
|
|
* Permission to use, copy, modify, and distribute this
|
|
|
|
* software is freely granted, provided that this notice
|
|
|
|
* is preserved.
|
|
|
|
* ====================================================
|
|
|
|
*/
|
|
|
|
|
|
|
|
/*
|
|
|
|
* wrapper powf(x,y) return x**y
|
|
|
|
*/
|
|
|
|
|
|
|
|
#include "fdlibm.h"
|
New expf, exp2f, logf, log2f and powf implementations
Based on code from https://github.com/ARM-software/optimized-routines/
This patch adds a highly optimized generic implementation of expf,
exp2f, logf, log2f and powf. The new functions are not only
faster (6x for powf!), but are also smaller and more accurate.
In order to achieve this, the algorithm uses double precision
arithmetic for accuracy, avoids divisions and uses small table
lookups to minimize the polynomials. Special cases are handled
inline to avoid the unnecessary overhead of wrapper functions and
set errno to POSIX requirements.
The new functions are added under newlib/libm/common, but the old
implementations are kept (in newlib/libm/math) for non-IEEE or
pre-C99 systems. Targets can enable the new math code by defining
__OBSOLETE_MATH_DEFAULT to 0 in newlib/libc/include/machine/ieeefp.h,
users can override the default by defining __OBSOLETE_MATH.
Currently the new code is enabled for AArch64 and AArch32 with VFP.
Targets with a single precision FPU may still prefer the old
implementation.
libm.a size changes:
arm: -1692
arm/thumb/v7-a/nofp: -878
arm/thumb/v7-a+fp/hard: -864
arm/thumb/v7-a+fp/softfp: -908
aarch64: -1476
2017-05-25 23:41:38 +08:00
|
|
|
#if __OBSOLETE_MATH
|
2000-02-18 03:39:52 +08:00
|
|
|
#include <errno.h>
|
|
|
|
|
|
|
|
#ifdef __STDC__
|
|
|
|
float powf(float x, float y) /* wrapper powf */
|
|
|
|
#else
|
|
|
|
float powf(x,y) /* wrapper powf */
|
|
|
|
float x,y;
|
|
|
|
#endif
|
|
|
|
{
|
|
|
|
#ifdef _IEEE_LIBM
|
|
|
|
return __ieee754_powf(x,y);
|
|
|
|
#else
|
|
|
|
float z;
|
|
|
|
struct exception exc;
|
|
|
|
z=__ieee754_powf(x,y);
|
2009-07-10 01:04:56 +08:00
|
|
|
if(_LIB_VERSION == _IEEE_|| isnan(y)) return z;
|
|
|
|
if(isnan(x)) {
|
2000-02-18 03:39:52 +08:00
|
|
|
if(y==(float)0.0) {
|
|
|
|
/* powf(NaN,0.0) */
|
|
|
|
/* error only if _LIB_VERSION == _SVID_ & _XOPEN_ */
|
|
|
|
exc.type = DOMAIN;
|
|
|
|
exc.name = "powf";
|
|
|
|
exc.err = 0;
|
|
|
|
exc.arg1 = (double)x;
|
|
|
|
exc.arg2 = (double)y;
|
2007-05-02 02:42:15 +08:00
|
|
|
exc.retval = 1.0;
|
2000-02-18 03:39:52 +08:00
|
|
|
if (_LIB_VERSION == _IEEE_ ||
|
|
|
|
_LIB_VERSION == _POSIX_) exc.retval = 1.0;
|
|
|
|
else if (!matherr(&exc)) {
|
|
|
|
errno = EDOM;
|
|
|
|
}
|
|
|
|
if (exc.err != 0)
|
|
|
|
errno = exc.err;
|
|
|
|
return (float)exc.retval;
|
|
|
|
} else
|
|
|
|
return z;
|
|
|
|
}
|
|
|
|
if(x==(float)0.0){
|
|
|
|
if(y==(float)0.0) {
|
|
|
|
/* powf(0.0,0.0) */
|
|
|
|
/* error only if _LIB_VERSION == _SVID_ */
|
|
|
|
exc.type = DOMAIN;
|
|
|
|
exc.name = "powf";
|
|
|
|
exc.err = 0;
|
|
|
|
exc.arg1 = (double)x;
|
|
|
|
exc.arg2 = (double)y;
|
|
|
|
exc.retval = 0.0;
|
|
|
|
if (_LIB_VERSION != _SVID_) exc.retval = 1.0;
|
|
|
|
else if (!matherr(&exc)) {
|
|
|
|
errno = EDOM;
|
|
|
|
}
|
|
|
|
if (exc.err != 0)
|
|
|
|
errno = exc.err;
|
|
|
|
return (float)exc.retval;
|
|
|
|
}
|
|
|
|
if(finitef(y)&&y<(float)0.0) {
|
|
|
|
/* 0**neg */
|
|
|
|
exc.type = DOMAIN;
|
|
|
|
exc.name = "powf";
|
|
|
|
exc.err = 0;
|
|
|
|
exc.arg1 = (double)x;
|
|
|
|
exc.arg2 = (double)y;
|
|
|
|
if (_LIB_VERSION == _SVID_)
|
|
|
|
exc.retval = 0.0;
|
|
|
|
else
|
|
|
|
exc.retval = -HUGE_VAL;
|
|
|
|
if (_LIB_VERSION == _POSIX_)
|
|
|
|
errno = EDOM;
|
|
|
|
else if (!matherr(&exc)) {
|
|
|
|
errno = EDOM;
|
|
|
|
}
|
|
|
|
if (exc.err != 0)
|
|
|
|
errno = exc.err;
|
|
|
|
return (float)exc.retval;
|
|
|
|
}
|
|
|
|
return z;
|
|
|
|
}
|
|
|
|
if(!finitef(z)) {
|
|
|
|
if(finitef(x)&&finitef(y)) {
|
2009-07-10 01:04:56 +08:00
|
|
|
if(isnan(z)) {
|
2000-02-18 03:39:52 +08:00
|
|
|
/* neg**non-integral */
|
|
|
|
exc.type = DOMAIN;
|
|
|
|
exc.name = "powf";
|
|
|
|
exc.err = 0;
|
|
|
|
exc.arg1 = (double)x;
|
|
|
|
exc.arg2 = (double)y;
|
|
|
|
if (_LIB_VERSION == _SVID_)
|
|
|
|
exc.retval = 0.0;
|
|
|
|
else
|
|
|
|
exc.retval = 0.0/0.0; /* X/Open allow NaN */
|
|
|
|
if (_LIB_VERSION == _POSIX_)
|
|
|
|
errno = EDOM;
|
|
|
|
else if (!matherr(&exc)) {
|
|
|
|
errno = EDOM;
|
|
|
|
}
|
|
|
|
if (exc.err != 0)
|
|
|
|
errno = exc.err;
|
|
|
|
return (float)exc.retval;
|
|
|
|
} else {
|
|
|
|
/* powf(x,y) overflow */
|
|
|
|
exc.type = OVERFLOW;
|
|
|
|
exc.name = "powf";
|
|
|
|
exc.err = 0;
|
|
|
|
exc.arg1 = (double)x;
|
|
|
|
exc.arg2 = (double)y;
|
|
|
|
if (_LIB_VERSION == _SVID_) {
|
|
|
|
exc.retval = HUGE;
|
|
|
|
y *= 0.5;
|
2017-12-13 03:38:01 +08:00
|
|
|
if(x<0.0f&&rintf(y)!=y) exc.retval = -HUGE;
|
2000-02-18 03:39:52 +08:00
|
|
|
} else {
|
|
|
|
exc.retval = HUGE_VAL;
|
|
|
|
y *= 0.5;
|
2017-12-13 03:38:01 +08:00
|
|
|
if(x<0.0f&&rintf(y)!=y) exc.retval = -HUGE_VAL;
|
2000-02-18 03:39:52 +08:00
|
|
|
}
|
|
|
|
if (_LIB_VERSION == _POSIX_)
|
|
|
|
errno = ERANGE;
|
|
|
|
else if (!matherr(&exc)) {
|
|
|
|
errno = ERANGE;
|
|
|
|
}
|
|
|
|
if (exc.err != 0)
|
|
|
|
errno = exc.err;
|
|
|
|
return (float)exc.retval;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
if(z==(float)0.0&&finitef(x)&&finitef(y)) {
|
|
|
|
/* powf(x,y) underflow */
|
|
|
|
exc.type = UNDERFLOW;
|
|
|
|
exc.name = "powf";
|
|
|
|
exc.err = 0;
|
|
|
|
exc.arg1 = (double)x;
|
|
|
|
exc.arg2 = (double)y;
|
|
|
|
exc.retval = 0.0;
|
|
|
|
if (_LIB_VERSION == _POSIX_)
|
|
|
|
errno = ERANGE;
|
|
|
|
else if (!matherr(&exc)) {
|
|
|
|
errno = ERANGE;
|
|
|
|
}
|
|
|
|
if (exc.err != 0)
|
|
|
|
errno = exc.err;
|
|
|
|
return (float)exc.retval;
|
|
|
|
}
|
|
|
|
return z;
|
|
|
|
#endif
|
|
|
|
}
|
|
|
|
|
|
|
|
#ifdef _DOUBLE_IS_32BITS
|
|
|
|
|
|
|
|
#ifdef __STDC__
|
|
|
|
double pow(double x, double y)
|
|
|
|
#else
|
|
|
|
double pow(x,y)
|
|
|
|
double x,y;
|
|
|
|
#endif
|
|
|
|
{
|
|
|
|
return (double) powf((float) x, (float) y);
|
|
|
|
}
|
|
|
|
|
|
|
|
#endif /* defined(_DOUBLE_IS_32BITS) */
|
New expf, exp2f, logf, log2f and powf implementations
Based on code from https://github.com/ARM-software/optimized-routines/
This patch adds a highly optimized generic implementation of expf,
exp2f, logf, log2f and powf. The new functions are not only
faster (6x for powf!), but are also smaller and more accurate.
In order to achieve this, the algorithm uses double precision
arithmetic for accuracy, avoids divisions and uses small table
lookups to minimize the polynomials. Special cases are handled
inline to avoid the unnecessary overhead of wrapper functions and
set errno to POSIX requirements.
The new functions are added under newlib/libm/common, but the old
implementations are kept (in newlib/libm/math) for non-IEEE or
pre-C99 systems. Targets can enable the new math code by defining
__OBSOLETE_MATH_DEFAULT to 0 in newlib/libc/include/machine/ieeefp.h,
users can override the default by defining __OBSOLETE_MATH.
Currently the new code is enabled for AArch64 and AArch32 with VFP.
Targets with a single precision FPU may still prefer the old
implementation.
libm.a size changes:
arm: -1692
arm/thumb/v7-a/nofp: -878
arm/thumb/v7-a+fp/hard: -864
arm/thumb/v7-a+fp/softfp: -908
aarch64: -1476
2017-05-25 23:41:38 +08:00
|
|
|
#endif /* __OBSOLETE_MATH */
|