Let me write some more, after pondering the problem for a few days. Regular powers n^k (mod m) can be efficiently computed based on repeated squaring and the basic rule that
n^{2k} = (n^2)^k = (n^2)^k (and without mod m, one would use the same algorithm, it's just that the size of the result grows linearly with k). First, let's note that there's unlikely to be any comparatively efficient way to compute falling factorial powers, for it there were, there would be an efficient way to compute n! (mod m), and that would be a fast factoring algorithm: if n is the product of two unknown primes, then the smallest factor is p = gcd (n, floor(sqrt n)! mod n ) The closest analog of a squaring rule for falling factorial powers I've been able to find is n^_{2k} = n^{_k} (n-k)^{_k} where n^{_k} denotes the falling factorial n (n-1) ... (n-k+1). And the right hand side doesn't look much like a square, I'd wish we could rewrite it to a formula like n^_{2k} = ((n-k/2)^{_k})^{_2} + some small number but I'm afraid we can't. So, back to binomial coefficients, n \over k = n^{_k} / k! and the case where n is larger than one limb, k a single limb, and k <= 2n. I guess k has to be reasonably small, for the result to fit in available memory (but I don't have any numbers handy). Then k! is going to be much smaller than n^{_k}, so that computing k! and doing the exact division should be a very small fraction of the total work. So for k!, I'd suggest just using mpz_oddfac_1 (and it would make sense to me to change it to return the exponent for the omitted power of two, which from the implementation of mpz_fac_ui seems to be k - popcount(k)). For n^_{k}, I would suggest a function to compute the odd part and return the exponent of the power of two, which would repeatedly split in even and odd numbers and compute products of the form n (n-2) (n-4) ... for odd n. Each of these products could be computed by grouping four terms at a time and use > n (n-2)(n-4)(n-6) = [(n-3)^2 - 5]^2 - 25 I like this trick, replacing three multiplies by two squares (but it's unclear how much it will save compared to the total running time). For the reasons given above, I wouldn't expect the trick to generalize much further. Then multiply all the resulting odd factors together in some roughly balanced tree. And then we could do something like void mpz_bin_ui (mpz_ptr r, mpz_srcptr n, unsigned long int k) { mpz_t k_fac; mp_bitcnt_t k_exp, n_exp; mpz_init (k_fac); k_exp = mpz_oddfac_1 (k_fac, k); n_exp = mpz_odd_falling_fac (r, n, k); ASSERT (k_exp <= n_exp); mpz_divexact (r, r, k_fac); mpz_mul_2exp (r, r, n_exp - k_exp); mpz_clear (k); } Would you like to give it a try? One could try some sieving and instead produce n^_{k} as a product of prime powers for all primes <= k and a residual, to be able to easily subtract the corresonding prime-power representation of k!. It would be interesting with some analysis of the size of n \over k after dividing out all prime factors smaller than k. If the prime powers of primes <= k makes up a large part, sieving could pay off, since that part can be computed efficiently. Will n \over k be more or less "smooth" than a random number of similar size? I'd guess slightly more smooth, since it's going to be huge compared to n, but still can't have any prime factors larger than n. But we can still have lots of prime factors in the range k < p <= n. Regards, /Niels -- Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677. Internet email is subject to wholesale government surveillance. _______________________________________________ gmp-devel mailing list gmp-devel@gmplib.org https://gmplib.org/mailman/listinfo/gmp-devel