Fast mpz_bin_ui

[Niels Möller]

Fredrik Johansson <fredrik.johansson at gmail.com> writes:

> The current mpz_bin_ui is not asymptotically fast. I propose replacing it
> with a call to mpz_bin_uiui when n fits one word (as already suggested in a
> code comment) and a plain divide and conquer product otherwise.

If I get this right, the speedup comes from a more balanced
multiplication tree for n (n-1)...(n-k+1), right?

I wonder if one can do something more clever. mpz_bin_ui_ui seems to
count the multiplicities of primefactors, somewhat related to the
factorial code, I think.

It would be nice to avoid having to compute k! and divide it out at the
end.

For a start, consider the powers of two, and for concreteness, say n =
1000, k = 15. Let !! denote semifactorial, k (k-2) (k-4).... 

We can compute k! = 15!, taking out powers of two, as

  15! = 15!! 14!! = 2^7 15!! 7! =  2^7 15!! 7!! 6!!
      = 2^10 15!! 7!! 3!
      = 2^11 3^3 (5*7)^2 (15*13*11*9)

(Further sieving could take other primes into account too, the current
factorial code does that, I think)

To do the same splitting for de numerator n(n-1)...(n-k+1), we have

  1000*999*...*986 = (1000*998*...*986) (999*997*...987)
    = 2^7 (500*499*...*493) (999*997*...987)
    = 2^7 (500*498*496*494) (499*497*495*493) (999*997*...987)
    = 2^11 (250*249*248*247) (499*497*495*493) (999*997*...987)
    = 2^11 (250*248) (249*247) (499*497*495*493) (999*997*...987)
    = 2^13 (125*124) (249*247) (499*497*495*493) (999*997*...987)

We could go on to take out factors of two of 124, but above is enough to
cancel out the power of two in 15!. Unlike the factorial case, we don't
get overlapping ranges and no nice squares or powers of products

 So we would get 1000 over 15 as

    2^2 (125*124) * (249*247) * (499*497*495*493) * (999*997*...987)
  -------------------------------------------------------------------
                   3^3 (5*7)^2 (15*13*11*9)

We could do sieving to take out the other factors of k!, but it's not
clear to me that would be an improvement. Assuming n is a bignum,
computing k! and dividing by it will be a fairly small operation, and
taking of the factors early won't save much work in computing the
numerator, as far as I see.

> Related to this, it would be useful to have a public function that computes
> rising factorials (say mpz_rising_ui). Then mpz_bin_ui could be implemented
> very simply by calling either mpz_bin_uiui or
> mpz_rising_ui+mpz_fac_ui+mpz_divexact.

There are a couple of building blocks in oddfac.c and prodlimbs.c, but
all for single-limb inputs, it seems.

If we want to do something reasonably simple, it would make sense to use
mpz_oddfac (which computes the odd part of the factorial), and mpz_fall_fac or
mpz_odd_fall_fac for computing the falling factorial or the odd part
thereof.

I'm not familiar with the sophisticated things done for mpz_bin_uiui and
mpz_fac, but maybe parts of that could be applicable also to the falling
factorials of a bignum?

It might work to first take out powers of two, and then rewrite
multiplies as squares like

  n (n-2)           = (n-1)^2 - 1
  n (n-2)(n-4)(n-6) = [(n-3)^2 - 5]^2 - 25

Regards,
/Niels

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