【BZOJ 3907】网格（Catalan数）

题目链接
这个题推导公式跟\(Catalan\)数是一样的，可得解为\(C_{n+m}^n-C_{n+m}^{n+1}\)
然后套组合数公式\(C_n^m=\frac{n!}{m!(n-m)!}\)
用阶乘分解的方法对分子和分母分解质因数然后指数相减，最后把剩下的高精度乘起来就行了，这样就避免了高精除法。可以用快速幂，但我太lan了，就直接暴力乘起来了。
说一下怎么阶乘分解，直接对每个数分解质因数的时间复杂度是\(O(n\sqrt{n})\)，这显然是不可忍受的。
于是，考虑先用线筛求出\(1~n\)之间所有质数，然后枚举所有质数\(p\)，\(1~n\)中所有\(p\)的倍数包含一个质因子\(p\)，所有\(p^2\)的倍数包含一个质因子\(p\)，...
所以，\(n!\)中\(p\)的指数为\(\sum_{i=1}^{\left\lfloor \log_p n\right\rfloor}\left\lfloor n/p^i \right\rfloor\)
枚举\(p\)求就行了，时间复杂度\(O(n\ log\ n)\)。

Code:

#include <cstdio>
#include <cstring>
const int MAXN = 10010;
int prime[MAXN], v[MAXN], c[MAXN], d[MAXN];
int n, cnt, m;
const int MOD = 10000;
struct Int{
    int s[MAXN];
    Int(int x){ memset(s, 0, sizeof s); s[++s[0]] = x % 10, x /= 10; } 
    void mul(int x){
        s[1] *= x;
        for(int i = 2; i <= s[0]; ++i){
           s[i] = s[i] * x + s[i - 1] / MOD;
           s[i - 1] %= MOD;
        }
        while(s[s[0]] >= MOD){
          s[++s[0]] = s[s[0] - 1] / MOD;
          s[s[0] - 1] %= MOD;
        }
    }
    void cut(Int x){
        for(int i = 1; i <= x.s[0]; ++i)
           s[i] -= x.s[i];
        for(int i = 1; i <= x.s[0]; ++i){
           if(s[i] < 0){
             s[i] += MOD;
             --s[i + 1];
           }
        }
        if(!s[s[0]]) --s[0];
    }
    void print(){
        printf("%d", s[s[0]]);
        for(int i = s[0] - 1; i; --i)
           printf("%04d", s[i]);
    }
}a(1), b(1);
int main(){
    scanf("%d%d", &n, &m);
    for(int i = 2; i <= n + m; ++i){
       if(!v[i]){
         v[i] = i;
         prime[++cnt] = i;
       }  
       for(int j = 1; j <= cnt; ++j){
          if(prime[j] > v[i] || prime[j] * i > n + m) break;
          v[prime[j] * i] = prime[j];
       }
    }
    for(int i = 1; i <= cnt; ++i)
       for(int j = prime[i]; j <= n + m; j *= prime[i])
          c[i] += (n + m) / j;
    for(int i = 1; i <= cnt && prime[i] <= n; ++i)
       for(int j = prime[i]; j <= n; j *= prime[i])
          c[i] -= n / j;
    for(int i = 1; i <= cnt && prime[i] <= m; ++i)
       for(int j = prime[i]; j <= m; j *= prime[i])
          c[i] -= m / j;
    for(int i = 1; i <= cnt; ++i)
       for(int j = 1; j <= c[i]; ++j)
          a.mul(prime[i]);
          
    for(int i = 1; i <= cnt; ++i)
       for(int j = prime[i]; j <= n + m; j *= prime[i])
          d[i] += (n + m) / j;
    for(int i = 1; i <= cnt && prime[i] <= n + 1; ++i)
       for(int j = prime[i]; j <= n + 1; j *= prime[i])
          d[i] -= (n + 1) / j;
    for(int i = 1; i <= cnt && prime[i] <= m - 1; ++i)
       for(int j = prime[i]; j <= m - 1; j *= prime[i])
          d[i] -= (m - 1) / j;
    for(int i = 1; i <= cnt; ++i)
       for(int j = 1; j <= d[i]; ++j)
          b.mul(prime[i]);
    a.cut(b);
    a.print();
    return 0;
}