使用 pow 函数时如何避免 -nan(ind) 错误是否可以使用带有非整数指数的负底数的 pow

问题描述

我在使用 pow 函数 -nan(ind) 打印到屏幕时出错。想知道是否有一种方法可以将 pow 用于具有负底数和非整数指数的数字。

目前的 pow 函数是 pow(-12.4112021858,0.2)。并给我 -nan(ind) 错误。 如果我将基数更改为 12.41，它似乎计算得很好。

编辑 - 我将指数设置为与 0.2 不同的数字

int main() {



    double a = 1.83;
    double v = 1.25;
    double r = 0;

    double sum = 0;
 
    double exponent = 0.2;

    double result = 1;
    while (true)
    {
    printf("Enter Radius \n");
    scanf_s("%lf",&r);
    sum = 1 - r/ a;
    printf("%lf\n",sum);
    sum = sum * v;
    printf("%lf\n",sum);
    sum=  pow(sum,exponent);
    printf("%lf\n",sum);

}
}

解决方法

想知道是否有办法将 pow 用于具有负底数和非整数指数的数字。

首先，检查文档。 C 2018 7.12.7.4 指定了 powf、pow 和 powl：

pow 函数计算 x 的 y 次幂。如果 x 是有限且负数且 y 是有限且不是整数值，则会发生域错误......

您的 x，大约为 −12.4112021858，是有限的和负数，而您的 y，大约为 0.2，是有限的而不是整数值。所以出现域错误。

这意味着您不能期望得到结果，除非您使用的 pow 专门为超出 C 标准要求的其他情况提供支持，即使 .2 在 double 中完全表示。

（当出现域错误时，返回一个实现定义的结果。这可能是一个 NaN，一个有效的数学结果，例如 -2 表示 pow(-32,.2) 如果使用基于十进制的浮点数或其他内容。实现也可能通过 errno 或浮点异常报告错误。有关详细信息，请参阅 C 2018 7.12.1。）

其次，.2 在大多数 C 实现的 double 格式中无法表示。 C 实现通常使用 IEEE-754 binary64 格式。在这种格式中，最接近 .2 的可表示值是 0.200000000000000011102230246251565404236316680908203125。对于源代码pow(-12.4112021858,.2)，首先将数字转换为double，然后使用参数-12.41120218580000056363132898695766925811767578125和参数调用pow 0.200000000000000011102230246251565404236316680908203125。因此，您不是在请求具有实数结果的操作。

如果您的 C 实现使用基于十进制的浮点数，则 .2 将是可表示的，并且 pow(-12.4112021858,.2) 返回 x 的第五个根是合理的，因为第五个根是负实数。 （这将是对 pow 标准规范的扩展，如上所述。）

如果您知道 y 应该是五分之一或有理数 p/q 其中 q 是奇数，如果 p 是偶数，您可以将所需结果计算为 pow(fabs(x),y)，如果 p 是奇数，则计算为 copysign(pow(fabs(x),y),x)。

建议使用 cpow 的评论之一，但这不会产生您想要的结果。 cpow(-12.4112021858,.2) 将返回大约 1.3388 + .9727 i。 （复数幂“函数”是多值的，但定义了 cpow 以产生该结果。）

,

通过使用复杂的数学是可能的...问题是 pow 使用的是 ln 操作，它是在复杂域多值上，这意味着它具有无限数量的有效结果（如结果的虚部加上 k*2*PI，其中 k 是任何整数）。因此，当使用的结果没有添加正确的 k 时，结果将是复杂的，而不是您想要的。

但是对于生根（1.0/指数 -> 整数），可以直接获得 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,+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,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;
    }
//-------------------------------------------------------------------------

我正在使用 GLSL 之类的数学 vec2 您可以轻松地将其重写为 x,y 组件，例如：

//-------------------------------------------------------------------------
float mypow(float 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 我们需要将 pow 分解为 2 部分：

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