HDU 2254 奥运

/*
*  [题意]
*   给出n条道路，k个询问，每个询问包括起点v1、终点v2、t1天、t2天
*   问从v1到v2走了i天一共有多少走法(mod 2008)？(t1<=i<=t2)
*  [解题方法]
*   设B = A^i;
*   则A[u][v] 表示 从u到v走了i天(等价于走了i条边)的走法有多少
*   那么题目就转化为求：C = (A^t1+A^(t1+1)+...+A^t2) % 2008;
*   C[v1][v2]即为所求
*/
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <map>
using namespace std;
#define LL long long
#define FF(i, n) for(int i = 0; i < n; i++)
#define M 31

int mod = 2008, n;

struct mat {
    int x[M][M];
};

mat matadd(const mat &a, const mat &b)
{
    mat res;
    FF(i, n) FF(j, n)
        res.x[i][j] = (a.x[i][j]+b.x[i][j]) % mod;
    return res;
}

mat matmul(const mat &a, const mat &b)
{
    mat res;
    FF(i, n) FF(j, n) res.x[i][j] = 0;
    FF(i, n) FF(k, n) if(a.x[i][k]) FF(j, n) if(b.x[k][j])
        res.x[i][j] = (res.x[i][j]+(LL)a.x[i][k]*b.x[k][j]%mod) % mod;
    return res;
}

mat qmod(mat a, int b)
{
    mat res;
    FF(i, n) FF(j, n) res.x[i][j] = (i==j);
    for ( ; b; b >>= 1)
    {
        if (b & 1) res = matmul(res, a);
        a = matmul(a, a);
    }
    return res;
}

mat cal(const mat &a, int k)    //分治求(a^0+a^1+a^2+...+.a^k)%mod
{
    if (k <= 0) return qmod(a, 0);
    int b = (k+1)/2;
    mat o = cal(a, b-1);
    mat res = matadd(o, matmul(qmod(a, b), o));
    if (k % 2 == 0) res = matadd(res, qmod(a, k));
    return res;
}

map<int, int> Map;

int main()
{
    mat A;
    int t1, t2, a, b, t;
    while (~scanf("%d", &t))
    {
        FF(i, M) FF(j, M) A.x[i][j] = 0;
        Map.clear();
        n = 0;
        while (t--)
        {
            scanf("%d%d", &a, &b);
            if (!Map[a]) a = Map[a] = ++n;
            else a = Map[a];
            if (!Map[b]) b = Map[b] = ++n;
            else b = Map[b];
            ++A.x[a-1][b-1];    //注意重边
        }
        scanf("%d", &t);
        while (t--)
        {
            scanf("%d%d%d%d", &a, &b, &t1, &t2);
            if (!Map[a] || !Map[b]) {
                puts("0");
                continue;
            }
            mat m1 = cal(A, t1-1);
            mat m2 = cal(A, t2);
            a = Map[a] - 1;
            b = Map[b] - 1;
            printf("%d\n", (m2.x[a][b]-m1.x[a][b]+mod)%mod);
        }
    }
    return 0;
}