利用Cayley-Hamilton theorem 优化矩阵线性递推

平时有关线性递推的题,很多都可以利用矩阵乘法来解决。 时间复杂度一般是O(K³logn)因此对矩阵的规模限制比较大。

下面介绍一种利用利用Cayley-Hamilton theorem加速矩阵乘法的方法。

 

Cayley-Hamilton theorem：

记矩阵A的特征多项式为f(x)。 则有f(A)=0.

证明可以看 维基百科 https://en.wikipedia.org/wiki/Cayley–Hamilton_theorem#A_direct_algebraic_proof

另外我在高等代数的课本上找到了 证明(和维基百科里的第一种证明方法是一样的)

 

 

 

下面介绍几个 利用可以这个定理解决的题目：

1. project euler 258    

显然可以用矩阵乘法来做。下面讲一下怎么利用Cayley-Hamilton theorem 来优化： 详细的论述可以参考这篇。

设M为K阶矩阵,主要思想就是将$M^n$表示为 $b_0 M^0\ +\ b_1 M^1\ +\ \cdots\ b_{K-1}M^{K-1}$这样的形式.

根据Cayley-Hamilton theorem $M^K\ =\ a_0 M^0\ +\ a_1 M^1\ +\ \cdots\ a_{K-1}M^{K-1}$    

由于转移矩阵的特殊性,不难证明$a_i$恰好是线性递推公式里的系数。

假设我们已经将$M^n$表示为 $b_0 M^0\ +\ b_1 M^1\ +\ \cdots\ b_{K-1}M^{K-1}$这样的形式,不难得到$M^{n+1}$的表示法。只要将$M^n$乘个M之后得到的项中$M^K$拆成小于K次的线性组合就好了。  这样我们可以预处理出$M^0\ M^1\ \cdots\ M^{2K-2}$的表示法。

对于次数更高的, $M^{i+j}=M^i*M^j$  可以看成是两个多项式的乘法。 利用快速幂 可以在O(K²logn)的时间求出$M^n$的表示法.

 

另外有一个优化常数的trick, 可以预处理出$M^1$ $M^2$  $M^4$  $M^8$....  $M^{2^r}$这些项, 对于$M^n$只要根据二进制位相应的乘上这些项就好了。 这样做比直接做快速幂快一倍(少了一半的多项式乘法操作)。

 

参考代码：

//ans=12747994 
#include <cstdio>
#include <iostream>
#include <queue>
#include <algorithm>
#include <cstring>
#include <set>
using namespace std;

#define N 2000
typedef long long ll;

const int Mod=20092010;
int a[N],f[N<<1];
int k=2000;



//基本思想是把A^n 表示成A^0  A^1 A^2 ... A^(k-1)的线性组合
//A^(p+q)可看成两个多项式相乘,只要实现预处理出A^0  A^1 A^2 ... A^(2k-2)的多项式表示法 
//A^k可以根据特征多项式的性质得到 ,A^(n+1)次可以从A^n次得到  根据这个来预处理 
struct Poly
{
    int b[N];
}P[N<<1];


Poly operator * (const Poly &A,const Poly &B)
{
    Poly ans; memset(ans.b,0,sizeof(ans.b));
    for (int i=0;i<=2*k-2;i++)
    {
        int res=0;
        for (int j=max(0,i-k+1);j<k && j<=i;j++)
        {
            res+=1ll*A.b[j]*B.b[i-j]%Mod;
            if (res>=Mod) res-=Mod;
        }
        if (i<k) {ans.b[i]=res; continue;}
        
        //把次数大于等于k的搞成小于k 
        for (int j=0;j<k;j++)   
        {
            ans.b[j]+=1ll*res*P[i].b[j]%Mod;
            if (ans.b[j]>=Mod) ans.b[j]-=Mod;
        }
    }
    return ans;
}

Poly Power_Poly(ll p)
{
    if (p<=2*k-2) return P[p];
    
    Poly ans=P[0],A=P[1];
    for (;p;p>>=1)
    {
        if (p&1) ans=ans*A;
        A=A*A;
    }
    return ans;
}

