Divide and conquer binomial

[Niels Möller]

Let me summarize my understanding of the binomial, and see if anybody
has any additions. I'm aware of a couple of different algorithms:

1. When the result is small (single limb, possibly with extensions to a
   few limbs):

   Tabulate the odd part of n! mod B and the corresponding mod B
   inverse. Then one can compute the odd part of bin(n, k) using n! /
   (k! (n-k)!) and these tables. To deal with the powers of two, use the
   cute formula

     ctz(n!) = n - popcount(n)

   (ctz is count trailing zeros, aka 2-adic valuation, often denoted
   v_2).

2. When k is small and n is large, just multiply together n (n-1) ...
   (n-k+1), preferably in an order using roughly balanced
   multiplications, and divide by k! using tabulated values and
   inverses.

3. When k and n are of moderate size, so that it is practical to iterate
   over all primes <= n, count multiplicity of each prime in dividend
   and divisor. I think Marco has implemented this.

4. When k is of moderate size, and n is too large for iterating over the
   primes, but n and k still small enough so that bin(n, k) fits in
   available storage, it would make sense to use some divide and conquer
   algorithm. For example, bin(2^80, 2^20), which is approximately
   2^{64427324.069447006302169138472946275466} (using pari/gp with
   default precision and the lngamma function) and needs 61 Mb of
   storage.

   First, one question is what complexity to expect. The result is of
   size N = O(k log n), so I don't think one can't hope for anything
   better than O(M(k log n)). I think one should be happy if one get
   within a log N factor of that.
   
   One possibly useful identity (assume for simplicity that k is
   even) is:

     bin(n, k) = bin(n, k/2) bin(n-k/2, k/2) / bin(k, k/2)

   Here, we halve the size of k, getting two recursive calls with
   essentially the same n (which means that we will unfortunately not
   get to use (3) as a base case, only (2). And one recursive call for
   the adjustment factor bin(k, k/2).

   I think Marco has tried this or something similar.

   * If the complexity is supposed to be close to linear in k, it seems
     right to divide into two pieces with size k/2. It should save work
     by both (i) getting a reasonably balanced multiplication tree, and
     (ii) divide out many common factors early.

   * What about bin(k, k/2)? I imagine that k is typically in range for
     method (3). Or are there any other shortcuts for this particular
     form?

   * The factor bin(k, k/2) is still pretty large, of size close to k
     (i.e., value on the order of 2^k, or 2^{2^20} in the example
     bin(2^80, 2^20)).

     If we can somehow divide out most of those factors earlier, that
     should save considerable work.

/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.