C++のコードですがこのようにすることもできます。

    //aとbの最大公約数
    int gcd(int a, int b) {
    	switch (b) {
    	case 0: return a;
    	default: return gcd(b, a % b);
    	}
    }
    long long int comb(int n, int r) {
    	//n!
    	std::vector<int> numerator;
    	for (int i = n; i > 0; --i)
    		numerator.push_back(i);
    
    	//numeratorをr!で割る
    	for (auto i = 1; i <= r; ++i) {
    		int d = i;
    		for (int j = 0; j < numerator.size(); ++j) {
    			int g = gcd(d, numerator[j]);
    			d /= g;
    			numerator[j] /= g;
    		}
    	}
    	//numeratorを(n-r)!で割る
    	for (int i = 1; i <= n - r; ++i) {
    		int d = i;
    		for (int j = 0; j < numerator.size(); ++j) {
    			int g = gcd(d, numerator[j]);
    			d /= g;
    			numerator[j] /= g;
    		}
    	}
    
    	//numeratorを掛け合わせる。
    	long long int result = 1;
    	for (int i = 0; i < numerator.size(); ++i) {
    		result = (result * numerator[i]) % MOD;
    	}
    	return result;
    }

基本的な考えは、分子の階乗部分と分母の階乗部分を掛け合わせる前に約分していくというものです。


----------
今回の問題に関しては明らかに不要な効率化ですが、2と5の素因数を含まない数についてオイラーの公式等を使って計算して、後で2と5を掛けるという方法もあります。

    long long int power(long long int base, long long int exp) {
    	switch (exp) {
    	case 0: return 1;
    	case 1: return base;
    	default: return power(base * base % MOD, exp / 2) * power(base, exp % 2) % MOD;
    	}
    }
    //n!をfactorで割ることができる回数
    int count_factor(int n, int factor) {
    	switch (n) {
    	case 0: return 0;
    	default: return n / factor + count_factor(n / factor, factor);
    	}
    }
    int remove2_5(int n) {
    	if (n % 2 == 0) return remove2_5(n / 2);
    	if (n % 5 == 0) return remove2_5(n / 5);
    	return n;
    }
    long long int comb(int n, int r) {
    	long long int numerator = 1, denominator = 1;
    	for (auto i = 1; i <= n; ++i) {
    		numerator = numerator * remove2_5(i) % MOD;
    	}
    	for (auto i = 1; i <= r; ++i) {
    		denominator = denominator * remove2_5(i) % MOD;
    	}
    	for (auto i = 1; i <= n - r; ++i) {
    		denominator = denominator * remove2_5(i) % MOD;
    	}
    	return numerator * power(denominator, 400000000 - 1) % MOD *
    		power(2, count_factor(n, 2) - count_factor(r, 2) - count_factor(n - r, 2)) % MOD *
    		power(5, count_factor(n, 5) - count_factor(r, 5) - count_factor(n - r, 5)) % MOD;
    }

これだとn = 10^6ぐらいでも高速に求められます。

    comb(1234567, 5000) #=> 772988800