int main()
{
    freopen("in.in","r",stdin);
    freopen("out.out","w",stdout);
    
    //f[n]=a[k-1]f[n-1]....a[0]f[n-k]
    a[0]=a[1]=1; ll n; n=1e18;
    for (int i=0;i<k;i++) P[i].b[i]=1;
    
    //P[k]=a[0]P[0]+a[1]P[1]+....a[k-1]P[k-1]
    for (int i=0;i<k;i++) P[k].b[i]=a[i]; 
    
    //Calculate P[k+1]...P[2k-2]
    //using P[n+1]=a[0]*b[k-1]+ (a[1]*b[k-1]+b[0]) + (a[2]*b[k-1]+b[1]) +...(a[k-1]*b[k-1]+b[k-2])
    for (int j=k+1;j<=2*k-2;j++)
    {
        P[j].b[0]=1ll*a[0]*P[j-1].b[k-1]%Mod;
        for (int i=1;i<k;i++)
            P[j].b[i]=(1ll*a[i]*P[j-1].b[k-1]%Mod+P[j-1].b[i-1])%Mod;
    }
    
    Poly tmp=Power_Poly(n-k+1); int ans=0;
    
    for (int i=0;i<k;i++) f[i]=1;
    for (int i=k;i<=2*k-2;i++) f[i]=(f[i-1999]+f[i-2000])%Mod;
    
    //A^n*X=b[0]*A^0*X+b[1]*A^1*X+...b[k-1]*A^(k-1)*X      A^i*X= {f[k-1+i]  f[k-2+i]...  f[0+i]}
    for (int i=0;i<k;i++)
    {
        ans+=1ll*tmp.b[i]*f[k-1+i]%Mod;
        if (ans>=Mod) ans-=Mod;
    }
    printf("%d\n",ans);
    return 0;
}

 

 

---

 

2.设，求的值。其中，和

     。

题目链接：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1229 

 

首先这个题可以用扰动法 搞出一个关于k的递推公式,在O(k²)的时间解决。

具体可以参考这篇。 虽然不是同一个式子,但是方法是一样的,扰动法在《具体数学》上也有介绍。因此本文不再赘述。

给个AC代码供参考：

 

#include <cstdio>
#include <iostream>
#include <queue>
#include <algorithm>
#include <cstring>
#include <set>
using namespace std;

#define N 2010
typedef long long ll;
const int Mod=1000000007;

int k;
ll n,r;
int f[N],inv[N];
int fac[N],fac_inv[N];

int C(int x,int y)
{
    if (y==0) return 1;
    if (y>x) return 0;
    
    int res=1ll*fac[x]*fac_inv[y]%Mod;
    return 1ll*res*fac_inv[x-y]%Mod;
}

int Power(ll a,ll p)
{
    int res=1; a%=Mod;
    for (;p;p>>=1)
    {
        if (p&1) res=1ll*res*a%Mod;
        a=a*a%Mod;
    }
    return res;
}

int Solve1()
{
    f[0]=n;  int t=n+1;
    for (int i=1;i<=k;i++)
    {
        f[i]=t=1ll*t*(n+1)%Mod;
        for (int j=0;j<i;j++)
        {
            f[i]+=Mod-1ll*C(i+1,j)*f[j]%Mod;
            if (f[i]>=Mod) f[i]-=Mod;
        }
        f[i]--; if (f[i]<0) f[i]+=Mod;
        f[i]=1ll*f[i]*inv[i+1]%Mod;
    }
    return f[k];
}


int Solve2()
{
    f[0]=Power(r,n+1)-r%Mod;
    if (f[0]<0) f[0]+=Mod;
    f[0]=1ll*f[0]*Power(r-1,Mod-2)%Mod;
    
    for (int i=1;i<=k;i++)
    {
        f[i]=1ll*Power(n+1,i)*Power(r,n+1)%Mod;
        f[i]-=r%Mod; if (f[i]<0) f[i]+=Mod;
        
        int tmp=0;
        for (int j=0;j<i;j++)
        {
            tmp+=1ll*C(i,j)*f[j]%Mod;
            if (tmp>=Mod) tmp-=Mod;
        }
        f[i]-=1ll*(r%Mod)*tmp%Mod;
        if (f[i]<0) f[i]+=Mod;
        f[i]=1ll*f[i]*Power(r-1,Mod-2)%Mod;
        //cout<<i<<" "<<f[i]<<endl;
    }
    return f[k];
}


            
int main()
{
    //freopen("in.in","r",stdin);
    //freopen("out.out","w",stdout);
    
    inv[1]=1; for (int i=2;i<N;i++) inv[i]=1ll*(Mod-Mod/i)*inv[Mod%i]%Mod;
    fac[0]=1; for (int i=1;i<N;i++) fac[i]=1ll*fac[i-1]*i%Mod;
    fac_inv[0]=1; for (int i=1;i<N;i++) fac_inv[i]=1ll*fac_inv[i-1]*inv[i]%Mod;
    
    int T; scanf("%d",&T);
    while (T--)
    {
        cin >> n >> k >> r;
        if (r==1) printf("%d\n",Solve1());
        else printf("%d\n",Solve2());
    }
    
    return 0;
}

 

 

 

