Compute (x+yCx)%mod

I want to compute, as quickly as possible,

$$\dfrac{(x+y)!}{x!y!} \mod m$$

with $x,y \le 10^6$ and a prime $m=10^9+7$ (called mod in my code).

My current approach, computing

$$\left(\mathtt{Prod}(y+1, y+x) \mod m\right) \cdot \left(\mathtt{InverseMod}(\mathtt{factMod}(x) \mod m)\right)$$

is too slow. Here is the implementation:

static inline ll powmod(unsigned a, unsigned b){
register ll x=1,y=a;
while(b){
if(b&1){
x*=y; if(x>=MOD)x%=MOD;
}
y*=y; if(y>=MOD)y%=MOD;
b>>=1;
}
return x;
}

static inline ll InverseMod(ll n){
return powmod(n,MOD-2);
}

static inline ll prodMod(ll minx,ll maxx){
for(unsigned i=minx+1; i<=maxx; i++){
minx*=i; if(minx>=MOD)minx%=MOD;
}
return minx;
}

static inline ll factMOD(unsigned n){
register ll ans0=1,ans1=1;register unsigned i,m;
if(n&1){
for(i=1,m=(n+1)>>1; i<m; i++){
ans0*=(i<<1); if(ans0>=MOD)ans0%=MOD;
ans1*=((i<<1)+1); if(ans1>=MOD)ans1%=MOD;
}
return ans0*ans1%MOD;
}else{
for(i=1; i<=n; i++){
ans0*=i; if(ans0>=MOD)ans0%=MOD;
}
return ans0;
}
}

static inline ll chemin(ll x, ll y){
return (prodMod(x+1,x+y) *
InverseMod(factMOD(y)))%MOD;
}

I've found and implemented Wilson's theorem. The factMODWilson function is the one to call to compute (n!) % MOD when MOD-n is little against n, but here it's not the case.

Does someone know an efficient approach when n=1e6 and MOD=1e9+7?

ll factMODWilson(ll n){ //n! % MOD efficient when MOD-n<n
ll res=1,i;
for(i=1; i<MOD-n; i++){
res*=i;
if(res>=MOD)res%=MOD;
}
res=InverseMod(res);
if(!(n&1))
res= -res +MOD;
}
return res%MOD;
}
• Using assembly is faster than using plain C++. – user3725053 Nov 26 '14 at 15:08
• Maybe you can generalize the prime powers approach: cs.stackexchange.com/questions/14456/… – Yves Daoust Nov 26 '14 at 15:15
• @user3725053 : In this case, using the right algorithm for the job is faster. Translating the code presented here into assembly by hand will not speed up its execution by a factor that's relevant to the problem. – Daniel Kamil Kozar Nov 26 '14 at 15:17
• @YvesDaoust interesting approach for factorials, but it doesn't seems simple to fit this approach to a factorial modular. – bilbo Nov 26 '14 at 16:33
• If x >= y, (x+y)!/(x!)(y)! = (x+1)(x+2)...(x+y)/y! If x < y, (x+y)!/(x!)(y)! = (y+1)(y+2)...(y+x)/x! – rcgldr Nov 26 '14 at 18:44

On my system, Windows 7 64 bit, Intel 2600K 3.4ghz, in 64 bit mode, this example is about 3 times faster. This code was compiled with Visual Studio Express 2013. Update - in prodMod and factMod I changed prod = 1ull; to prod = i++; for one less loop. 2nd update - made function names consistent, ProdMod constant good for <= 2e6, and FactMod constant (this was changed) good for <= 1e6. No significant change to overall time.

#include <ctime>
#include <iostream>
typedef unsigned long long uint64_t;
#define MOD (1000000007ull)
static clock_t dwTimeStart;     // clock values
static clock_t dwTimeStop;

static inline uint64_t PowMod(uint64_t a, uint64_t b){
register uint64_t x=1,y=a;
while(b){
if(b&1){
x*=y; if(x>=MOD)x%=MOD;
}
y*=y; if(y>=MOD)y%=MOD;
b>>=1;
}
return x;
}

static inline uint64_t InverseMod(uint64_t n){
return PowMod(n,MOD-2);
}

static inline uint64_t ProdMod(uint64_t minx, uint64_t maxx){
uint64_t fact = minx;
uint64_t prod;
uint64_t i = (minx+1);
while(i<=maxx){
prod = i++;
while(prod <= 9223372036854ull && i <= maxx)
prod*=i++;
prod%=MOD;
fact=(fact*prod)%MOD;
}
return fact;
}

static inline uint64_t FactMod(uint64_t n){
register uint64_t fact=1;
register uint64_t prod;
register uint64_t i;
i = 1;
while(i <= n){
prod = i++;
while(prod <= 18446744073708ull && i <= n)
prod*=i++;
prod%=MOD;
fact=(fact*prod)%MOD;
}
return fact;
}

static inline uint64_t CheMin(uint64_t x, uint64_t y){
return (ProdMod(x+1,x+y) * InverseMod(FactMod(y)))%MOD;
}

int main()
{
uint64_t i, j;
dwTimeStart = clock();
for(j = 999800; j <= 1000000; j++)
i = CheMin(j, j);
dwTimeStop = clock();
std::cout << "Number of ticks    " << (dwTimeStop-dwTimeStart) << std::endl;
std::cout << "Answer             " << i << std::endl;
return 0;
}
• Comments are not for extended discussion; this conversation has been moved to chat. – Jamal Nov 29 '14 at 22:49

Just some further optimization attempts: The important thing is to mimimize the number of modular reductions. rcgldr's approach tries to do it dynamically, which doesn't seem to be worth it. With numbers somewhere below 2e6, you can multiply exactly three of them; no need to test anything.

