10.04 T4 方程+矩阵快速幂

城堡之路

(road.pas/c/cpp)

题目描述

BB要到SS的城堡去玩了。我们可以认为两座城堡位于同一条数轴上，BB的城堡坐标是0，SS的城堡坐标是N。正常情况下，BB会朝着同一个方向（即SS的城堡相对于BB的城堡的方向）走若干步之后来到SS的城堡，而且步长都为1或2。可是，今天BB在途中遇见了来自SPB的小猫KK，惊奇之下，居然有一步走反了方向！不过，BB并没有神智不清，它只有一步走反了方向，而且这一步的步长也是1或2. 同时，BB并不会路过SS的城堡而不停下来。当然，BB是在途中遇到KK的，所以它不会在自己家门口就走错方向。
 举个例子，如果SS的城堡坐标是3，那么下面两条路径是合法的：
 0->1->2->1->3
 0->1->-1->1->3
 当然，还有其它的合法路径。下面这些路径则是不合法的：
 0->-1->1->3 （BB不可能第一步就走错方向）
 0->1->3（BB一定是有一步走错方向的）
 0->2->1->0->2->3（BB只有一步是走错方向的）
 0->-1->0->3（BB每步的长度一定是1或2）
 0->1->2->4->3（BB不会越过SS的城堡再回来）
 0 -> 1 -> 2 -> 3 -> 2 -> 3（BB一旦到达了SS的城堡，就会停下来）
你现在需要帮助BB求出，它一共有多少种方法能够到达SS的城堡呢？

输入格式

 一行一个整数N，表示SS城堡的坐标。

输出格式

 一行一个整数，表示BB到SS城堡的不同路径数。由于这个数字可能很大，你只需要输出它mod 1000000007的结果。

样例输入

2

样例输出

5

样例说明

对于样例，如下5条路径是合法的：
 0->1->0->2
 0->1->-1->0->1->2
 0->1->-1->0->2
 0->1->0->1->2
 0->1->-1->1->2
数据范围与约定

对于10%的数据，N<=20.
对于70%的数据，N<=1000.
对于90%的数据，N<=1000000.
对于100%的数据，N<=10^15.

 

算法一

直接暴力搜索，搜索每步要如何走，同时判断方案的合法性。
时间复杂度：O(N!)
空间复杂度：O(N)
期望得分：10分

算法二

我们注意到，我们可以枚举如下两个内容：
1、BB在哪里返回了一步
2、BB返回的步子是1还是2
确定这两个内容之后，我们可以方便地由Fibonacci数列算出满足这两个条件的方案总数 = Sigma(F[x] * F[N-x+b])
其中，F表示Fibonacci数列，F[0]=F[1]=1；x表示我们枚举的BB的返回时坐标；b表示我们枚举的BB的步长；N表示SS城堡的坐标。
时间复杂度：O(N²)
空间复杂度：O(N)
期望得分：70分

算法三

我们发现，算法二的时间主要是浪费在了枚举上面。有没有方法可以不用枚举返回的位置呢？
答案是肯定的。我们可以把SS城堡坐标为1..5时的答案Ans[N]展开来看：
Ans[1] = 0 + 0 + 0 + 0
Ans[2] = F[1] * F[2] + F[1] * F[3] = F[1] * (F[2] + F[3]) = F[1] * F[4] + 0 + 0 + 0
Ans[3] = F[1] * F[3] + F[1] * F[4] + F[2] * F[2] + F[2] * F[3] = F[1] * F[4] +0 + F[4] * F[2] + F[3] * F[1]
Ans[4] = F[1] * F[4] + F[1] * F[5] + F[2] * F[3] + F[2] * F[4] + F[3] * F[2] + F[3] * F[3] = F[1] * F[4] + F[4] * F[2] + F[3] * F[1] + F[1] * F[4] + F[4] * F[3] + F[3]* F[2]

