HDU6172 Array Challenge

传送门
佛辣。
推出前几项后可以用$BM$算法得出$\lfloor \sqrt{a_n}\rfloor$的最短线性递推式。
注意是直接算出来不取模，然后开根，并不是$Cipolla$。嗯就是这个东西。。很神奇。
它居然可以线性递推。
所以要勇敢大胆打表蒙结论啊。

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const
typedef std::vector<int> poly;
cs int mod=1e9+7,N=13;

inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
inline int mul(int x,int y){return 1ll*x*y%mod;}

inline void Add(int &x,int y){x=(x+y>=mod)?x+y-mod:x+y;}
inline void Dec(int &x,int y){x=(x-y<0)?x-y+mod:x-y;}
inline void Mul(int &x,int y){x=1ll*x*y%mod;}

inline int quickpow(int a,int b,int ret=1){for(;b;b>>=1,Mul(a,a))if(b&1)Mul(ret,a);return ret;}
inline int inv(int x){return quickpow(x,mod-2);}

inline void print(cs poly &A){
	for(int re i=0;i<A.size();++i)
		printf("%lld ",A[i]);
	puts("");
}

inline poly BM(poly A){
	int delta[N],fail[2],cur=0;poly ans[2];
	for(int re i=0;i<A.size();++i){
		delta[i]=A[i];
		for(int re j=0;j<ans[cur].size();++j)
			Dec(delta[i],mul(ans[cur][j],A[i-j-1]));
		if(!delta[i]) continue;fail[cur]=i;
		if(ans[cur].empty()){cur^=1,ans[cur].resize(i+1);continue;}
		int k=mul(delta[i],inv(delta[fail[cur^1]]));poly now;
		now.resize(i-fail[cur^1]-1),now.push_back(k);
		for(int re j=0;j<ans[cur^1].size();++j)
			now.push_back(mul(dec(0,k),ans[cur^1][j]));
		if(now.size()<ans[cur].size()) now.resize(ans[cur].size());
		for(int re j=0;j<ans[cur].size();++j) Add(now[j],ans[cur][j]);
		cur^=1,ans[cur]=now;
	}return ans[cur];
}

ll H[N],B[N];
poly A;

struct matrix{
	int a[3][3];
	matrix(int t=0){
		memset(a,0,sizeof(a));
		a[1][1]=a[2][2]=t;
	}
	friend inline matrix operator*(cs matrix &a,cs matrix &b){
		matrix ret;
		for(int re i=1;i<=2;++i)
			for(int re k=1;k<=2;++k)
				for(int re j=1;j<=2;++j)
					Add(ret.a[i][j],mul(a.a[i][k],b.a[k][j]));
		return ret;
	}
	friend inline matrix operator^(matrix a,ll b){
		matrix ret(1);
		for(;b;b>>=1,a=a*a) if(b&1) ret=ret*a;
		return ret;
	}
	inline void print(){
		for(int re i=1;i<=2;++i){
			for(int re j=1;j<=2;++j)
				printf("%d ",a[i][j]);
			puts("");
		}puts("");
	}
}T,G;
int Case;ll x;
int main(){
	H[0]=2,H[1]=3,H[2]=6;
	for(int re i=3;i<N;++i) H[i]=H[i-1]*4+H[i-2]*17-H[i-3]*12-16;
	for(int re i=1;i<N-1;++i) B[i]=H[i+1]*H[i]*3+H[i+1]*H[i-1]*9+H[i]*H[i]*9+H[i]*H[i-1]*27-H[i+1]*18-H[i]*126-H[i-1]*81+192;
	for(int re i=2;i<N-1;++i) A.push_back((int)(((ll)sqrt(B[i]+(ll)pow(4,i)))%mod));
	//print(A);poly trans=BM(A);print(trans);
	T.a[1][1]=7,T.a[1][2]=mod-4,T.a[2][1]=1,T.a[2][2]=0;
	G.a[1][1]=A[1],G.a[2][1]=A[0],G.a[1][2]=G.a[2][2]=0;
	scanf("%d",&Case);
	while(Case--){
		scanf("%lld",&x);
		printf("%d\n",((T^(x-2))*G).a[2][1]);
	}
}