通过使用复杂的数学是可能的......问题是pow 正在使用ln 操作，它是在复域多值上，这意味着它有无限数量的有效结果（作为结果的虚部与k*2*PI 相加，其中k 是任意整数）。因此，当使用的结果未添加正确的k 时，结果将是复杂的，而不是您想要的。
但是对于求根（1.0/exponent -> integer），k 可以这样直接获得：
int k=int(floor((1.0/exponent)+0.5))/2;
我刚才凭经验发现的。当我重写我的 complex math cpow 来使用它时，我得到了这个 C++ 代码：
//---------------------------------------------------------------------------
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;
}
//-------------------------------------------------------------------------
我正在使用像数学vec2 这样的 GLSL，您可以轻松地重写为 x,y 组件，例如：
//-------------------------------------------------------------------------
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
}
//-------------------------------------------------------------------------
返回复域pow 的实部，并选择正确的多值cln 子结果。我还在代码中添加了域测试作为注释，以防您需要/想要实现它。
您的案例的结果如下：
mypow(-12.411202,0.200000) = -1.654866
已经测试了更多的数字和1/odd number 指数，看起来它正在工作。
[Edit1] 处理混合指数
现在，如果我们有a^(b0+1/b1) 形式的指数，其中b0,b1 是整数，a<0 我们需要将战俘分解为两部分：
a^(b0+1/b1) = a^b0 * a^(1/b1)
所以当添加到上述函数时：
//-------------------------------------------------------------------------
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)
}
//-------------------------------------------------------------------------
和样本输出：
mypow(-12.41120243,3.200000) = 3163.76635742
mypow(3163.76635742,0.312500) = 12.41120243