本文要介绍的是利用优化后的矩阵乘法来解决本题(也许是我写的太丑,常数巨大,极限数据要跑6s,不能AC本题,但是可以拿来作为练习)

首先要知道怎么构造矩阵。

 

设$F(n,j)=n^j r^n$ 列向量 $X=[S_{n-1}\ ,\ F(n,k)\ ,\ F(n,k-1),\ \cdots\ ,\ F(n,0)]^T$

 

 

更加详细的题解可以参考这篇. 上面的推导过程的出处(矩阵实在不会用latex公式打,就只好copy了。)

我的方法和他略有不同, 我的列向量第一个元素是$S_{n-1}$ ,因此我的转移矩阵的第一行是1,1,0,0...0  其他都是一样的。

 

另外我猜测这题的这个式子和国王奇遇记那题一样,应该也是一个什么玩意乘上一个多项式,可以用多项式插值的办法来求(排行榜前面的代码跑的都很快,目测是O(K)的)。

比较懒,懒得去想了。。有兴趣的朋友可以去搞一搞？

我的代码（最后几个点TLE,用来练习 优化矩阵乘法）：

 

//http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1229
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <vector>
#include <map>
#include <cstdlib>
#include <set>
using namespace std;

#define X first
#define Y second
#define Mod 1000000007
#define N 2010
typedef long long ll;
typedef pair<int,int> pii;

int K;
ll n,r;
int fac[N],fac_inv[N],S[N];

struct Poly
{
    int cof[N];
    void print()
    {
        for (int i=0;i<K;i++) printf("%d ",cof[i]);
        printf("\n");
    }
}M[N<<1],R[62];

int a[N],b[N];

int C(int x,int y)
{
    if (y==0) return 1;
    
    int res=1ll*fac[x]*fac_inv[y]%Mod;
    return 1ll*res*fac_inv[x-y]%Mod;
}

int Power(ll x,ll p)
{
    int res=1; x%=Mod;
    for (;p;p>>=1)
    {
        if (p&1) res=1ll*res*x%Mod;
        x=x*x%Mod;
    }
    return res;
}


Poly operator * (const Poly &A,const Poly &B)
{
    Poly ans; memset(ans.cof,0,sizeof(ans.cof));
    
    for (int i=0;i<2*K-1;i++)
    {
        int res=0;
        for (int j=max(0,i-K+1);j<=i && j<K;j++)
        {
            res+=1ll*A.cof[j]*B.cof[i-j]%Mod;
            if (res>=Mod) res-=Mod;
        }
        if (i<K) {ans.cof[i]=res; continue;}
        
        for (int j=0;j<K;j++)
        {
            ans.cof[j]+=1ll*res*M[i].cof[j]%Mod;
            if (ans.cof[j]>=Mod) ans.cof[j]-=Mod;
        }
    }
    return ans;
}

Poly Poly_Power(ll p)
{
    if (p<=2*K-2) return M[p];
    
    Poly res=M[0];
    //p&(1ll<<j)  千万别忘了1后面的ll　！！！ 
    for (int j=0;j<=60;j++) if (p&(1ll<<j)) res=res*R[j];
    return res;
}


