HDU 5226 Tom and matrix（组合数学+Lucas定理）

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=5226

题意：给一个矩阵a，a[i][j] = C（i，j）(i>=j) or 0(i < j)，求(x1，y1)，(x2，y2)这个子矩阵里面的所有数的和。

思路：首先可以推导出一个公式C（n，i）+C（n + 1，i）+...+C（m，i） = C（m + 1，i + 1）

知道了这个公式，就可以将子矩阵里每行（或每列）的和值表示成组合数的差值，现在的关键是求出C（n，m）（mod p）.

由于n和m可能很大，p很小，不能直接求，要借助Lucas定理。关于Lucas定理，可参考：http://www.cnblogs.com/zxndgv/archive/2011/09/17/2179591.html。

code：

 1 #include <cstdio>
 2 typedef __int64 LL;
 3 const int MAXN = 100005;
 4 int p;
 5 LL fac[MAXN];
 6
 7 // 得到阶乘  fac[i] = i! % p
 8 void GetFact()
 9 {
10     fac[0] = 1;
11     for (LL i = 1; i < MAXN; ++i)
12         fac[i] = fac[i - 1] * i % p;
13 }
14
15 // 快速模幂 a^b % p
16 LL Pow(LL a, LL b)
17 {
18     LL temp = a % p;
19     LL ret = 1;
20     while (b)
21     {
22         if (b & 1) ret = ret * temp % p;
23         temp = temp * temp % p;
24         b >>= 1;
25     }
26     return ret;
27 }
28
29 /*
30 欧拉定理求逆元
31 (a / b) (mod p) = (a * x) (mod p) x表示b的逆元 并且 b * x = 1 (mod p) 只有b和p互质才存在逆元
32
33 b * x = 1 (mod p) x是b关于p的逆元
34
35 b^phi(p) = 1 (mod p)
36
37 b * b^(phi(p) - 1) (mod p) = b * x (mod p)
38
39 x = b^(phi(p) - 1) = b^(p - 2)
40
41 (a / b) (mod p) = (a * x) (mod p) = (a * b^(p - 2)) (mod p)
42
43 经过上面的推导，得出：
44
45 (a / b) (mod p) = (a * b^(p - 2)) (mod p) (b 和 p互质)
46
47 */
48 LL Cal(LL n, LL m)
49 {
50     if (m > n) return 0;
51     return fac[n] * Pow(fac[m] * fac[n - m], p - 2) % p;
52 }
53
54 LL Lucas(LL n, LL m)
55 {
56     if (m == 0) return 1;
57     return Cal(n % p, m % p) * Lucas(n / p, m / p) % p;
58 }
59
60 int main()
61 {
62     int x1, y1, x2, y2;
63     while (scanf("%d %d %d %d %d", &x1, &y1, &x2, &y2, &p) == 5)
64     {
65         if (x2 < y1)    // 预判 子矩阵全部0值区域
66         {
67             printf("0\n");
68             continue;
69         }
70         if (x2 == y1)   // 预判 子矩阵只有右上角值为1，其余为0
71         {
72             printf("1\n");
73             continue;
74         }
75         GetFact();
76         if (x1 < y1) x1 = y1;
77         if (y2 > x2) y2 = x2;
78         LL ans = 0;
79         for (int i = y1; i <= y2; ++i)
80         {
81             if (i > x1)
82                 ans = (ans + Lucas(x2 + 1, i + 1)) % p;
83             else
84                 ans = (ans + Lucas(x2 + 1, i + 1) - Lucas(x1 + 1, i + 1) + Lucas(x1, i)) % p;
85         }
86         printf("%I64d\n", ans);
87     }
88     return 0;
89 }