HDU 6395 Sequence 矩阵幂,下取整,分段

题意:F[1]=A,F[2]=B,  F[i]=C*F[i-2]+D*F[i-1]+[p/i]
1<=P,n<=1e9.  求出第n项的值.

[p/i]随i不同而变化, 如果没有这一项或者这一项是一个常数的话,显然就可以用矩阵幂来推.
[p/i]下取整的值最多只有2*sqrt(p)项. (i<sqrt(p)显然结果最多sqrt(p)项, i>=sqrt(p) 结果最大为sqrt(p),则[sqrt(p),...1]).

初始前一段的两个值为F1,F2,然后用矩阵幂快速求出下一段最后两个值即可.O(sqrt(p) * log n).

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+7;
int T,A,B,C,D,P,n;
struct mat{
	ll b[3][3];
	mat(){memset(b,0,sizeof(b));}
	mat operator *(mat &tmp){
		mat ans;
		for(int i=0;i<3;i++)
			for(int j=0;j<3;j++)
				for(int k=0;k<3;k++)
					ans.b[i][j]=(ans.b[i][j]+(b[i][k]*tmp.b[k][j])%mod)%mod;
		return ans;
	}
}unit,num;
ll calc(ll i){
    ll val=P/i;
    ll l=i,r=n;
    while(l<r){
        int mid=r-(r-l)/2;
        if(P/mid==val)    l=mid;
        else if(P/mid<val)    r=mid-1;
        else    l=mid+1;
    }
    return l;
}
mat powmod(mat a,ll n){
	mat ans=unit;
	while(n){
		if(n&1)	ans=ans*a;
		n>>=1;
		a=a*a;
	}
	return ans;
}
int main()
{ 
	memset(unit.b,0,sizeof(unit.b));
	unit.b[0][0]=unit.b[1][1]=unit.b[2][2]=1;
    cin>>T;
    while(T--){
        cin>>A>>B>>C>>D>>P>>n;
        ll f1=B,f2=A;
        if(n<=2){
            printf("%d\n",n==1?f2:f1);
            continue;
        } 
	//	b={{D,1,0},{C,0,0},{0,0,1}};//{fi,fi-1,1}*B^n
		memset(num.b,0,sizeof(num.b));
      	num.b[0][0]=D,num.b[0][1]=1;
      	num.b[1][0]=C,num.b[2][2]=1;
	    for(ll i=3;i<=n;){    
            ll j=calc(i),val=P/i;// [i,i+1...j]結果都相同.都為p/i.
        //  cout<<i<<' '<<j<<' '<<val<<'\n';
          	num.b[2][0]=val;	
			mat aa=powmod(num,j-i+1);
			ll F1=((f1*aa.b[0][0])%mod+(f2*aa.b[1][0])%mod+aa.b[2][0])%mod;
            ll F2=((f1*aa.b[0][1])%mod+(f2*aa.b[1][1])%mod+aa.b[2][1])%mod;
            i=j+1;
            f1=F1,f2=F2;
        }
      	printf("%lld\n",f1);
    }
    return 0;
}