php - 计算浮点幂 (PHP/BCMath)

Question

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

我开始寻找允许我计算数字的 n 次方根的算法，我的 found this answer 看起来非常可靠，实际上我使用 WolframAlpha 的 expanded the formula 并且我能够将它的速度提高大约 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))