每日算法--矩阵乘法优化递推

　　今日开坑每日算法，慢慢填吧，争取省选前每日一更。

　　---把学过的讲一遍，就相当于巩固一遍（伟大的SymenYang。。。。）

　　矩阵大家应该比较熟悉，至少熟悉其基本的运算规律。比如乘法，加法等（不懂的递归调用矩阵），对于矩阵乘法，Cij=∑（aik*bkj） k∈[1,n];因为矩阵乘法满足结合律，所以当求矩阵的幂时，可以使用快速幂加快速度，从而达到logN。。

　　矩阵可以用来优化递推，先来看一个简单的例子：斐波拉契数列!

　　递推式：f[i]为数列第i项 f[i] = f[i-1] + f[i-2]; （定义域。。。）

　　观察递推式，f[i] = 1*某个数 + 1*某个数，是否有点像矩阵乘法的定义，仔细推理一下，如何构造这个矩阵呢（略过，直接给出答案）

　　

|f[i],f[i-1]|             |1     ,    1|                |f[i]+f[i-1]  ,f[i]|
|0   ,0     |      *      |1     ,    0|        ==      |0            ,   0|

　　也许会问，为什么下面两个数都为0还要开一个2*2的矩阵呢，主要是为了快速幂的方便，一个可以和自己乘上许多次（>=2）的矩阵只有可能是正方形的，所以要开这样一个矩阵。

　　接下来是一道进阶题目——BZOJ1297 SCOI2009 迷路

　　首先设出状态：f[i][j]表示在i时到达j的方案数，很容易得知 f[i][j]=∑f[i-k][m] dis[m][j]==k;

　　但是这个状态似乎是有点麻烦，因为f[i]可能由许多层以前推出，不能实现矩阵乘法。

　　考虑到边权<=10，我们可以把一个点拆成多个（最多10个）,然后使相连点之间的距离都变为了1，如此f[i][j]=f[i-1][m] （存在m->j的路径），将每一刻的方案数作为矩阵即可，代码：

