Integer exponentiation in OCaml

integer ocaml floating-accuracy exponentiation

39,007

Solution 1

Regarding the floating-point part of your question: OCaml calls the underlying system's pow() function. Floating-point exponentiation is a difficult function to implement, but it only needs to be faithful (that is, accurate to one Unit in the Last Place) to make 2. ** 3. = 8. evaluate to true, because 8.0 is the only float within one ULP of the mathematically correct result 8.

All math libraries should(*) be faithful, so you should not have to worry about this particular example. But not all of them actually are, so you are right to worry.

A better reason to worry would be, if you are using 63-bit integers or wider, that the arguments or the result of the exponentiation cannot be represented exactly as OCaml floats (actually IEEE 754 double-precision numbers that cannot represent 9_007_199_254_740_993 or 2⁵³ + 1). In this case, floating-point exponentiation is a bad substitute for integer exponentiation, not because of a weakness in a particular implementation, but because it is not designed to represent exactly integers that big.

(*) Another fun reference to read on this subject is “A Logarithm Too Clever by Half” by William Kahan.

Solution 2

Not in the standard library. But you can easily write one yourself (using exponentiation by squaring to be fast), or reuse an extended library that provides this. In Batteries, it is Int.pow.

Below is a proposed implementation:

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)

If there is a risk of overflow because you're manipulating very big numbers, you should probably use a big-integer library such as Zarith, which provides all sorts of exponentiation functions.

(You may need the "modular exponentiation", computing (a^n) mod p; this can be done in a way that avoids overflows by applying the mod in the intermediary computations, for example in the function pow above.)

Solution 3

Here's another implementation which uses exponentiation by squaring (like the one provided by @gasche), but this one is tail-recursive

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

39,007

Author by

user2258552

SOreadytohelp

Updated on October 07, 2020

Comments

user2258552 over 3 years

Is there a function for integer exponentiation in OCaml? ** is only for floats. Although it seems to be mostly accurate, isn't there a possibility of precision errors, something like 2. ** 3. = 8. returning false sometimes? Is there a library function for integer exponentiation? I could write my own, but efficiency concerns come into that, and also I'd be surprised if there isn't such a function already.
user2258552 almost 11 years

Great answer. Unfortunately, I could only pick one best answer :/. Also, I'm not convinced that this is always the fastest way to implement integer exponentiation. In fact, I think there's a Project Euler question (which I have yet to solve) along these lines. I really think integer exponentiation should be added to the standard library. Even if it is not any more efficient (which I'm not sure is true), it's a pretty common thing to need and having to convert and deconvert from floats is annoying. Of course importing a library isn't hard, but there's no reason this shouldn't be standard.
user2258552 almost 11 years

Is floating-point exponentiation as fast as integer exponentiation given that the arguments are the same? Also, just to be clear, is your statement that floating-point exponentiation should be faithful true for all integers a,b for which -2^30 ≤ a^b < 2^30, or just for my particular example of 2 3 and 8?
Pascal Cuoq almost 11 years

@user2258552 Regarding speed: floating-point exponentiation is probably slower than a well-written integer one. Regarding the meaning of “should”: a faithful elementary function is accurate to one ULP across its definition domain. All libms should be faithful because it is a reasonable compromise between computational cost and accuracy that would make pretty much everyone happy. Accuracy to 0.5 ULP is a bit too much to expect of all libms because of the table-maker's dilemma, but accuracy to 1 ULP is achievable at reasonable cost. (but, again, pow() is one of the hardest elementary functions)
user2258552 almost 11 years

Regarding speed: in light of that, then it really makes very little sense to not include an integer exponentiation function in the standard library...
gasche almost 11 years

Well, if you have better ideas of how to implement integer implementation in the general case, feel free to suggest an implementation.
The Show that never ends almost 8 years

@user2258552 I do not agree with your premise that integer exponentiation is so common. In practice, you almost always work with small fixed exponents, or you need an arbitrary precision arithmetic as suggested by gasche. TL;DR: Stop to believe you need integer exponentiation on fixed precision ints and realize you need an arbitrary precision arithmetic library.
gasche almost 8 years

Note that tail-recursivity does not matter for a function logarithmic in its input. How could you ever blow the stack? Of course, if tail-recursivity gives a different viewpoint that reveals something interesting about the code, or makes it easier to read, then it can still be interesting.
Vladimir Keleshev over 7 years

@gasche you're right. This code does not make sense for 63 or 31 bit integers. An algorithm like that would make sense for arbitrary-precision numbers.