BZOJ 3782: 上学路线

3782: 上学路线

Time Limit: 10 Sec  Memory Limit: 128 MB
Submit: 211  Solved: 82
[Submit][Status][Discuss]

Description

小C所在的城市的道路构成了一个方形网格，它的西南角为(0,0)，东北角为(N,M)。小C家住在西南角，学校在东北角。现在有T个路口进行施工，小C不能通过这些路口。小C喜欢走最短的路径到达目的地，因此他每天上学时都只会向东或北行走；而小C又喜欢走不同的路径，因此他问你按照他走最短路径的规则，他可以选择的不同的上学路线有多少条。由于答案可能很大，所以小C只需要让你求出路径数mod P的值。

Input

第一行，四个整数N、M、T、P。

接下来的T行，每行两个整数，表示施工的路口的坐标。

Output

一行，一个整数，路径数mod P的值。

 

 

Sample Input

3 4 3 1019663265
3 0
1 1
2 2

Sample Output

8

HINT

1<=N,M<=10^10

0<=T<=200

p=1000003或p=1019663265

Source

 

[Submit][Status][Discuss]

Lucas定理 + 中国剩余定理CRT

#include <cstdio>
#include <algorithm>

#define int long long

const int mxt = 205;

int n, m, t, p, f[mxt];

struct point
{
    int x, y;
}pt[mxt];

inline bool cmp(const point &a, const point &b)
{
    if (a.x != b.x)
        return a.x < b.x;
    else
        return a.y < b.y;
}

namespace case1    // p = 1000003
{
    const int siz = 1000005;
    
    int fac[siz], inv[siz];
    
    int lucas(int a, int b)
    {
        if (a < b)
            return 0;
        
        if (a < p && b < p)
            return fac[a] * inv[b] % p * inv[a - b] % p;
        
        return lucas(a / p, b / p) * lucas(a % p, b % p) % p;
    }
    
    inline void prework(void)
    {
        fac[0] = fac[1] = 1;
        inv[0] = inv[1] = 1;
        
        for (int i = 2; i <= p; ++i)
        {
            fac[i] = fac[i - 1] * i % p;
            inv[i] = inv[p % i] * (p - p/i) % p;
        }
        
        for (int i = 2; i <= p; ++i)
            inv[i] = inv[i] * inv[i - 1] % p;
    }
    
    inline int calc(int a, int b)
    {
        return lucas(a + b, a);
    }
    
    inline void main(void)
    {
        prework();
        
        int ans = calc(n, m);
        
        for (int i = 1; i <= t; ++i)
        {
            f[i] = calc(pt[i].x, pt[i].y);
            
            for (int j = 1; j < i; ++j)
                if (pt[j].x <= pt[i].x && pt[j].y <= pt[i].y)
                    f[i] = (f[i] - f[j] * calc(pt[i].x - pt[j].x, pt[i].y - pt[j].y) % p + p) % p;
        }
        
        for (int i = 1; i <= t; ++i)
            if (pt[i].x <= n && pt[i].y <= m)
                ans = (ans - f[i] * calc(n - pt[i].x, m - pt[i].y) % p + p) % p;
        
        printf("%lld
", ans);
    }
}

namespace case2    // p = 3 * 5 * 6793 * 10007
{
    const int siz = 10010;
    
    const int mod[4] = 
    {
        3, 5, 6793, 10007
    };
    
    int fac[4][siz], inv[4][siz];
    
    inline void prework(void)
    {
        for (int i = 0; i < 4; ++i)
        {
            fac[i][0] = fac[i][1] = 1;
            inv[i][0] = inv[i][1] = 1;
            
            for (int j = 2; j <= mod[i]; ++j)
            {
                fac[i][j] = fac[i][j - 1] * j % mod[i];
                inv[i][j] = inv[i][mod[i] % j] * (mod[i] - mod[i]/j) % mod[i];
            }
            
            for (int j = 2; j <= mod[i]; ++j)
                inv[i][j] = inv[i][j - 1] * inv[i][j] % mod[i];
        }
    }
    
    int lucas(int a, int b, int k)
    {
        if (a < b)
            return 0;
            
        if (a < mod[k] && b < mod[k])
            return fac[k][a] * inv[k][b] % mod[k] * inv[k][a - b] % mod[k];
        
        return lucas(a / mod[k], b / mod[k], k) * lucas(a % mod[k], b % mod[k], k) % mod[k];
    }
    
    inline int calc(int a, int b)
    {
        return (
            lucas(a + b, a, 0) * 339887755 % p
        +    lucas(a + b, a, 1) * 407865306 % p
        +    lucas(a + b, a, 2) * 673070820 % p
        +    lucas(a + b, a, 3) * 618502650 % p
        ) % p;
    }
    
    inline void main(void)
    {
        prework();
        
        int ans = calc(n, m);
        
        for (int i = 1; i <= t; ++i)
        {
            f[i] = calc(pt[i].x, pt[i].y);
            
            for (int j = 1; j < i; ++j)
                if (pt[j].x <= pt[i].x && pt[j].y <= pt[i].y)
                    f[i] = (f[i] - f[j] * calc(pt[i].x - pt[j].x, pt[i].y - pt[j].y) % p + p) % p;
        }
        
        for (int i = 1; i <= t; ++i)
            if (pt[i].x <= n && pt[i].y <= m)
                ans = (ans - f[i] * calc(n - pt[i].x, m - pt[i].y) % p + p) % p;
        
        printf("%lld
", ans);
    }
}

signed main(void)
{
    scanf("%lld%lld%lld%lld", &n, &m, &t, &p);
    
    for (int i = 1; i <= t; ++i)
        scanf("%lld%lld", &pt[i].x, &pt[i].y);
    
    std::sort(pt + 1, pt + t + 1, cmp);
    
    if (p == 1000003)
        case1::main();
    else
        case2::main();
}

@Author: YouSiki