void Init()
{
    memset(a,0,sizeof(a));
    memset(b,0,sizeof(b));
    
    
    //求出特征多项式 M^(K+2)=a[0]*M^0 +  a[1]*M^1 + ... a[K+1]*M^(K+1) 
    int op;
    if (r==1)   //特征多项式f(x)= (x-1)^(K+2)
    {
        for (int i=0;i<K+2;i++)
        {
             op=(K-i)&1? -1:1;
             a[i]=op*C(K+2,i);
             a[i]=-a[i];
             if (a[i]<0) a[i]+=Mod;
        }
    }
    else   //特征多项式f(x)= (x-1) * (x-r)^(K+1)
    {
        for (int i=0;i<=K+1;i++)
        {
             op=(K+1-i)&1? -1:1;
             b[i]=1ll*op*C(K+1,i)*Power(r,K+1-i)%Mod;
             if (b[i]<0) b[i]+=Mod;
        }
        for (int i=1;i<K+2;i++) a[i]=b[i-1];
        for (int i=0;i<=K+1;i++) 
        {
            a[i]-=b[i];
            a[i]=-a[i];
            if (a[i]<0) a[i]+=Mod;
        }
    }
    
    K+=2;   //矩阵的规格
    
    
    //预处理M^0...M^(2K-2) 
    
    memset(M,0,sizeof(M));
    
    for (int i=0;i<K;i++) M[i].cof[i]=1;
    for (int i=0;i<K;i++) M[K].cof[i]=a[i];
    
    for (int i=K+1;i<=2*K-2;i++)
    {
        M[i].cof[0]=1ll*a[0]*M[i-1].cof[K-1]%Mod;
        for (int j=1;j<K;j++)
        {
            M[i].cof[j]=1ll*a[j]*M[i-1].cof[K-1]%Mod+M[i-1].cof[j-1];
            if (M[i].cof[j]>=Mod) M[i].cof[j]-=Mod;
        }
    }
     
    //预处理M^1 M^2 M^4 M^8 M^16... 
    R[0]=M[1];
    for (int i=1;i<=60;i++) R[i]=R[i-1]*R[i-1];
    
    S[0]=0;
    for (int i=1;i<K;i++) 
    {
        S[i]=1ll*Power(i,K-2)*Power(r,i)%Mod;
        S[i]+=S[i-1];
        if (S[i]>=Mod) S[i]-=Mod;
    }
}


int Solve()
{
    Poly tmp=Poly_Power(n);
    
    int ans=0,t;
    
    for (int i=0;i<K;i++)
    {
        ans+=1ll*tmp.cof[i]*S[i]%Mod;
        if (ans>=Mod) ans-=Mod;
    }
    return ans;
}

int main()
{
    //freopen("in.in","r",stdin);
    //freopen("out.out","w",stdout);    
    
    fac[0]=1;
    for (int i=1;i<N;i++) fac[i]=1ll*fac[i-1]*i%Mod;
    fac_inv[N-1]=Power(fac[N-1],Mod-2);
    for (int i=N-2;i>=0;i--) fac_inv[i]=1ll*fac_inv[i+1]*(i+1)%Mod;

    
    int T; scanf("%d",&T);
    while (T--)
    {
        cin >> n >> K >> r;
        Init();
        printf("%d\n",Solve());
    }

    
    return 0;
}

 

 

 

 

3. hdu 3483   

   福利,和上面那题几乎完全是一样的,就是 范围小了很多。可以用来测试上面未通过的代码。   坑点就是 模数 最大是2e9,  最好开个long long, 否则int范围做加法就会爆。

　参考代码：

　

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <vector>
#include <map>
#include <cstdlib>
#include <set>
using namespace std;

#define X first
#define Y second
#define N 70
typedef long long ll;
typedef pair<int,int> pii;

int K,Mod;
ll n,r;
ll S[N];
ll cob[N][N];

struct Poly
{
    ll cof[N];
    void print()
    {
        for (int i=0;i<K;i++) printf("%d ",cof[i]);
        printf("\n");
    }
}M[N<<1],R[62];

ll a[N],b[N];

ll C(int x,int y)
{
    return cob[x][y];
}

ll Power(ll x,ll p)
{
    ll res=1; x%=Mod;
    for (;p;p>>=1)
    {
        if (p&1) res=1ll*res*x%Mod;
        x=x*x%Mod;
    }
    return res;
}


Poly operator * (const Poly &A,const Poly &B)
{
    Poly ans; memset(ans.cof,0,sizeof(ans.cof));
    
    for (int i=0;i<2*K-1;i++)
    {
        ll res=0;
        for (int j=max(0,i-K+1);j<=i && j<K;j++)
        {
            res+=1ll*A.cof[j]*B.cof[i-j]%Mod;
            if (res>=Mod) res-=Mod;
        }
        if (i<K) {ans.cof[i]=res; continue;}
        
        for (int j=0;j<K;j++)
        {
            ans.cof[j]+=1ll*res*M[i].cof[j]%Mod;
            if (ans.cof[j]>=Mod) ans.cof[j]-=Mod;
        }
    }
    return ans;
}

