在 C 中计算分数指数

Question

我正在尝试评估 a^n，其中 a 和 n 是有理数。
我不想使用任何预定义的函数，例如 sqrt()或 pow()
所以我试图使用牛顿法来使用这种方法得到一个近似解:

3^0.2 = 3^(1/5) , so if x = 3^0.2, x^5 = 3.

Probably the best way to solve that (without a calculator but still using the basic arithmetic operations) is to use "Newton's method". Newton's method for solving the equation f(x)= 0 is to set up a sequence of numbers xn defined by taking x0 as some initial "guess" and then xn+1= xn- f(xn/f '(xn) where f '(x) is the derivative of 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，因为 (a power of) y 在除法中用作分母；你会得到零误差的除法。就我个人而言，我会检查结果应该是正数还是负数，并且数量级小于或大于 1；两个规则来选择一个安全的初始 y。

问题？