alternative power function

How about using fabs?

Using fabs won't cut it (try (a,b) = (-1,1)). What walkingcarcas wants is a sign-preserving power function.

I'm pretty sure there isn't a name for this. At least not one that I've encountered (as a maths graduate). Mathematicians tend to rely on even powers being positive-definite (over the reals), so p(a,b)=sgn(a)|a^b| is about as concise as it gets (although personally I'd prefer p(a,b)=a|a^b|/|a| as it extends meaningfully to complex numbers :)).

Regards
Admiral

Quote:

Original post by LessBread How about using fabs? *** Source Snippet Removed ***

Well you could write it as :

double alt2(double a, double b)
{
return (a<0 ? -1 : 1) * pow(fabs(a),b);
}

but this is actually a little bit slower it seems on my microsoft C++ 2003 compiler.

Sorry for the hijack, but since we're on the topic of pow(), can anyone point me to a fast integer/fixed-point approximation of this function? Any tips & pointers would be greatly appreciated. (I'm working on a Minkowski implementation of the Worley cell noise... the pow() function is quite taxing.)

I'm curently reading this article, which is somewhat helpful.
http://www.gamasutra.com/features/20020801/beaudoin_01.htm

Tachikoma, do you have any constraints on the domain of a or b? The real performance boosts in function approximation come from identifying limitations and exploiting them. If you need a good approximation when a and b wander over the whole real line, you'll have trouble finding it.

Regards
Admiral

If the number raised to a power doesn't change often but the exponent does, precalculate the natural log of the number reducing your pow to exponential functions.

If y = x^n
log y = log(x^n)
=> log y = n log x
=> y = exp( n log x )

so you can precalc log x and just do exp( nK ) where K = log x.

I don't know if this will help your case though.

I doubt this would be any faster, but if your exponents are always naturals, you can use a fast exponentiation algorithm. That's something like this:

double pow_int(double x, unsigned n){  if(n<=1)    return n==0?1.0:x;  if(n%2==0)    return pow_int(x*x,n/2);  return pow_int(x,n-1)*x;}

Well, non-recursive versions of this idea are easy to write and are probably faster.

You can easily adapt that function to work with fixed-point numbers if you need to, while I don't see an obvious way of doing it with non-natural exponents.

Sorry for the late response, but I've actually "lost" this thread and took me a while to find it on these massive forums. :S

Quote:

Original post by TheAdmiral Tachikoma, do you have any constraints on the domain of a or b? The real performance boosts in function approximation come from identifying limitations and exploiting them. If you need a good approximation when a and b wander over the whole real line, you'll have trouble finding it. Regards Admiral

The constraints are:
a is a 24.8 fixed-point number, always in the positive range of (0x00000000 <= a <= 0x000001FF).

b is also positive, typically anywere between (0.0f < b <= 20.0f). The range of greatest interest is between 0.0f and 6.0f, so I'd like to keep b in fractional increments in that range; anything beyond 6.0f can be integer increments. I'm quite flexible with b's constraints, so I'm willing to change them in favour of speed. At the moment b is a float, but I'd like to re-implement it as a 24.8 fixed-point number if possible.

The basic equation for the Minkowski distance is as follows:

Distance = (dU^b + dV^b)^(1/b)

...where dU, dV is the absolute difference between the cell's centroid and the current rendering point. b is an arbitrary exponent value (which is in the range specified above), but it remains constant while the texture is generated. In code it looks like this:

//Some arbitrary exponent b and it's reciprocal. float b = 2.4f;float ib = 1.0f / b;//U, V = cell coord //X, Y = render point coord //dU, dV = 24.8 fixed point difference, max range 0x0 ~ 0x1FFdU = labs(U - X);dV = labs(V - Y);//Minkowski distance function, result should be in 22.10 fixed point int Distance = (int)pow(pow((float)dU, b) + pow((float)dV, b), ib);

The fist two pow() functions that accept dU, dV and b is not a big deal, because I can use a look-up table with 512 precomputed values. (This table can be generated just before rendering the cell noise texture.)

