在 C 中计算小数指数

Question

我正在尝试评估 a^n，其中 a 和 n 是有理数。 我不想使用任何预定义的函数，例如 sqrt() 或 pow()

所以我尝试使用牛顿法来得到一个近似解：

3^0.2 = 3^(1/5) ，所以如果 x = 3^0.2，x^5 = 3。

可能是解决这个问题的最佳方法（没有计算器但仍然 使用基本算术运算）就是使用“牛顿法”。 牛顿求解方程 f(x)= 0 的方法是建立一个 通过将 x0 作为一些初始“猜测”定义的数字序列 xn 然后 xn+1= xn- f(xn/f '(xn) 其中 f '(x) 是 f 的导数。

发表于physicsforums

该方法的问题在于，如果我想计算5.2^0.33333，我需要找到这个方程x^10000 - 5.2^33333 = 0 的根。我最终得到了巨大的数字，并且大部分时间都得到inf 和nan 错误。

有人可以就如何解决这个问题给我建议吗？或者，有人可以提供另一种算法来计算 a^n 吗？

Answer 1

看来你的任务是计算

⎛ xN ⎞(aN / aD)
⎜⎼⎼⎼⎼⎟           where xN,xD,aN,aD ∈ ℤ,  xD,aD ≠ 0
⎝ xD ⎠

仅使用乘法、除法、加法和减法，建议使用Newton's method 来实现。

我们试图求解的方程（对于y）是

             (aN / aD)
y = (xN / xD)            where y ∈ ℝ

牛顿法找到一个函数的根。如果我们想用它来解决上面的问题，我们从左边减去右边，得到一个函数，其零给我们想要的 y：

                  (aN/aD)
f(y) = y - (xN/xD)        = 0

帮助不大。我想这就是你得到的？这里的重点是暂时不要形成该函数，因为我们没有办法计算有理数的有理幂！

首先，让我们确定 aD 和 xD 都是正数。如果 aD 是负数，我们可以简单地通过否定 aN 和 aD 来做到这一点（所以 aN/aD 不变），如果 xD 为负数，则同时否定 xN 和 xD。请记住，根据定义，xD 或 aD 都不是零。然后，我们可以简单地将两边都提高到 aD 次方：

 aD            aN     aN     aN
y   = (xN / xD)   = xN   / xD

我们甚至可以通过将两边都乘以最后一项来消除除法：

 aD     aN     aN
y   × xD   = xN

现在，这看起来很有希望！我们从中得到的函数是

        aD   aN     aN
f(y) = y   xD   - xN

牛顿法也需要导数，这显然是

f(y)            aD   aN
⎼⎼⎼⎼ = df(y) = y   xD   y / aD
 dy

牛顿的方法本身就依赖于迭代

             f(y)
y    = y  - ⎼⎼⎼⎼⎼⎼
 i+1    i    df(y)

如果你算算，你会发现迭代只是

                                 aD
                y[i]      y[i] xN
y[i+1] = y[i] - ⎼⎼⎼⎼ + ⎼⎼⎼⎼⎼⎼⎼⎼⎼⎼⎼⎼⎼⎼
                 aD           aD   aN
                       aD y[i]   xD

您不需要将所有 y 值保存在内存中；记住最后一个就足够了，当它们的差异足够小时就停止迭代。

上面还有求幂，但现在它们只是整数求幂，即

  aD
xN   = xN × xN × .. × xN
       ╰───────┬───────╯
              aD times

您可以非常简单地做到这一点，例如只需将参数乘以所需的次数，例如在 C 语言中，

double ipow(const double base, const int exponent)
{
    double result = 1.0;
    int    i;
    for (i = 0; i < exponent; i++)
        result *= base;
    return result;
}

有更有效的方法来做integer exponentiation，但是上面的函数应该是完全可以接受的。

最后一个问题是选择初始的y，这样你就可以收敛。您不能使用 0，因为（的幂）y 用作除法中的分母；你会得到零误差除法。就个人而言，我会检查结果是否应该是正数或负数，以及小于或大于 1 的数量级；两个规则来选择一个安全的初始y。

问题？

Answer 2

您可以使用generalized binomial theorem。替换 y=1 和 x=a-1。您可能希望根据所需的精度在足够的项之后截断无限级数。为了能够将术语数量与准确性联系起来，您需要确保x^r 术语的绝对值正在减少。因此，根据a 和n 的值，您应该应用公式来计算a^n 和a^(-n) 之一并使用它来获得您想要的结果。

Answer 3

整数次幂的解法是：

int poweri (int x, unsigned int y)
{
    int temp;
    if (y == 0)
        return 1;

    temp = poweri (x, y / 2);
    if ((y % 2) == 0)
        return temp * temp;
    else
        return x * temp * temp;
}

但是，平方根并不能提供干净的封闭解决方案。在wikipedia-square root 和Wolfram Mathworks Square Root Algorithms 可以找到一些很好的背景知识。两者都提供了几种可以满足您需求的方法，您只需要选择一种适合您的目的。

稍作修改，来自维基百科的这个例程（修改为返回平方根并提高准确性）返回一个令人惊讶的准确平方根。是的，联合的使用会有一些抱怨，而且它只在整数和浮点存储等效的情况下才有效，但是如果你在破解自己的平方根，这是相对有效的：

float sqrt_f (float x)
{
        float xhalf = 0.5f*x;
        union
        {
            float x;
            int i;
        } u;
        u.x = x;
        u.i = 0x5f3759df - (u.i >> 1);
        /* The next line can be repeated any number of times to increase accuracy */
        // u.x = u.x * (1.5f - xhalf * u.x * u.x);
        int i = 10;
        while (i--)
            u.x *= 1.5f - xhalf * u.x * u.x;

        return 1.0f / u.x;
}