Their product must be reduced, then it gets multiplied to the result, which in turn has to be reduced in order to be ready for the next multiplication. This leaves us with 2 reductions every 3 numbers.

I believe, we can do better, namely just 1 reduction for 2 numbers (warning, Java code ahead!):

private long prodMod(int i, int j) {
long result = 1;
long p = mulMod(i, j);
while (i<j-2) {
assert p == mulMod(i, j);
result = mulMod(result, p);
p += j-- - i++ - 1;

if (p>MOD) p -= MOD;
assert p < MOD;
}

assert i <= j;
while (i<=j) result = mulMod(result, i++);
return result;
}

Here, mulMod does exactly what its name says. In every step, result *= i++ * j-- gets computed using modular arithmetic. Once in a blue moon, p must be reduced.

This looks fine, but is not fast, probably due to pipelining issues. The value of result gets change in every iteration, the division takes ages, and time gets wasted due to data dependencies. With some accumulators for intermediate results it got way faster

private long prodMod(int i, int j) {
long acc0 = 1;
long acc1 = 1;
long acc2 = 1;
long acc3 = 1;
long acc4 = 1;
long acc5 = 1;
long p = mulMod(i, j);
while (i<j-20) { // 20 is surely big enough for this unrolling
assert p == mulMod(i, j);
acc0 = mulMod(acc0, p);
p += j-- - i++ - 1;

assert p == mulMod(i, j);
acc1 = mulMod(acc1, p);
p += j-- - i++ - 1;

assert p == mulMod(i, j);
acc2 = mulMod(acc2, p);
p += j-- - i++ - 1;

assert p == mulMod(i, j);
acc3 = mulMod(acc3, p);
p += j-- - i++ - 1;

assert p == mulMod(i, j);
acc4 = mulMod(acc4, p);
p += j-- - i++ - 1;

assert p == mulMod(i, j);
acc5 = mulMod(acc5, p);
p += j-- - i++ - 1;

if (p>MOD) p -= MOD;
assert p < MOD;
}

assert i <= j;
long result = 1;
result = mulMod(result, acc0);
result = mulMod(result, acc1);
result = mulMod(result, acc2);
result = mulMod(result, acc3);
result = mulMod(result, acc4);
result = mulMod(result, acc5);
while (i<=j) result = mulMod(result, i++);
return result;
}

On my Core i5-2400 it takes 4.5 s for 2000 iterations (Java needs some warmup, measuring just 200 iterations would make no sense).

Caching

For caching every 1024 steps, add

private static final int CACHE_SCALE = 10;
private final long[] cache = new long[(2_000_000 >> CACHE_SCALE) + 1];

rename the old method to uncachedProdMod, and add

private long prodMod(int min, int max) {
final int hi = (max >> CACHE_SCALE) << CACHE_SCALE;
long result = uncachedProdMod(hi, max);
int lo = (min >> CACHE_SCALE) << CACHE_SCALE;
if (lo < min) lo += 1 << CACHE_SCALE;
result = mulMod(result, uncachedProdMod(min, lo-1));
result = mulMod(result, cachedProdMod(lo, hi-1));
return result;
}

private long cachedProdMod(int min, int max) {
long result = 1;
while (min < max) {
final int next = min + (1 << CACHE_SCALE);
final int index = min >> CACHE_SCALE;
long prod = cache[index];
if (prod == 0) {
prod = uncachedProdMod(min, next-1);
cache[index] = prod;
}
result = mulMod(result, prod);
min = next;
}
return result;
}

At index i, the cache stores the product of all numbers between i * 2**CACHE_SCALE and (i+1) * 2**CACHE_SCALE - 1 (and 0 denotes an uninitialized slot).

"Memoization" would be probably a better term for this.

• Java=bytecode=1.5GB per process but of good use for arbitrary precision. Python is faster, C++ further more. why using it for modular arithmetics??? – bilbo Nov 30 '14 at 4:13
• @bilbo What 1.5GB??? In general, Python is way slower, but Java's BigInteger is not its strongest part. I wrote it in Java as it's my strongly preferred language. I'd bet a C programmer can read it easily. I guess, my algorithm is faster, given that it makes more iterations per second on a weaker CPU. – maaartinus Nov 30 '14 at 7:11
• On some profilers 1.5GB is the RAM consumption one can see at startup of any Java software. At my stupefaction your code recoded in C++ is faster than all the previous I've tested for the problem at this time: 1.88s congratulations for your optimization. – bilbo Nov 30 '14 at 13:27
• @bilbo I guess that's the virtual memory initially allocated, you can tune it down with something like -Xmx500M. The default is based on total memory available and Java is memory hungry, but I've never cared. I'm glad to hear that the optimization is good, but is it good enough? – maaartinus Dec 1 '14 at 18:09
• unfortunately, it's not enough (time required is 0.9s) – bilbo Dec 1 '14 at 22:00

First you can get your math to work a little so you have less work to do.

If we have $x<y$, and we know that:

$(x+y)! = 1•2•3•...•(y-1)(y)(y+1)...(y+x-1)(y+x) = y!(y+1)...(y+x)$

Then you have:

$\frac{y!(y+1)...(y+x)}{x!y!}=\frac{(y+1)...(y+x)}{x!}$

So you have a lot less work yo do, given that $x<y$.

• If you read my code attentively your idea is what I've implemented in my code. I think the improvement to find is in math modular arithmetics in the computation of (n!)%mod and of 𝙿𝚛𝚘𝚍(y+1,y+x)%mod – bilbo Nov 27 '14 at 17:04