Re: Compiler has no way to predict how pow() is implemented
Was poking around in my numerics code, and found the code from the Cephes math library below (attribution included) to illustrate the point I wanted to make - that pow(x,y) is exp(y log x). Note that FORTRAN has an exponentiation operator ** and that it works differently for small integers as exponents - so there would be alternate entry points.
/*\t\t\t\t\t\t\tpowf.c
*
*\tPower function
*
*
*
* SYNOPSIS:
*
* float x, y, z, powf();
*
* z = powf( x, y );
*
*
*
* DESCRIPTION:
*
* Computes x raised to the yth power. Analytically,
*
* x**y = exp( y log(x) ).
*
* Following Cody and Waite, this program uses a lookup table
* of 2**-i/16 and pseudo extended precision arithmetic to
* obtain an extra three bits of accuracy in both the logarithm
* and the exponential.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,10 100,000 1.4e-7 3.6e-8
* 1/10 < x < 10, x uniformly distributed.
* -10 < y < 10, y uniformly distributed.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* powf overflow x**y > MAXNUMF MAXNUMF
* powf underflow x**y < 1/MAXNUMF 0.0
* powf domain x<0 and y noninteger 0.0
*
*/
\x0c
/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include "mconf.h"
static char fname[] = {"powf"};
/* 2^(-i/16)
* The decimal values are rounded to 24-bit precision
*/
static float A[] = {
1.00000000000000000000E0,
9.57603275775909423828125E-1,
9.17004048824310302734375E-1,
8.78126084804534912109375E-1,
8.40896427631378173828125E-1,
8.05245161056518554687500E-1,
7.71105408668518066406250E-1,
7.38413095474243164062500E-1,
7.07106769084930419921875E-1,
6.77127778530120849609375E-1,
6.48419797420501708984375E-1,
6.20928883552551269531250E-1,
5.94603538513183593750000E-1,
5.69394290447235107421875E-1,
5.45253872871398925781250E-1,
5.22136867046356201171875E-1,
5.00000000000000000000E-1
};
/* continuation, for even i only
* 2^(i/16) = A[i] + B[i/2]
*/
static float B[] = {
0.00000000000000000000E0,
-5.61963907099083340520586E-9,
-1.23776636307969995237668E-8,
4.03545234539989593104537E-9,
1.21016171044789693621048E-8,
-2.00949968760174979411038E-8,
1.89881769396087499852802E-8,
-6.53877009617774467211965E-9,
0.00000000000000000000E0
};
/* 1 / A[i]
* The decimal values are full precision
*/
static float Ainv[] = {
1.00000000000000000000000E0,
1.04427378242741384032197E0,
1.09050773266525765920701E0,
1.13878863475669165370383E0,
1.18920711500272106671750E0,
1.24185781207348404859368E0,
1.29683955465100966593375E0,
1.35425554693689272829801E0,
1.41421356237309504880169E0,
1.47682614593949931138691E0,
1.54221082540794082361229E0,
1.61049033194925430817952E0,
1.68179283050742908606225E0,
1.75625216037329948311216E0,
1.83400808640934246348708E0,
1.91520656139714729387261E0,
2.00000000000000000000000E0
};
#ifdef DEC
#define MEXP 2032.0
#define MNEXP -2032.0
#else
#define MEXP 2048.0
#define MNEXP -2400.0
#endif
/* log2(e) - 1 */
#define LOG2EA 0.44269504088896340736F
extern float MAXNUMF;
#define F W
#define Fa Wa
#define Fb Wb
#define G W
#define Ga Wa
#define Gb u
#define H W
#define Ha Wb
#define Hb Wb
#ifdef ANSIC
float floorf( float );
float frexpf( float, int *);
float ldexpf( float, int );
float powif( float, int );
#else
float floorf(), frexpf(), ldexpf(), powif();
#endif
/* Find a multiple of 1/16 that is within 1/16 of x. */
#define reduc(x) 0.0625 * floorf( 16 * (x) )
#ifdef ANSIC
float powf( float x, float y )
#else
float powf( x, y )
float x, y;
#endif
{
float u, w, z, W, Wa, Wb, ya, yb;
/* float F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
int e, i, nflg;
nflg = 0;\t/* flag = 1 if x<0 raised to integer power */
w = floorf(y);
if( w < 0 )
\tz = -w;
else
\tz = w;
if( (w == y) && (z < 32768.0) )
\t{
\ti = w;
\tw = powif( x, i );
\treturn( w );
\t}
if( x <= 0.0F )
\t{
\tif( x == 0.0 )
\t\t{
\t\tif( y == 0.0 )
\t\t\treturn( 1.0 ); /* 0**0 */
\t\telse
\t\t\treturn( 0.0 ); /* 0**y */
\t\t}
\telse
\t\t{
\t\tif( w != y )
\t\t\t{ /* noninteger power of negative number */
\t\t\tmtherr( fname, DOMAIN );
\t\t\treturn(0.0);
\t\t\t}
\t\tnflg = 1;
\t\tif( x < 0 )
\t\t\tx = -x;
\t\t}
\t}
/* separate significand from exponent */
x = frexpf( x, &e );
/* find significand in antilog table A[] */
i = 1;
if( x <= A[9] )
\ti = 9;
if( x <= A[i+4] )
\ti += 4;
if( x <= A[i+2] )
\ti += 2;
if( x >= A[1] )
\ti = -1;
i += 1;
/* Find (x - A[i])/A[i]
* in order to compute log(x/A[i]):
*
* log(x) = log( a x/a ) = log(a) + log(x/a)
*
* log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a
*/
x -= A[i];
x -= B[ i >> 1 ];
x *= Ainv[i];
/* rational approximation for log(1+v):
*
* log(1+v) = v - 0.5 v^2 + v^3 P(v)
* Theoretical relative error of the approximation is 3.5e-11
* on the interval 2^(1/16) - 1 > v > 2^(-1/16) - 1
*/
z = x*x;
w = (((-0.1663883081054895 * x
+ 0.2003770364206271) * x
- 0.2500006373383951) * x
+ 0.3333331095506474) * x * z;
w -= 0.5 * z;
/* Convert to base 2 logarithm:
* multiply by log2(e)
*/
w = w + LOG2EA * w;
/* Note x was not yet added in
* to above rational approximation,
* so do it now, while multiplying
* by log2(e).
*/
z = w + LOG2EA * x;
z = z + x;
/* Compute exponent term of the base 2 logarithm. */
w = -i;
w *= 0.0625; /* divide by 16 */
w += e;
/* Now base 2 log of x is w + z. */
/* Multiply base 2 log by y, in extended precision. */
/* separate y into large part ya
* and small part yb less than 1/16
*/
ya = reduc(y);
yb = y - ya;
F = z * y + w * yb;
Fa = reduc(F);
Fb = F - Fa;
G = Fa + w * ya;
Ga = reduc(G);
Gb = G - Ga;
H = Fb + Gb;
Ha = reduc(H);
w = 16 * (Ga + Ha);
/* Test the power of 2 for overflow */
if( w > MEXP )
\t{
\tmtherr( fname, OVERFLOW );
\treturn( MAXNUMF );
\t}
if( w < MNEXP )
\t{
\tmtherr( fname, UNDERFLOW );
\treturn( 0.0 );
\t}
e = w;
Hb = H - Ha;
if( Hb > 0.0 )
\t{
\te += 1;
\tHb -= 0.0625;
\t}
/* Now the product y * log2(x) = Hb + e/16.0.
*
* Compute base 2 exponential of Hb,
* where -0.0625 <= Hb <= 0.
* Theoretical relative error of the approximation is 2.8e-12.
*/
/* z = 2**Hb - 1 */
z = ((( 9.416993633606397E-003 * Hb
+ 5.549356188719141E-002) * Hb
+ 2.402262883964191E-001) * Hb
+ 6.931471791490764E-001) * Hb;
/* Express e/16 as an integer plus a negative number of 16ths.
* Find lookup table entry for the fractional power of 2.
*/
if( e < 0 )
\ti = -( -e >> 4 );
else
\ti = (e >> 4) + 1;
e = (i << 4) - e;
w = A[e];
z = w + w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
z = ldexpf( z, i ); /* multiply by integer power of 2 */
if( nflg )
\t{
/* For negative x,
* find out if the integer exponent
* is odd or even.
*/
\tw = 2 * floorf( (float) 0.5 * w );
\tif( w != y )
\t\tz = -z; /* odd exponent */
\t}
return( z );
}
-drl
