By using complex math it is possible ... The problem is that the pow is using ln operation which is on complex domain multi-valued which means it has infinite number of valid result (as imaginary part of result is added with k*2*PI where k is any integer). So when the used result is not added with correct k then the result will be complex and not the one you want.
However for rooting (1.0/exponent -> integer) the k can be obtained directly as this:
int k=int(floor((1.0/exponent)+0.5))/2;
Which I discovered just now empirically. When I rewrite my complex math cpow to use this I got this C++ code:
//---------------------------------------------------------------------------
vec2 cadd(vec2 a,vec2 b) // a+b
{
return a+b;
}
vec2 csub(vec2 a,vec2 b) // a-b
{
return a-b;
}
vec2 cmul(vec2 a,vec2 b) // a*b
{
return vec2((a.x*b.x)-(a.y*b.y),(a.x*b.y)+(a.y*b.x));
}
vec2 cdiv(vec2 a,vec2 b) // a/b
{
float an=atan2(-a.y,-a.x)-atan2(-b.y,-b.x);
float r=length(a)/length(b);
return r*vec2(cos(an),sin(an));
}
vec2 csqr(vec2 a) // a^2
{
return cmul(a,a);
}
vec2 cexp(vec2 a) // e^a
{
// e^(x+y*i)= e^x * e^(y*i) = e^x * ( cos(y) + i*sin(y) )
return exp(a.x)*vec2(cos(a.y),sin(a.y));
}
vec2 cln(vec2 a) // ln(a) + i*pi2*k where k={ ...-1,0,+1,... }
{
return vec2(log(length(a)),atan2(a.y,a.x));
}
vec2 cpow(vec2 a,vec2 b) // a^b
{
return cexp(cmul(cln(a),b));
}
vec2 ctet(vec2 a,int b) // a^^b
{
vec2 c=vec2(1.0,0.0);
for (;b>0;b--) c=cpow(a,c);
return c;
}
//-------------------------------------------------------------------------
vec2 cln(vec2 a,int k) // ln(a) + i*pi2*k where k={ ...-1,0,+1,... }
{
return vec2(log(length(a)),atan2(a.y,a.x)+float(k+k)*M_PI);
}
float mypow(float a,float b) // a^b
{
if (b<0.0) return 1.0/mypow(a,-b);
int k=0;
if ((a<0.0)&&(b<1.0)) // rooting with negative base
{
k=floor((1.0/b)+0.5);
k/=2;
}
return cexp(cmul(cln(vec2(a,0.0),k),vec2(b,0.0))).x;
}
//-------------------------------------------------------------------------
I am using GLSL like math vec2 you can easily rewrite to x,y components instead for example to this:
//-------------------------------------------------------------------------
float mypow(float a,float b) // a^b
{
if (b<0.0) return 1.0/mypow(a,-b);
int k=0;
if ((a<0.0)&&(b<1.0)) // rooting with negative base
{
k=floor((1.0/b)+0.5);
k/=2;
}
float x,y;
x=log(fabs(a));
y=atan2(0.0,a)+(float(k+k)*M_PI);
x*=b; y*=b;
// if (fabs(exp(x)*sin(y))>1e-6) throw domain error; // abs imaginary part is not zero
return exp(x)*cos(y); // real part of result
}
//-------------------------------------------------------------------------
Which returns real part of complex domain pow and also select the correct multi-valued cln sub-result. I also added domain test as a comment in the code in case you need/want to implement it.
Here result for your case:
mypow(-12.411202,0.200000) = -1.654866
Have tested more number and 1/odd number exponents and looks like it is working.
[Edit1] handling mixed exponents
Now if we have exponents in form a^(b0+1/b1) where b0,b1 are integers and a<0 we need to dissect the pow to 2 parts:
a^(b0+1/b1) = a^b0 * a^(1/b1)
so when added to above function:
//-------------------------------------------------------------------------
float mypow(float a,float b) // a^b
{
if (b<0.0) return 1.0/mypow(a,-b); // handle negative exponents
int k;
float x0,x,y,e,b0;
k=0; // for normal cases k=0
b0=floor(b); // integer part of exponent
b-=b0; // decimal part of exponent
// integer exponent (real domain power |a|^b0 )
x0=pow(fabs(a),b0);
if ((a<0.0)&&((int(b0)&1))==1) x0=-x0; // just add sign if odd exponent and negative base
// decimal exponent (complex domain rooting a^b )
if ((a<0.0)&&(b>0.0)) // rooting with negative base
{
k=floor((1.0/b)+0.5);
k&=0xFFFFFFFE;
}
x=b*(log(fabs(a))); e=exp(x);
y=b*(atan2(0.0,a)+(float(k)*M_PI));
x=e*cos(y);
// y=e*sin(y); if (fabs(y)>1e-6) throw domain error;
return x*x0; // full complex result is x0*(x+i*y)
}
//-------------------------------------------------------------------------
and sample output:
mypow(-12.41120243,3.200000) = 3163.76635742
mypow(3163.76635742,0.312500) = 12.41120243
