HDU 6395 Sequence 矩陣快速冪+分塊
阿新 • • 發佈:2018-12-26
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≤200≤A,B,C,D≤1091≤P,n≤109
Sample Input
2
3 3 2 1 3 5
3 2 2 2 1 4
Sample Output
36
24
題意:求Fn
思路:n很大,顯然要使用矩陣快速冪,觀察有P/n向下取整,所以可以利用分塊+快速冪來計算
#include <cstdio> #include <iostream> #include <cstring> #include <algorithm> #include <map> #include <set> #include <queue> #include <cmath> using namespace std; #define Matr 4 typedef long long ll; const int mmax = 1e5 + 5; const ll mod = 1e9 + 7; int num[mmax]; struct mat//矩陣結構體,a表示內容,size大小 矩陣從1開始 { ll a[Matr][Matr],size; mat() { size=0; memset(a,0,sizeof(a)); } }; void print(mat m)//輸出矩陣資訊,debug用 { int i,j; printf("%d\n",m.size); for(i=0;i<m.size;i++) { for(j=0;j<m.size;j++)printf("%d ",m.a[i][j]); printf("\n"); } } mat multi(mat m1,mat m2)//兩個相等矩陣的乘法,對於稀疏矩陣,有0處不用運算的優化 { mat ans=mat(); ans.size=m1.size; for(int i=1;i<=m1.size;i++) for(int j=1;j<=m2.size;j++) if(m1.a[i][j])//稀疏矩陣優化 for(int k=1;k<=m1.size;k++) ans.a[i][k]=(ans.a[i][k]+m1.a[i][j]*m2.a[j][k])%mod; return ans; } mat quickmulti(mat m,mat re,int n)//二分快速冪 { mat ans=mat(); int i; for(i=1;i<=m.size;i++) for(int j=1;j<=3;j++) ans.a[i][j]=re.a[i][j]; ans.size=m.size; while(n) { if(n&1)ans=multi(m,ans); m=multi(m,m); n>>=1; } return ans; } int main() { int A,b,c,d,p,n,t; scanf("%d",&t); while(t--){ scanf("%d%d%d%d%d%d",&A,&b,&c,&d,&p,&n); mat re,k; k.a[1][1]=d; k.a[1][2]=c; k.a[1][3]=p; k.a[2][1]=1; k.a[2][2]=0; k.a[2][3]=0; k.a[3][1]=0; k.a[3][2]=0; k.a[3][3]=1; for (int i=1;i<=3;i++) re.a[i][i]=1; re.size=3; k.size=3; int mid1=0,mid2,len=0; for(int i=1;i<=p;i++){ //剛開始都出現一次,出現相同即停止 mid2=p/i; if(mid1==mid2) break; num[len++]=mid2; mid1=mid2; } for(int i=mid2-1;i>0;i--) num[len++]=i; num[len]=0; while(num[len-1]<3) len--; num[len]=2; int i; for(i=len-1;i>=0;i--) { if(n<num[i]) break; k.a[1][3]=p/num[i]; re=quickmulti(k,re,num[i]-num[i+1]); } k.a[1][3]=p/n; re=quickmulti(k,re,n-num[i+1]); int sum=(re.a[1][1]*b+re.a[1][2]*A+re.a[1][3])%mod; printf("%d\n",sum); } return 0; }