OCaml 中的整数取幂

Posted 2023-02-21

【中文标题】OCaml 中的整数取幂

【英文标题】：Integer exponentiation in OCaml

【发布时间】：2013-06-01 18:55:12

【问题描述】：

OCaml 中是否有整数求幂函数？ ** 仅适用于花车。虽然它看起来大部分是准确的，但是否存在精度错误的可能性，例如 2. ** 3. = 8. 有时返回 false ？是否有用于整数幂运算的库函数？我可以自己编写，但会涉及到效率问题，而且如果还没有这样的功能，我会感到惊讶。

【参考方案1】：

不在标准库中。但是您可以轻松地自己编写一个（使用exponentiation by squaring 更快），或者重用提供此功能的扩展库。在Batteries 中是Int.pow。

下面是一个建议的实现：

let rec pow a = function
  | 0 -> 1
  | 1 -> a
  | n -> 
    let b = pow a (n / 2) in
    b * b * (if n mod 2 = 0 then 1 else a)

如果由于处理非常大的数字而存在溢出风险，则可能应该使用大整数库，例如 Zarith，它提供了各种求幂函数。

（您可能需要“模幂运算”，计算 (a^n) mod p；这可以通过在中间计算中应用 mod 来避免溢出，例如在上面的函数 pow 中。）

【讨论】：

很好的答案。不幸的是，我只能选择一个最佳答案：/。此外，我不相信这始终是实现整数求幂的最快方法。事实上，我认为沿着这些思路存在一个 Project Euler 问题（我还没有解决）。我真的认为应该将整数幂运算添加到标准库中。即使它不再有效（我不确定这是不是真的），它也是一件很常见的事情，并且必须从浮点数转换和反转换是很烦人的。当然，导入库并不难，但没有理由这不应该是标准的。

好吧，如果您对如何在一般情况下实现整数实现有更好的想法，请随时提出实现建议。

@user2258552 我不同意你的假设，即整数取幂是如此普遍。在实践中，您几乎总是使用小的固定指数，或者您需要按照 gasche 的建议进行任意精度算术。 TL;DR：不要再相信你需要对固定精度整数进行整数取幂，并意识到你需要一个任意精度的算术库。

【参考方案2】：

关于您问题的浮点部分：OCaml 调用底层系统的 pow() 函数。浮点求幂是一个难以实现的函数，但它只需要忠实（即精确到一个Unit in the Last Place）就可以使2. ** 3. = 8. 计算为true，因为8.0 是唯一的float在数学上正确结果的一个 ULP 内 8.

所有数学库都应该（*）忠实，因此您不必担心这个特定的示例。但是not all of them actually are，所以你的担心是对的。

---

如果您使用 63 位整数或更宽的整数，那么担心的更好理由是，参数或求幂的结果不能完全表示为 OCaml 浮点数（实际上 IEEE 754 双精度数不能表示 @ 987654329@ 或 2⁵³ + 1)。在这种情况下，浮点取幂不能很好地替代整数取幂，这不是因为特定实现的弱点，而是因为它的设计目的不是精确地表示那么大的整数。

---

(*) 关于这个主题的另一个有趣的参考资料是 William Kahan 的“A Logarithm Too Clever by Half”。

【讨论】：

在参数相同的情况下，浮点取幂是否与整数取幂一样快？另外，为了清楚起见，您是否声明浮点取幂对于所有整数 a,b 应该是正确的，其中 -2^30 ≤ a^b

@user2258552 关于速度：浮点取幂可能比写得很好的整数慢。关于“应该”的含义：忠实的初等函数在其定义域中精确到一个 ULP。所有的库都应该是忠实的，因为它是计算成本和准确性之间的合理折衷，这将使几乎每个人都感到高兴。由于制表者的困境，所有 libms 的 0.5 ULP 精度有点过高，但以合理的成本可以实现 1 ULP 的精度。 （但同样，pow() 是最难的基本函数之一）

关于速度：鉴于此，那么在标准库中不包含整数求幂函数真的没有什么意义......

【参考方案3】：

这是另一种使用平方求幂的实现（就像@gasche 提供的那样），但这个实现是尾递归

let is_even n = 
  n mod 2 = 0

(* https://en.wikipedia.org/wiki/Exponentiation_by_squaring *)
let pow base exponent =
  if exponent < 0 then invalid_arg "exponent can not be negative" else
  let rec aux accumulator base = function
    | 0 -> accumulator
    | 1 -> base * accumulator
    | e when is_even e -> aux accumulator (base * base) (e / 2)
    | e -> aux (base * accumulator) (base * base) ((e - 1) / 2) in
  aux 1 base exponent

【讨论】：

请注意，尾递归对于其输入中的对数函数无关紧要。你怎么能炸掉堆栈？当然，如果尾递归给出了一个不同的观点来揭示代码的一些有趣之处，或者让它更容易阅读，那么它仍然很有趣。

@gasche 你是对的。此代码对 63 或 31 位整数没有意义。这样的算法对于任意精度的数字都是有意义的。