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=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=MOD)ans0%=MOD;
      ans1*=((i=MOD)ans1%=MOD;
    }
    return ans0*ans1%MOD;
  }else{
    for(i=1; i=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=MOD)res%=MOD;
   }
   res=InverseMod(res);   
    if(!(n&1))
      res= -res +MOD;
  }
  return res%MOD;
}

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 
#include 
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;
}

StackExchange Code Review Q#70929, answer score: 2