Poly Poly_Power(ll p)
{
    if (p<=2*K-2) return M[p];
    
    Poly res=M[0];
    //p&(1ll<<j)  千万别忘了1后面的ll　！！！ 
    for (int j=0;j<=60;j++) if (p&(1ll<<j)) res=res*R[j];
    return res;
}


void Init()
{
    memset(a,0,sizeof(a));
    memset(b,0,sizeof(b));
    
    cob[0][0]=1;
    for (int i=1;i<N;i++)
    {
        cob[i][0]=1;
        for (int j=1;j<N;j++)
            cob[i][j]=(cob[i-1][j-1]+cob[i-1][j])%Mod;
    }
    
    //求出特征多项式 M^(K+2)=a[0]*M^0 +  a[1]*M^1 + ... a[K+1]*M^(K+1) 
    int op;
    if (r==1)   //特征多项式f(x)= (x-1)^(K+2)
    {
        for (int i=0;i<K+2;i++)
        {
             op=(K-i)&1? -1:1;
             a[i]=op*C(K+2,i);
             a[i]=-a[i];
             if (a[i]<0) a[i]+=Mod;
        }
    }
    else   //特征多项式f(x)= (x-1) * (x-r)^(K+1)
    {
        for (int i=0;i<=K+1;i++)
        {
             op=(K+1-i)&1? -1:1;
             b[i]=1ll*op*C(K+1,i)*Power(r,K+1-i)%Mod;
             if (b[i]<0) b[i]+=Mod;
        }
        for (int i=1;i<K+2;i++) a[i]=b[i-1];
        for (int i=0;i<=K+1;i++) 
        {
            a[i]-=b[i];
            a[i]=-a[i];
            if (a[i]<0) a[i]+=Mod;
        }
    }
    
    K+=2;   //矩阵的规格
    
    
    //预处理M^0...M^(2K-2) 
    
    memset(M,0,sizeof(M));
    
    for (int i=0;i<K;i++) M[i].cof[i]=1;
    for (int i=0;i<K;i++) M[K].cof[i]=a[i];
    
    for (int i=K+1;i<=2*K-2;i++)
    {
        M[i].cof[0]=1ll*a[0]*M[i-1].cof[K-1]%Mod;
        for (int j=1;j<K;j++)
        {
            M[i].cof[j]=1ll*a[j]*M[i-1].cof[K-1]%Mod+M[i-1].cof[j-1];
            if (M[i].cof[j]>=Mod) M[i].cof[j]-=Mod;
        }
    }
     
    //预处理M^1 M^2 M^4 M^8 M^16... 
    R[0]=M[1];
    for (int i=1;i<=60;i++) R[i]=R[i-1]*R[i-1];
    
    S[0]=0;
    for (int i=1;i<K;i++) 
    {
        S[i]=1ll*Power(i,K-2)*Power(r,i)%Mod;
        S[i]+=S[i-1];
        if (S[i]>=Mod) S[i]-=Mod;
    }
}


ll Solve()
{
    Poly tmp=Poly_Power(n);
    
    ll ans=0;
    
    for (int i=0;i<K;i++)
    {
        ans+=1ll*tmp.cof[i]*S[i]%Mod;
        if (ans>=Mod) ans-=Mod;
    }
    return ans;
}

int main()
{
    //freopen("in.in","r",stdin);
    //freopen("out.out","w",stdout);    

    int T; 
    while (true)
    {
        cin >> n >> K >> Mod; r=K;
        if (n<0) break;
        Init();
        printf("%I64d\n",Solve());
    }

    
    return 0;
}

 

　

4.codechef Dec Challenge  POWSUMS  https://www.codechef.com/problems/POWSUMS

官方题解： https://discuss.codechef.com/questions/86250/powsums-editorial

参考代码：

//https://www.codechef.com/problems/POWSUMS
//https://discuss.codechef.com/questions/86250/powsums-editorial
//https://discuss.codechef.com/questions/49614/linear-recurrence-using-cayley-hamilton-theorem
 
 
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <vector>
#include <map>
#include <cstdlib>
#include <set>
using namespace std;
 
#define X first
#define Y second
#define Mod 1000000007
#define N 310
typedef long long ll;
typedef pair<int,int> pii;
 
inline int Mul(int x,int y){return 1ll*x*y%Mod;}
inline int Add(int x,int y){return ((x+y)%Mod+Mod)%Mod;}
 
