计算浮点幂 (PHP/BCMath)

Question

我正在为bcmath 扩展编写一个包装器，而关于bcpow() 的bug #10116 特别烦人——它将$right_operand ($exp) 转换为(本机 PHP，不是任意长度）整数，因此当您尝试计算一个数字的平方根（或任何其他高于 1 的根）时，您总是以 1 而不是正确的结果结束。

我开始寻找可以让我计算数字的 n 次根的算法，我 found this answer 看起来很可靠，我实际上 expanded the formula 使用 WolframAlpha 并且我能够将它的速度提高大约 5%，同时保持结果的准确性。

这是一个模仿我的 BCMath 实现及其局限性的纯 PHP 实现：

function _pow($n, $exp)
{
    $result = pow($n, intval($exp)); // bcmath casts $exp to (int)

    if (fmod($exp, 1) > 0) // does $exp have a fracional part higher than 0?
    {
        $exp = 1 / fmod($exp, 1); // convert the modulo into a root (2.5 -> 1 / 0.5 = 2)

        $x = 1;
        $y = (($n * _pow($x, 1 - $exp)) / $exp) - ($x / $exp) + $x;

        do
        {
            $x = $y;
            $y = (($n * _pow($x, 1 - $exp)) / $exp) - ($x / $exp) + $x;
        } while ($x > $y);

        return $result * $x; // 4^2.5 = 4^2 * 4^0.5 = 16 * 2 = 32
    }

    return $result;
}

上述seems to work great 除非1 / fmod($exp, 1) 不产生整数。例如，如果$exp 是0.123456，它的倒数将是8.10005，而pow() 和_pow() 的结果会有点不同（demo）：

• pow(2, 0.123456) = 1.0893412745953
• _pow(2, 0.123456) = 1.0905077326653
• _pow(2, 1 / 8) = _pow(2, 0.125) = 1.0905077326653

如何使用“手动”指数计算达到同样的准确度？

Answer 1

寻找（正）数a的第n^th根的算法是牛顿算法寻找零的

f(x) = x^n - a.

这仅涉及以自然数为指数的幂，因此很容易实现。

用指数0 < y < 1 计算幂，其中y 的形式不是1/n 的整数n 更复杂。做类比，解决

x^(1/y) - a == 0

将再次涉及计算具有非整数指数的幂，这正是我们要解决的问题。

如果y = n/d是小分母d的有理数，那么问题很容易通过计算解决

x^(n/d) = (x^n)^(1/d),

但是对于大多数理性的0 < y < 1，分子和分母都相当大，而中间的x^n 会很大，因此计算会占用大量内存并且需要（相对）较长的时间。 （以0.123456 = 1929/15625 的指数为例，还不错，但0.1234567 会比较费力。）

计算一般理性0 < y < 1 的幂的一种方法是写

y = 1/a ± 1/b ± 1/c ± ... ± 1/q

使用整数a < b < c < ... < q 并乘/除单个x^(1/k)。 （每一个有理的0 < y < 1都有这样的表示，而最短的这样的表示一般不会涉及很多术语，例如

1929/15625 = 1/8 - 1/648 - 1/1265625;

在分解中仅使用加法会导致具有更大分母的更长表示，例如

1929/15625 = 1/9 + 1/82 + 1/6678 + 1/46501020 + 1/2210396922562500,

这样会涉及更多的工作。）

通过混合这些方法可以进行一些改进，首先通过y 的连分数展开找到具有小分母的y 的近似有理逼近 - 例如指数1929/15625 = [0;8,9,1,192] 并使用前四个部分商产生近似10/81 = 0.123456790123... [注意10/81 = 1/8 - 1/648，最短分解成纯分数的部分和是收敛的] - 然后将余数分解成纯分数。

但是，一般来说，这种方法会导致计算大 n 的 n^th 个根，如果最终结果的所需精度很高，这也会很慢且占用大量内存。

总而言之，实现exp和log并使用可能更简单、更快

x^y = exp(y*log(x))