HUNAN 11562 The Triangle Division of the Convex Polygon(大卡特兰数)

  1 #include <cstdio>
  2 #include <cmath>
  3 #include <stdlib.h>
  4 #include <memory.h>
  5 typedef int typec;
  6 typec GCD(typec a, typec b)
  7 {
  8     return b? GCD(b, a % b) : a;
  9 }
 10 typec extendGCD(typec a, typec b, typec& x, typec& y)
 11 {
 12     if(!b) return x = 1, y = 0, a;
 13     typec res = extendGCD(b, a % b, x, y), tmp = x;
 14     x = y, y = tmp -(a/b)*y;
 15     return res;
 16 }
 17 typec power(typec x, typec k)
 18 {
 19     typec res = 1;
 20     while(k)
 21     {
 22         if(k&1) res *= x;
 23         x *= x, k >>= 1;
 24     }
 25     return res;
 26 }
 27 typec powerMod(typec x, typec k, typec m)
 28 {
 29     typec res = 1;
 30     while(x %= m, k)
 31     {
 32         if(k&1) res *= x, res %= m;
 33         x *= x, k >>= 1;
 34     }
 35     return res;
 36 }
 37 typec inverse(typec a, typec p, typec t = 1)
 38 {
 39     typec pt = power(p, t);
 40     typec x, y;
 41     y = extendGCD(a, pt, x, y);
 42     return x < 0? x += pt : x;
 43 }
 44 typec linearCongruence(typec a, typec b, typec p, typec q)
 45 {
 46     typec x, y;
 47     y = extendGCD(p, q, x, y);
 48     x *= b - a, x = p * x + a, x %= p * q;
 49     if(x < 0) x += p * q;
 50     return x;
 51 }
 52 const int PRIMEMAX = 1000;
 53 int prime[PRIMEMAX + 1];
 54 int getPrime()
 55 {
 56     memset(prime, 0, sizeof(int) * (PRIMEMAX + 1));
 57     for(int i = 2; i <= PRIMEMAX; i++)
 58     {
 59         if(!prime[i]) prime[++prime[0]] = i;
 60         for(int j = 1; j <= prime[0] && prime[j] <= PRIMEMAX/i; j++)
 61         {
 62             prime[prime[j]*i] = 1;
 63             if(i % prime[j] == 0) break;
 64         }
 65     }
 66     return prime[0];
 67 }
 68 int factor[100][3], facCnt;
 69 int getFactors(int x)
 70 {
 71     facCnt = 0;
 72     int tmp = x;
 73     for(int i = 1; prime[i] <= tmp / prime[i]; i++)
 74     {
 75         factor[facCnt][1] = 1, factor[facCnt][2] = 0;
 76         if(tmp % prime[i] == 0)
 77             factor[facCnt][0] = prime[i];
 78         while(tmp % prime[i] == 0)
 79             factor[facCnt][2]++, factor[facCnt][1] *= prime[i], tmp /= prime[i];
 80         if(factor[facCnt][2]) facCnt++;
 81     }
 82     if(tmp != 1) factor[facCnt][0] = tmp, factor[facCnt][1] = tmp, factor[facCnt++][2] = 1;
 83     return facCnt;
 84 }
 85 typec combinationModPt(typec n, typec k, typec p, typec t = 1)
 86 {
 87     if(k > n) return 0;
 88     if(n - k < k) k = n - k;
 89     typec pt = power(p, t);
 90     typec a = 1, b = k + 1, x, y;
 91     int pcnt = 0;
 92     while(b % p == 0) pcnt--, b /= p;
 93     b %= pt;
 94     for(int i = 1; i <= k; i++)
 95     {
 96         x = n - i + 1, y = i;
 97         while(x % p == 0) pcnt++, x /= p;
 98         while(y % p == 0) pcnt--, y /= p;
 99         x %= pt, y %= pt, a *= x, b *= y;
100         a %= pt, b %= pt;
101     }
102     if(pcnt >= t) return 0;
103     extendGCD(b, pt, x, y);
104     if(x < 0) x += pt;
105     a *= x, a %= pt;
106     return a * power(p, pcnt) % pt;
107 }
108 const typec PTMAX = 45000;
109 typec facmod[PTMAX];
110 void initFacMod(typec p, typec t = 1)
111 {
112     typec pt = power(p, t);
113     facmod[0] = 1 % pt;
114     for(int i = 1; i < pt; i++)
115     {
116         if(i % p) facmod[i] = facmod[i - 1] * i % pt;
117         else facmod[i] = facmod[i - 1];
118     }
119 }
120 typec factorialMod(typec n, typec &pcnt, typec p, typec t = 1)
121 {
122     typec pt = power(p, t), res = 1;
123     typec stepCnt = 0;
124     while(n)
125     {
126         res *= facmod[n % pt], res %= pt;
127         stepCnt += n /pt, n /= p, pcnt += n;
128     }
129     res *= powerMod(facmod[pt - 1], stepCnt, pt);
130     return res %= pt;
131 }
132 typec combinationModPtFac(typec n, typec k, typec p, typec t = 1)
133 {
134     if(k > n || p == 1) return 0;
135     if(n - k < k) k = n - k;
136     typec pt = power(p, t), pcnt = 0, pmcnt = 0;
137     if(k < pt) return combinationModPt(n, k, p, t);
138     initFacMod(p, t);
139     typec a = factorialMod(n, pcnt, p, t);
140     typec b = factorialMod(k, pmcnt, p, t);
141     b *= b, pmcnt <<= 1, b %= pt;
142     typec tmp = k + 1;
143     while(tmp % p == 0) tmp /= p, pmcnt++;
144     b *= tmp % pt, b %= pt;
145     pcnt -= pmcnt;
146     if(pcnt >= t) return 0;
147     a *= inverse(b, p, t), a %= pt;
148     return a * power(p, pcnt) % pt;
149 }
150 typec combinationModFac(typec n, typec k, typec m)
151 {
152     getFactors(m);
153     typec a, b, p, q;
154     for(int i = 0; i < facCnt; i++)
155     {
156         if(!i) a = combinationModPtFac(n, k, factor[i][0], factor[i][2]), p = factor[i][1];
157         else b = combinationModPtFac(n, k, factor[i][0], factor[i][2]), q = factor[i][1];
158         if(!i) continue;
159         a = linearCongruence(a, b, p, q), p *= q;
160     }
161     return a;
162 }
163 int main()
164 {
165     getPrime();
166     typec n, k;
167     while(scanf("%d %d", &n, &k) != EOF)
168         printf("%d\n", combinationModFac(2 * (n-2), n-2, k));
169     return 0;
170 }