题解:
大胆猜想有递推式,然后BM立水之。
#include <bits/stdc++.h>
using namespace std;
const int RLEN=1<<18|1;
inline char nc() {
static char ibuf[RLEN],*ib,*ob;
(ib==ob) && (ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
return (ib==ob) ? -1 : *ib++;
}
inline int rd() {
char ch=nc(); int i=0,f=1;
while(!isdigit(ch)) {if(ch=='-')f=-1; ch=nc();}
while(isdigit(ch)) {i=(i<<1)+(i<<3)+ch-'0'; ch=nc();}
return i*f;
}
const int mod=1e9+7;
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 (long long)x*y%mod;}
inline int power(int a,int b,int rs=1) {for(;b;b>>=1,a=mul(a,a)) if(b&1) rs=mul(rs,a); return rs;}
const int N=3e3+50, M=6e3+50;
vector <int> coef;
struct poly {
vector <int> a;
poly(int d=0,int t=0) {a.resize(d+1); a[d]=t;}
inline int& operator [](const int &b) {return a[b];}
inline const int& operator [](const int &b) const {return a[b];}
inline int deg() const {return a.size()-1;}
friend inline poly operator *(const poly &a,const poly &b) {
poly c(a.deg()+b.deg(),0);
for(int i=0;i<=a.deg();i++)
for(int j=0;j<=b.deg();j++)
c[i+j]=add(c[i+j],mul(a[i],b[j]));
return c;
}
friend inline poly operator +(const poly &a,const poly &b) {
poly c(max(a.deg(),b.deg()));
for(int i=0;i<=a.deg();i++) c[i]=add(c[i],a[i]);
for(int i=0;i<=b.deg();i++) c[i]=add(c[i],b[i]);
return c;
}
inline void opt() {
int d=coef.size()-1;
for(int i=deg();i>=d;i--) if(a[i])
for(int j=1;j<=d;j++)
a[i-j]=add(a[i-j],mul(a[i],coef[j]));
a.resize(d);
}
};
namespace bm {
int L,cnt;
int a[M],fail[M],delta[M];
poly R[M];
inline vector <int> solve() {
for(int i=1;i<=L;i++) {
int d=a[i];
if(cnt) for(int j=1;j<=R[cnt].deg();++j) d=dec(d,mul(R[cnt][j],a[i-j]));
if(!d) continue;
fail[cnt]=i; delta[cnt]=d;
if(!cnt) {
++cnt;
} else {
int coef=mul(delta[cnt],power(delta[cnt-1],mod-2));
for(int j=fail[cnt-1]+1;j<i;j++) R[cnt+1].a.push_back(0);
R[cnt+1].a.push_back(coef);
for(int j=1;j<=R[cnt-1].deg();j++) R[cnt+1].a.push_back(mul(mod-coef,R[cnt-1][j]));
R[cnt+1]=R[cnt+1]+R[cnt]; ++cnt;
}
} return R[cnt].a;
}
}
long long h[N],b[N],a[N];
int main() {
bm::L=9;
h[0]=2, h[1]=3, h[2]=6;
for(int i=3;i<=bm::L+2;i++) h[i]=4*h[i-1]+17*h[i-2]-12*h[i-3]-16;
for(int i=1;i<=bm::L+1;i++)
b[i]=3*h[i+1]*h[i]+9*h[i+1]*h[i-1]+9*h[i]*h[i]+27*h[i]*h[i-1]-18*h[i+1]-126*h[i]-81*h[i-1]+192;
for(int i=1;i<=bm::L;i++) bm::a[i]=sqrt(b[i+1]+(1ll<<((i+1)*2)));
coef=bm::solve();
for(int i=(cin>>i,i);i;i--) {
long long k; cin>>k; k--;
poly f(0,1), g(1,1);
for(k--;k;k>>=1,g=g*g,g.opt())
if(k&1) f=f*g,f.opt();
int ans=0;
for(int i=0;i<=f.deg();i++)
ans=add(ans,mul(f[i],bm::a[i+1]));
cout<<ans<<'\n';
}
}