Week14 矩阵快速幂初阶（HDU - 1757 ）

矩阵快速幂简述

矩阵快速幂基础知识

在学习矩阵快速幂之前，要了解两个基础知识：
1、矩阵的乘法
2、整数的快速幂原理
其中矩阵乘法是线性代数的知识，这里不做概述；
值得一说的是整数的快速幂原理：
整数快速幂要解决的问题正如字面意思：某个整数的n次幂。
在计算一个幂级数a^x的结果时，如果使用暴力的方法，显然需要乘n次，在n很大时浪费掉大量时间；且还存在另外一个问题：如果计算的是a^x%p，如果a^x可能会直接超出整数范围导致溢出。
为了解决这个问题，引入一种快速幂算法：
1、快速幂的计算结果初始化为1
2、引入快速幂基数初始化为a，进行运算的幂次为x
按照如下的算法来进行计算：

while (x>0)
    {
        if(x & 1) ret=ret*a;
        a=a*a;
        x>>=1;
    }

计算过程中使用到了位运算，如果你对位运算不熟悉，可以先看这篇了解位运算的原理和运算定理：
位运算专讲
快速幂其实利用的是保证快速幂基数的值与a的每一位上权值相对，即

1位  2位     3位
a 	a的2次幂  a的4次幂

可以看出来，位数-幂次的对应关系是
幂次=2^位数
正是因为这种对应关系，利用对a不断的平方使其等于权值，对x不断地右移并与1做与运算，如果某位为1说明结果需要乘上这位上的权值。利用这种思路一，可以将求幂级数的问题优化到O（logn）级别

了解了常数的快速幂，接下来我们将矩阵乘法运算和常数的快速幂进行组合，将常数快速幂中的底数修改为要进行幂次计算的矩阵，将乘法重载为矩阵乘法就可以完成矩阵的快速幂，矩阵的快速幂模板如下：

struct Matrix
{
    int x[N][N];
    Matrix operator*(const Matrix& t) const 
    {
        Matrix res;
        for(int i=0;i<N;i++)
        {
            for(int j=0;j<N;j++)
            {
                res.x[i][j]=0;
                for(int k=0;k<N;k++)
                {
                    res.x[i][j]+=x[i][k]*t.x[k][j] % p;
                    res.x[i][j]%=p;
                }
            }
        }
        return res;
    }
    Matrix()
    {
        memset(x,0,sizeof(x));
    }
    Matrix(const Matrix & t)
    {
        memcpy(x,t.x,sizeof(x));
    }
};
Matrix quick_pow(Matrix a,int x)
{
    Matrix ret;
    memset(ret.x,0,sizeof(ret.x));
    for(int i=0;i<N;i++)
    {
        ret.x[i][i]=1;
    }
    while (x>0)
    {
        if(x & 1) ret=ret*a;
        a=a*a;
        x>>=1;
    }
    return ret;
}

有了矩阵快速幂，我们只是拥有了工具，还需要知道应该如何使用矩阵快速幂这个强大的算法：

矩阵快速幂----优化线性DP

在线性DP中，我们经常会求出一个线性递推式，例如：
𝑓𝑛=𝑎1∗𝑓𝑛−1+𝑎2∗𝑓𝑛−2+⋯+𝑎k∗𝑓𝑛−𝑘+𝑏
通常来说我们使用枚举递推求解，要使用O（n*枚举内部运算时间），在某些题目中如果转换次数很多，就会导致TLE，这里就是矩阵快速幂的最优应用场合了。这里为了方便讲述，以求自然数幂和为例：

由于该表达式仍可以展开，我们可以将其展开为这样的形式：

根据二项式展开定理（线性代数内容）：可以将上述的表达式展开为
这是一个头矩阵的第一行为Sn和Sn-1的表达式，如果将等式右继续化简直至每一项都已知，可以将等式化简为一个
A=a^n-kB的形式，其中a和B均已知，因为要求解矩阵幂次调用矩阵快速幂模板即可解决。计算出矩阵后我们会发现，结果矩阵A的第一行即为结果。
但是实际上我们只用了O（logn）的时间就完成了之前O（n）时间完成的任务，这就是矩阵快速幂的妙处。
这只是矩阵快速幂的简单应用，矩阵快速幂的使用博大精深，我们会分为三个专题来仔细分析各种类型的矩阵快速幂问题。

题目概述