My main issue is the final pow() function that uses 1/b as the exponent. A look-up table is more tricky here, because the base (i.e.. the result of the previous pow() functions)is non-linear and it's more difficult to use as a table index.

It would be nice if I could devise a single, compact look-up table for the Minkowski operation, by using dU and dV as a magical hash value for the table index. So instead of using a massive 512x512 2D table, i.e..
Distance = Minkowski_Table[dU][dV];
I prefer something more compact, such as...
Distance = Minkowski_Table[hash(dU, dV)];

Perhaps a partial fast analytical method combined with tables would be more realistic. I'm still investigating.

a bit off topic but...

it seems that for a faster power function for a floating point type you could make long* buffer pointing to the float, copy the value stored at that long* into a long, then copy the new exponent into the exponant field, then cast that back to a float and return it.

never tried that, just a random thought.

[edit] maybe it isn't a bit off topic since the above would preserve the sign of the number.

Your distance is of the form f(x,y) = (x^b+y^b)^{1/b}. Usually, b >= 1 is the constraint, not b > 0. When b < 1, you really do not even have a metric. From here on, I assume b >= 1. Your x and y values are in the interval [0,M] for some positive M. Define u = x/(2*M) and v = y/(2*M), so f(x,y) = 2*M*(u^b+v^b)^{1/b} with u and v in the interval [0,1/2]. The problem is efficiently computing (u^b+v^b)^{1/b}. Define z = u^b+v^b. The constraints on u and v and the condition b >= 1 guarantee that z <= 1.

The problem reduces to efficiently computing g(z) = z^{1/b}, where 0 <= z <= 1. Notice that g(z) = exp((1/b)*log(z)) = exp(-w), where w = -(1/b)*log(z) >= 0.

Now to compute h(w) = exp(-w) for w >= 0. Pulling these out of the air:

This approximation gives at most approximately 0.0002 absolute error:
exp(-w) = 1/(1+a1*w+a2*w^2+a3*w^3)^4
a1 = 0.2507213
a2 = 0.0292732
a3 = 0.0038278

This approximation gives at most approximately 0.00002 absolute error:
exp(-w) = 1/(1+a1*w+a2*w^2+a3*w^3+a4*w^4)^4
a1 = 0.24991035
a2 = 0.03158565
a3 = 0.00227723
a4 = 0.00026695

This approximation gives at most approximately 0.000002 absolute error:
exp(-w) = 1/(1+a1*w+a2*w^2+a3*w^3+a4*w^4+a5*w^5)^4
a1 = 0.250010936
a2 = 0.031198056
a3 = 0.002673255
a4 = 0.000127992
a5 = 0.000014876


This approximation gives at most approximately 0.0000002 absolute error:
exp(-w) = 1/(1+a1*w+a2*w^2+a3*w^3+a4*w^4+a5*w^5+a6*w^6)^4
a1 = 0.2499986842
a2 = 0.0312575832
a3 = 0.0025913712
a4 = 0.0001715620
a5 = 0.0000054302
a6 = 0.0000006906

Wasting Time, thanks for the info! Your suggestion looks very promising. :) I'll report back when I have a play around with it.

And to add something walkingcarcass' initial question, here is another way, without using branching nor using an extra multiplication:

//i=3 for 16-bit int, i=1 for 32-bit int, i=0 for 64-bit int. // You're in trouble if this is a 128-bit machine!#define i (sizeof(double)/sizeof(int) - 1)#define signmask (1 << (sizeof(int)*8 - 1))double alt_pow(double a, double b){double c = pow(a, b);((int*)&c) |= ( ((int*)&a) & signmask);return c;}

The trick here is to copy the sign bit in a to c. If you're on a 32-bit machine, this requires two 32-bit integer operations on the most significant 32-bit word of the 64-bit double variable.

The catch is, you must make sure the machine is big-endian, and that sizeof(double) >= sizeof(int)... so be careful.