【POJ2118】Firepersons-递推+矩阵乘法

测试地址：Firepersons
题目大意：求解常系数线性递推关系：$a_n=\sum_{i=1}^ka_{n-i}b_i$的第$i$项（$a_i$），其中$a_0,...,a_{k-1}$和$b_1,...,b_k$已给出。
做法：本题需要用到递推+矩阵乘法。
这题$i\le 10^9,k\le 100$，很容易想到用矩阵乘法做到$O(k^3\log i)$的复杂度，那么这题就做完了。
……但是有没有更好的做法呢？
令转移矩阵为$M$，我们最后要求的是$M^{i-k+1}X$（其中$X$为列向量，从下往上排列着$a_0,...,a_{k-1}$）。根据一些奇妙的线性代数结论，由特征多项式以及一些奇妙定理可以得到：
对于$i\ge 0$，有$M^{i+k}=\sum_{j=1}^kb_jM^{i+k-j}$。
有了这个东西，有什么用呢？我们可以得到一个重要的结论：所有$M_i$都可以表示成$M^0,...,M^{k-1}$的线性组合。这个可以用数学归纳法简单证明。
于是我们在原先矩阵快速幂的基础上，将两个$M$的幂乘起来的过程原先是需要$O(k^3)$的，但现在我们只需要做一个多项式乘法即可，时间复杂度为$O(k^2)$。但这样得出来的是$M^0,...,M^{2k-2}$的线性组合，那么只要用上面那个结论把高于$k-1$次的项往下规约即可，也是$O(k^2)$的，因此整个算法的时间复杂度为$O(k^2\log i)$，可以做到$k=2000$左右的水平。
据说用FFT优化多项式乘法等可以优化到$O(k\log k\log i)$……有待学习。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=10000;
int k;
ll a[210],b[110],n;
ll tmp[210],s[110],ss[110];

void mult(ll *A,ll *B,ll *C)
{
    memset(tmp,0,sizeof(tmp));
    for(int i=0;i<k;i++)
        for(int j=0;j<k;j++)
            tmp[i+j]=(tmp[i+j]+A[i]*B[j])%mod;
    for(int i=2*k-2;i>=k;i--)
    {
        for(int j=1;j<=k;j++)
            tmp[i-j]=(tmp[i-j]+tmp[i]*b[j])%mod;
    }
    for(int i=0;i<k;i++)
        C[i]=tmp[i];
}

void solve(ll n)
{
    int nowx=1;
    memset(s,0,sizeof(s));
    s[0]=1;
    memset(ss,0,sizeof(ss));
    if (k>1) ss[1]=1;
    else ss[0]=b[1];

    n=n-k+1;
    while(n)
    {
        if (n&1) mult(s,ss,s);
        n>>=1;nowx<<=1; 
        if (nowx<k)
        {
            for(int i=0;i<k;i++)
                ss[i]=(i==nowx);
        }
        else mult(ss,ss,ss);
    }

    for(int i=k;i<=2*k-2;i++)
    {
        a[i]=0;
        for(int j=1;j<=k;j++)
            a[i]=(a[i]+a[i-j]*b[j])%mod;
    }
    ll ans=0;
    for(int i=0;i<k;i++)
        ans=(ans+a[i+k-1]*s[i])%mod;
    printf("%lld\n",ans);
}

int main()
{
    while(scanf("%d",&k)&&k)
    {
        for(int i=0;i<k;i++)
            scanf("%lld",&a[i]);
        for(int i=1;i<=k;i++)
            scanf("%lld",&b[i]);
        scanf("%lld",&n);

        if (n<k) printf("%lld\n",a[n]);
        else solve(n);
    }

    return 0;
}