计算浮点数幂 (PHP/BCMath)

Question

我正在编写一个bcpow()的bug #10116特别恼人--它将$right_operand ($exp)转换为(本地PHP，而不是任意长度的)整数，所以当你尝试计算一个数字的平方根（或任何高于1的其他根）时，你总是得到1而不是正确的结果。

我开始寻找算法，以便能够计算一个数字的第n个根，我发现了这个答案，看起来非常可靠，我实际上扩展了这个公式使用WolframAlpha，我能够提高它的速度约5%，同时保持结果的准确性。

这是一个纯PHP实现，模仿了我的BCMath实现及其限制：

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)不产生整数。例如，如果$exp是0.123456，它的倒数将会是8.10005，pow()和_pow()的结果会有些不同（演示）：

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

如何使用“手动”指数计算实现相同级别的精度？

- Alix Axel

---

它的工作方式与广告完全一致。 _pow将小数部分“舍入”到最近的1/n。您可以递归地使其工作。因此，在计算_pow(2,0.125)之后，您可以计算_pow(2,0.125-123456)等等。 - Jeffrey Sax

啊，现在我明白了。那么bcmath没有exp和log吗？还是有其他原因导致a^b = exp(b*log(a))不可行？Jeffrey建议的递归当然可以工作，但如果需要许多1/k来表示指数，则其速度可能不令人满意。将指数写成有理数n/d并计算(a^n)^(1/d)是否可行，或者预计n和d太大了？也许值得研究的是用小分母（连分数展开）的有理数近似指数，并通过递归完成其余部分。 - Daniel Fischer

@JeffreySax：啊，我明白了...这很糟糕，但似乎还是不起作用（http://codepad.org/eI4ykyQU），或者我漏掉了什么？ - Alix Axel

@DanielFischer：感谢您回复我！=）嗯，bcmath API相当糟糕，除了*/+-之外，我们还有sqrt和残缺不全的pow：http://www.php.net/manual/en/ref.bc.php。我看到计算`(a^n)^(1/d)`的一个问题是`1/d`也可能是无理数。无论如何，我问这个问题主要是因为我很好奇——我怀疑我不需要在这么大的数字上使用无理指数。=） - Alix Axel

我认为我们可以安全地忽略无理数。我们可以用有理数任意逼近它们。问题在于这种逼近的分子和分母可能非常大。你能指定要处理的输入类型以及结果的精度吗？你需要的位数越少，逼近中就可以使用更小的分子和分母。 - Daniel Fischer

@DanielFischer：实际上，我只需要能够计算根号，你的答案完美地解决了这个问题。=）这个问题主要是因为我的一些测试（http://ideone.com/KNztd）没有通过，我很好奇是否可以用简单的方法解决。既然不行，而且这主要是为了满足我的好奇心，我投票关闭，随意发表你的第一个评论作为答案，我会接受它。=）顺便说一句，谢谢你的时间。 - Alix Axel

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 是分母较小的有理数，那么这个问题可以通过以下方式轻松解决：

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的例子，使用前四个部分商得到逼近值10/81 = 0.123456790123...（注意，10/81 = 1/8 - 1/648，最短分解成纯分数的部分和是收敛的）。然后将余数分解为纯分数。

然而，一般来说，这种方法会导致计算大的n次根，如果期望结果的精度很高，则速度和内存都会很慢。

总的来说，实现exp和log并使用它们可能会更简单和更快。

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

- Daniel Fischer

---

非常好的详细回答！谢谢。 - Alix Axel