HDU - 6395 Sequence 分段矩阵快速幂 ——矩阵快速幂模板

最新推荐文章于 2021-06-02 15:00:02 发布

GreyBtfly王宝彤

最新推荐文章于 2021-06-02 15:00:02 发布

阅读量367

点赞数

分类专栏：数学类算法

本文链接：https://blog.youkuaiyun.com/GreyBtfly/article/details/81666441

版权

数学类算法专栏收录该内容

11 篇文章

订阅专栏

博客目录

题目传送门

原题

Time Limit: 4000/2000 MS (Java/Others) Memory Limit: 262144/262144 K (Java/Others)
Total Submission(s): 1159 Accepted Submission(s): 426

Problem Description

Let us define a sequence as below

Your job is simple, for each task, you should output Fn module 109+7 .

Input

The first line has only one integer T , indicates the number of tasks.

Then, for the next T lines, each line consists of 6 integers, A , B , C , D , P , n .

1≤T≤20 0≤A,B,C,D≤1e9 1≤P,n≤1e9

Sample Input

2 
3 3 2 1 3 5 
3 2 2 2 1 4

Sample Output

36
24

Source

2018 Multi-University Training Contest 7

题目大意

给定一个地推公式

其中ABCDP都是输入进来的，然后输入一个n，求Fn，由于Fn太大，对1e9+7取余后再输出。

1≤T≤20 0≤A,B,C,D≤1e9 1≤P,n≤1e9

解析

这个题n的范围是1e9，而且有20个数据，暴力肯定TLE。

看到递推公式我们想到斐波那契数列的矩阵快速幂优化方法。传送门：斐波那契递推的矩阵快速幂优化

但是这个题目难在递推公式后面加了一个p/n取整。常规是不能用矩阵快速幂做的，但这个p/n取整特殊在取整。

我们想到：P既然范围是1e9，那么P/n所有可能的结果是不是远远小于1e9呢？

所以我写了个程序跑了一下：大概跑了12s（想想如果是暴力求解肯定爆了）

#include<bits/stdc++.h>
using namespace std;
int main()
{
   int n=1e9;
   int cnt=0; //本程序计算1e9内 P/i有多少种取值
   int f=0;
   for(int i=1;i<=n;i++)
   {
       if(n/i!=f)
       {
           cnt++;
           f=n/i;
       }
   }
   cout<<cnt<<endl;
}

也就是说，虽然n的范围是1e9，但是p/n的取值个数只有6e4

也就是总有i∈[a,b]，P/i不变。那么在这一段上就可以矩阵快速幂优化求解。Fn最终是多段矩阵的乘积。

以下是公式推导过程：（~~请原谅我的书写~~）

设next(i)为使P/i发生变化的下一个i，next（i）可以用二分求解，（网上有人说可以直接一个公式,P/(P/i)求，~~反正我也没弄出来~~）

AC代码

#include<iostream>
#include<cstring>
#define regi register int
using namespace std;
int const mod=1e9+7;
//请交 G++
int const rank=3;
typedef long long ll;
class mat{
public:
   ll m[rank][rank];
   inline ll* operator[](int x){
       return this->m[x];
   }
};
mat inline operator*(mat a,mat b){
   mat ret;
   memset(&ret,0,sizeof(mat));
   for(regi i=0;i<rank;i++){
       for(regi k=0;k<rank;k++){
           for(regi j=0;j<rank;j++){
               ret[i][j]+=a[i][k]*b[k][j]%mod;
               ret[i][j]%=mod;
           }
       }
   }
   return ret;
}
mat inline mpow(mat a,int n){
   mat w,ret;
   memset(&ret,0,sizeof(mat));
   w=a;
   for(regi i=0;i<rank;i++)
       ret[i][i]=1;
   for(;n;){
       if(n&1){
           ret=ret*w;
       }
       w=w*w;
       n>>=1;
   }
   return ret;
}
mat m,temp;
ll n,p,A,B,C,D,x;
int inline next(int i){   //听说有公式可以直接算出来，但是我没考虑出来，还是老老实实二分
   regi l=i,r=n+1;
   while(l<r){
       regi m=(l+r)>>1;//二分查找第一个变化的
       if(p/m==p/i){
           //不行二分左边取不到最舒服了
           l=m+1;
       }
       else{
           r=m;
       }
   }
   return l;
}

int main(){
   m[1][0]=1;
   m[1][1]=0;
   m[1][2]=0;
   m[2][0]=0;
   m[2][1]=0;
   m[2][2]=1;
   int T;
   cin>>T;
   for(;T--;){
       cin>>A>>B>>C>>D>>p>>n;
       m[0][0]=D;
       m[0][1]=C;
       if(n==1){
           cout<<A<<endl;
           continue;
       }
       else if(n==2){
           cout<<B<<endl;
           continue;
       }

       memset(&temp,0,sizeof(mat));
       for(regi i=0;i<rank;i++)
           temp[i][i]=1;
       int i=3,t;
       for(;i<=n;)
       {
           t=next(i);
           m[0][2]=p/i;
           temp=mpow(m,t-i)*temp;
           i=t;
       }
       cout<<((temp[0][0]*B%mod+temp[0][1]*A%mod)%mod+temp[0][2])%mod<<endl;
   }
}