7

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

我开始寻找可以让我计算数字的第 n 根的算法，我发现这个答案看起来很可靠，我实际上使用 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;
}

上面的方法似乎很好用， 除非1 / fmod($exp, 1)没有产生 integer。例如，如果$expis 0.123456，它的倒数将会是并且8.10005和 的结果会有点不同（demo）：pow()_pow()

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

如何使用“手动”指数计算达到相同的精度水平？

4

1 回答 1

5

用于找到（正）数的第 n次根的a算法是用于找到零的牛顿算法

f(x) = x^n - a.

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

用不是整数形式0 < y < 1的指数计算幂更复杂。做模拟，解决y1/nn

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的连分数展开找到小分母的接近有理近似值，并使用前四个部分商产生近似值[注意，最短的部分和分解成纯分数是收敛的] - 然后将余数分解成纯分数。y1929/15625 = [0;8,9,1,192]10/81 = 0.123456790123...10/81 = 1/8 - 1/648

然而，一般来说，这种方法会导致计算 large 的 n次根，n如果最终结果的期望精度很高，这也是缓慢且内存密集型的。

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

x^y = exp(y*log(x))
于 2012-05-10T01:03:38.497 回答