Ans[5] = F[1] * F[5] + F[1] * F[6] + F[2] * F[4] * F[2] * F[5] + F[3] * F[3] + F[3] * F[4] + F[4] * F[2] + F[4] * F[3] = F[1] * F[4] + F[4] * F[2] + F[3] * F[1] + F[1] * F[4] + F[4] * F[3] + F[3]* F[2] + F[1] * F[4] + F[4] * F[2] + F[3] * F[1] + F[4] * F[4] + F[3] * F[3]
可以注意到，上面的变换清楚的表示了Ans[i]的递推公式。
我们仍然用F[i]表示Fibonacci数列。其中，红色的字表示Ans[i-1]，黄色的字表示Ans[i-2]，蓝色的字表示F[4]*F[i-1]，绿色的字表示F[3]*F[i-2]。
∴ F[i] = F[i-1] + F[i-2]
Ans[i] = Ans[i-1] + Ans[i-2] + F[4] * F[i-1] + F[3] * F[i-2]
时间复杂度：O(N)
空间复杂度：O(N)
期望得分：90分

算法四

注意到N最大可能给到10¹⁵，所以需要用矩阵乘法来加速递推。
通过算法三，我们已经明确了递推式。
所以，只要构造矩阵A ={
1 1 5 3
1 0 0 0
0 0 1 1
0 0 1 0 }
其中，横向四个数分别表示Ans[i-1], Ans[i-2], F[i-1], F[i-2]
纵向四个数分别表示Ans[i], Ans[i-1], F[i], F[i-1]
然后用这个矩阵的N-2次幂去乘向量B ={
5
0
2
1 }
这四个数分别表示Ans[2], Ans[1], F[2], F[1]。
相乘得到一个4 * 1的矩阵之后，答案就是第一行第一列的那个数字了。
计算矩阵A的N-2次幂的时候，我们可以用矩阵乘法快速幂来使时间复杂度达到O(logN)
时间复杂度：O(logN)
空间复杂度：O(logN)
期望得分：100分

code：

#include<iostream>
#include<cstdio>
#include<cstring> 
using namespace std;
long long chang[5][5]={
    {0,0,0,0,0},
    {0,1,1,5,0},
    {0,1,0,3,0},
    {0,0,0,1,1},
    {0,0,0,1,0},
};
long long ni[5][5]={
    {0,0,0,0,0},
    {0,0,0,0,3},
    {0,1,0,0,5},
    {0,0,1,0,1},
    {0,0,0,1,1},
}; 
const long long mod=1e9+7;
struct matrix{
    long long a[5][5];
};
matrix mul(matrix a,matrix b){
    matrix c;
    memset(c.a,0,sizeof c.a);
    for(long long i=1;i<=4;i++){
        for(long long j=1;j<=4;j++){
            for(long long k=1;k<=4;k++){
                c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%mod;
            }
        }
    }
    return c;
}
matrix ksm(matrix a,long long b){
    matrix ans;
    memset(ans.a,0,sizeof ans.a);
    ans.a[1][1]=ans.a[2][2]=ans.a[3][3]=ans.a[4][4]=ans.a[5][5]=1;
    for(;b;b>>=1){
        if(b&1)ans=mul(ans,a);
        a=mul(a,a);
    }
    return ans;
}
void prt(matrix base){
    for(long long i=1;i<=4;i++){
        for(long long j=1;j<=4;j++){
            cout<<base.a[i][j]<<" ";
        }
        cout<<endl;
    }
}
int main(){
    long long n;
    cin>>n;
    if(n==1){
        cout<<0;
        return 0;
    }
    if(n==2){
        cout<<5;
        return 0;
    }
    if(n==3){
        cout<<18;
        return 0;
    }
    matrix base;
    memset(base.a,0,sizeof base.a);
    for(long long i=1;i<=4;i++){
        for(long long j=1;j<=4;j++){
            base.a[i][j]=chang[i][j];
        }
    }
    matrix ans;
    memset(ans.a,0,sizeof ans.a);
    ans.a[1][1]=ans.a[2][2]=ans.a[3][3]=ans.a[4][4]=ans.a[5][5]=1;
    base=ksm(base,n-2);
    long long temp=((base.a[1][3]*2+base.a[2][3])%mod+base.a[3][3]*5)%mod;
    cout<<temp;
    return 0;
}

over