Q老师 对数列有一种非同一般的热爱，尤其是优美的斐波那契数列。
这一天，Q老师 为了增强大家对于斐波那契数列的理解，决定在斐波那契的基础上创建一个新的数列 f(x) 来考一考大家。数列 f(x) 定义如下：
当 x < 10 时，f(x) = x；
当 x ≥ 10 时，f(x) = a0 * f(x-1) + a1 * f(x-2) + a2 * f(x-3) + …… + a9 * f(x-10)，ai 只能为 0 或 1。
Q老师 将给定 a0～a9，以及两个正整数 k m，询问 f(k) % m 的数值大小。
聪明的你能通过 Q老师 的考验吗？

输入样例

输出文件包含多组测试用例，每组测试用例格式如下：
第一行给定两个正整数 k m。（k < 2e9, m < 1e5）
第二行给定十个整数，分别表示 a0～a9。

10 9999
1 1 1 1 1 1 1 1 1 1
20 500
1 0 1 0 1 0 1 0 1 0

输出样例

对于每一组测试用例输出一行，表示 f(k) % m 的数值大小。

45
104

思路概述

一道比较简单的矩阵快速幂板子题。根据题目给出的递推式我们可以发现，前10个数的数值是完全已知的，而>10的部分每一个答案的推导都需要前10个数进行递推。于是我们就可以显而易见的是：该题目可以直接套用快速幂的板子，要注意的是转移矩阵是
$\left[\begin{array}{cccccccccc}a0& a1& a2& a3& a4& a5& a6& a7& a8& a9\\ & & & & & & & & & \\ 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ & & & & & & & & & \\ 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ & & & & & & & & & \\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ & & & & & & & & & \\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ & & & & & & & & & \\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ & & & & & & & & & \\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ & & & & & & & & & \\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ & & & & & & & & & \\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ & & & & & & & & & \\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\\ & & & & & & & & & \end{array}\right]$
结果矩阵是：
$\left[\begin{array}{c}f(x)\\ \\ f(x-1)\\ \\ f(x-2)\\ \\ f(x-3)\\ \\ f(x-4)\\ \\ f(x-5)\\ \\ f(x-6)\\ \\ f(x-7)\\ \\ f(x-8)\\ \\ f(x-9)\\ \end{array}\right]$
我们可以看到结果矩阵的第一项是$f(x)$ ，所以在查询值是使用该结果即可。

实现源码（C++）

#include<iostream>
#include<stdio.h>
#include<string.h>

using namespace std;
const int N=10;
int p;
int k;
int ori[N];
int ans[N];

struct Matrix
{
    int x[N][N];
    Matrix operator*(const Matrix& t) const 
    {
        Matrix res;
        for(int i=0;i<N;i++)
        {
            for(int j=0;j<N;j++)
            {
                res.x[i][j]=0;
                for(int k=0;k<N;k++)
                {
                    res.x[i][j]+=x[i][k]*t.x[k][j] % p;
                    res.x[i][j]%=p;
                }
            }
        }
        return res;
    }
    Matrix()
    {
        memset(x,0,sizeof(x));
    }
    Matrix(const Matrix & t)
    {
        memcpy(x,t.x,sizeof(x));
    }
};
Matrix quick_pow(Matrix a,int x)
{
    Matrix ret;
    memset(ret.x,0,sizeof(ret.x));
    for(int i=0;i<N;i++)
    {
        ret.x[i][i]=1;
    }
    while (x>0)
    {
        if(x & 1) ret=ret*a;
        a=a*a;
        x>>=1;
    }
    return ret;
}

Matrix rix;

int main()
{
    while(scanf("%d %d",&k,&p)!=EOF)
    {
    	for(int i=0;i<N;i++)
        	for(int j=0;j<N;j++)
        		rix.x[i][i]=0;
        for(int i=0;i<N;i++)
        ori[9-i]=i;
        for(int i=0;i<N;i++)
        scanf("%d",&rix.x[0][i]);
        for(int i=1;i<N;i++)
        rix.x[i][i-1]=1;
        if(k<10)
        {
            printf("%d\n",k);
            continue;
        }
        Matrix res=quick_pow(rix,k-9);
        for(int i=0;i<N;i++)
        {
            ans[i]=0;
            for(int j=0;j<N;j++)
            {
                ans[i]+=res.x[i][j]*ori[j] % p;
                ans[i]%=p;
            }
        }
        printf("%d\n",ans[0]);
    }
    return 0;
}