int n,Q;
int a[N],e[N],inv[N],f[N<<1];
 
struct Poly
{
    int cof[N];
    void print()
    {
        for (int i=0;i<n;i++) printf("%d ",cof[i]);
        printf("\n");
    }
}M[N<<1],tt[63];
 
Poly operator * (const Poly &A,const Poly &B)
{
    Poly ans; memset(ans.cof,0,sizeof(ans.cof));
    
    for (int i=0;i<=2*n-2;i++)
    {
        int res=0;
        for (int j=max(0,i-n+1);j<n && j<=i;j++)
        {
            res+=1ll*A.cof[j]*B.cof[i-j]%Mod;
            if (res>=Mod) res-=Mod; 
        }
        if (i<n) {ans.cof[i]=res;continue;}
        for (int j=0;j<n;j++)
        {
            ans.cof[j]+=1ll*res*M[i].cof[j]%Mod;
            if (ans.cof[j]>=Mod) ans.cof[j]-=Mod;
        }
    }
    return ans;
}
 
void Init()
{
    scanf("%d%d",&n,&Q); 
    for (int i=1;i<=n;i++) scanf("%d",&f[i]);
 
    e[0]=1; 
    for (int i=1;i<=n;i++)
    {
        int op=1; e[i]=0;
        for (int j=1;j<=i;j++,op=-op)
        {
            e[i]+=1ll*op*e[i-j]*f[j]%Mod;
            if (e[i]<0) e[i]+=Mod;
            if (e[i]>=Mod) e[i]-=Mod;
        }
        e[i]=1ll*e[i]*inv[i]%Mod;
        //cout<<i<<" "<<e[i]<<endl;
    }
    
    for (int i=n+1;i<=2*n-1;i++)
    {
        f[i]=0;  int op=1;
        for (int j=1;j<=n;j++,op=-op)
        {
            f[i]+=1ll*op*e[j]*f[i-j]%Mod;
            if (f[i]>=Mod) f[i]-=Mod;
            if (f[i]<0) f[i]+=Mod;
        }
    }
    
    for (int i=0;i<n;i++) 
        for (int j=0;j<n;j++)
            M[i].cof[j]= (i==j);
    
    //M^n= sum (a[i]*A^i)
    for (int i=n-1,op=1;i>=0;i--,op=-op)
    {
        int tmp=op*e[n-i]; 
        if (tmp<0) tmp+=Mod;
        M[n].cof[i]=a[i]=tmp;
    }
         
    //Calc linear combination form of   M^(n+1)...M^(2n-2)
    for (int i=n+1;i<=2*n-2;i++) 
    {
        M[i].cof[0]=1ll*a[0]*M[i-1].cof[n-1]%Mod;
        for (int j=1;j<n;j++)
        {
            M[i].cof[j]=1ll*a[j]*M[i-1].cof[n-1]%Mod+M[i-1].cof[j-1];
            if (M[i].cof[j]>=Mod) M[i].cof[j]-=Mod;
        }
        
    }
    
    //预处理出M^2 M^4 M^8... 随机数据可以加速很多  
    tt[0]=M[1];
    for (int i=1;i<=60;i++) tt[i]=tt[i-1]*tt[i-1];
}
 
Poly Poly_Power(ll p)
{
    if (p<=2*n-2) return M[p];
    Poly res=M[0];
    for (int j=60;j>=0;j--) if (p&(1ll<<j)) res=res*tt[j];
    return res;
}
 
int Solve(ll x)
{
    if (x<=2*n-2) return f[x];
    
    Poly tmp=Poly_Power(x-n);
    
    int res=0;
    for (int i=0;i<n;i++)
    {
        res+=1ll*tmp.cof[i]*f[n+i]%Mod;
        if (res>=Mod) res-=Mod;
    }
 
    if (res<0) res=res%Mod+Mod;
    
    return res;
}
 
int main()
{
    //freopen("in.in","r",stdin);
    //freopen("out.out","w",stdout);
    
    //calc_inv
    inv[1]=1;
    for (int i=2;i<N;i++) inv[i]=1ll*(Mod-Mod/i)*inv[Mod%i]%Mod;
    int T; ll x; 
    
    scanf("%d",&T);
    while (T--)
    {
        Init();
        while (Q--)
        {
            scanf("%lld",&x);
            printf("%d ",Solve(x));
        }
        printf("\n");
    }
 
    return 0;
}