/**************************************************************
    Problem: 1297
    User: SymenYang
    Language: C++
    Result: Accepted
    Time:1872 ms
    Memory:5596 kb
    PS：无视矩阵英文的拼写错误吧，没有写juzhen就很不错了。。。
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstdlib>
#include <cstring>
#include <string>
using namespace std;
struct martix
{
    int data[110][110];
    martix()
    {
        memset(data,0,sizeof(data));
    }
};
const int mod=2009;
martix cf(martix a,martix b,int n)
{
    martix ret;
    for (int i=0;i<=n;i++)
    {
        for (int j=0;j<=n;j++)
        {
            for (int k=0;k<=n;k++)
            {
                ret.data[i][j]=(ret.data[i][j]+(a.data[i][k]*b.data[k][j])%mod)%mod;
            }
        }
    }
    return ret;
}
 
 
void pi(martix a,int n)
{
    for(int i=0;i<=n;i++)
    {
        for (int j=0;j<=n;j++)
        {
            printf("%d ",a.data[i][j]);
        }
        printf("
");
    }
    printf("-------------------------------------------------------------
");
}
 
martix qp(martix a,long long b,int n)
{
    if (b==1) return a;
    martix c=qp(a,b>>1,n);
//  pi(c,n);
    c=cf(c,c,n);
//  pi(c,n);
    if (b%2==1) c=cf(c,a,n);
//  pi(c,n);
    return c;
}
 
 
int map[10][10];
int main()
{
    martix ans;
    int n,t;
    scanf("%d%d",&n,&t);
    char ch[20];
    for (int i=1;i<=n;i++)
    {
        scanf("%s",ch);
        for (int j=0;j<n;j++)
        {
            int a=ch[j]-'0';
            map[i][j+1]=a;
        }
    }
    for (int i=1;i<=n;i++)
    {
        for (int j=1;j<=n;j++)
        {
            if (map[i][j])
            {
                for (int k=(i-1)*9+1;k<=(i-1)*9+map[i][j];k++)
                {
                    if (k==(i-1)*9+map[i][j])
                    {
                        ans.data[k][(j-1)*9+1]=1;
                    }
                    else
                    {
                        ans.data[k][k+1]=1;
                    }
                }
            }
        }
    }
     
    ans.data[0][1]=1;
    ans=qp(ans,t+1,9*n);
    printf("%d",ans.data[0][n*9-9+1]);
}

　　我们发现，当递推关系较大时（10^18） 很容易想到矩阵快速幂，但如果再大呢？比如10^1000000？

　　费尔马小定理+矩阵——BZOJ3240  NOI2013 矩阵游戏

　　首先，可以看出这是一个裸的矩阵题，是不是感觉很水，秒掉？ 再看数据范围，出现了10^1000000.。。。。，似乎写高精度?还是太慢，发现mod1000000007，质数，恩不错。

　　费尔马小定理递归学习

　　a^(p-1) ≡ 1 （mod p） p 是质数

　　那么 a^k ==  a^(k%(p-1)) * (a^(p-1)) ^ (k/(p-1));

　　又a^(p-1) ≡ 1 （mod p） 所以a^k ≡ a^(k%(p-1)) （mod p） p 是质数

　　似乎有点意思，我们每读入一个数字，就使n=（n*10+num）% p 即可。。

　　但是！！！！！！

　　如果乘的时候比例==1呢，费马小不适用了！！

　　判断一下，为1时将p-1改为p即可。

　　代码：

/**************************************************************
    Problem: 3240
    User: Life_is_Beautiful
    Language: C++
    Result: Accepted
    Time:400 ms
    Memory:3228 kb
    PS：看吧，叫juzhen了吧
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstdlib>
#include <cstring>
using namespace std;
struct juzhen
{
    long long data[3][3];
     
    juzhen(long long a,long long b,long long c,long long d)
    {
        data[1][1]=a;
        data[1][2]=b;
        data[2][1]=c;
        data[2][2]=d;
    }
     
    juzhen(long long a[3][3])
    {
        for (int i=1;i<=2;i++)
            for (int j=1;j<=2;j++)
            {
                data[i][j]=a[i][j];
            }
    }
    juzhen(){}
};
const long long mod=1000000007;
juzhen cf(juzhen a ,juzhen b)
{
    long long c[3][3];
    for (int i=1;i<=2;i++)
    {
        for (int j=1;j<=2;j++)
        {
            c[i][j]=0;
            for (int k=1;k<=2;k++)
            {
                c[i][j]=(c[i][j]+(a.data[i][k]*b.data[k][j])%mod)%mod;
            }
        }
    }
    return juzhen(c);
}
 
juzhen qp(juzhen a,long long b)
{
    if (b==0) return juzhen(1,0,0,1);
    juzhen c=qp(a,b>>1);
    c=cf(c,c);
    if ((b%2)==1) c=cf(c,a);
    return c;
}
char ch[1000010];
char ch2[1000010];
int main()
{
    long long n,m,a,b,c,d;
    scanf("%s",ch);
    int l=strlen(ch);
    n=0;
    for (int i=0;i<l;++i)
    {
            n=(n*10+ch[i]-'0')%(mod-1);
    }
    scanf("%s",ch2);
    l=strlen(ch2);
    m=0;
    for (int i=0;i<l;++i)
    {
            m=(m*10+ch2[i]-'0')%(mod-1);
    }
    cin>>a>>b>>c>>d;
    if (a==1&&c==1)
    {
        n=0;
        for (int i=0;i<l;++i)
        {
                n=(n*10+ch[i]-'0')%mod;
        }
        l=strlen(ch2);
        m=0;
        for (int i=0;i<l;++i)
        {
                m=(m*10+ch2[i]-'0')%mod;
        }
    }
    juzhen ans(a,0,b,1);
    ans=qp(ans,m-1);
    juzhen cns=ans;
    ans=cf(ans,juzhen(c,0,d,1));
    ans=qp(ans,n-1);
    ans=cf(ans,cns);
    ans=cf(juzhen(1,1,0,0),ans);
    cout<<ans.data[1][1]<<endl;
    return 0;
}

　　总结一下，基本形如∑ki*f[i]（ki为定值）的递推式均可用矩阵快速幂求解，只需多加推导即可。