杜教BM（解决线性递推式的模板）

一、模板

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
#include<vector>
#include<map>
#include<set>
#include<cmath>
#include<queue>
#include<bitset>
#define ll long long
#define llu unsigned ll
#define pr make_pair
#define pb push_back
//#define lc (cnt<<1)
//#define rc (cnt<<1|1)
using namespace std;
const int inf=0x3f3f3f3f;
const ll lnf=0x3f3f3f3f3f3f3f3f;
const double dnf=1e15;
const int mod=998244353;
const double eps=1e-8;
const int maxn=550;

namespace linear_seq
{
    #define SZ(x) ((int)(x).size())
    typedef vector<int> VI;
    typedef pair<int,int> PII;
    const int N=10010;
    ll res[N],base[N],c[N],md[N];
    vector<int> Md;
    const ll mod=998244353;

    ll powmod(ll a,ll b)
    {
        ll res=1;
        a%=mod;
        for(;b;b>>=1)
        {
            if(b&1)
                res=res*a%mod;
            a=a*a%mod;
        }
        return res;
    }

    void mul(ll *a,ll *b,int k)
    {
        for(int i=0;i<k+k;++i) c[i]=0;
        for(int i=0;i<k;++i)
            if(a[i])for(int j=0;j<k ;++j)
                    c[i+j]=(c[i+j]+a[i]*b[j])%mod;

        for (int i=k+k-1;i>=k;i--)
            if(c[i])for(int j=0;j<(int)Md.size();++ j)
                    c[i-k+Md[j]]=(c[i-k+Md[j]]-c[i]*md[Md[j]])%mod;

        for(int i=0;i<k;++i) a[i]=c[i];
    }

    int solve(ll n,VI a,VI b)
    {
        ll ans=0,pnt=0;
        int k=SZ(a);
        for(int i=0;i<k;++i) md[k-1-i]=-a[i];
        md[k]=1;Md.clear();
        for(int i=0;i<k;++i)
            if(md[i]!=0)Md.push_back(i);

        for(int i=0;i<k;++i) res[i]=base[i]=0;
        res[0]=1;

        while ((1ll<<pnt)<=n)
            pnt++;

        for (int p=pnt;p>=0;p--)
        {
            mul(res,res,k);
            if ((n>>p)&1)
            {
                for (int i=k-1;i>=0;i--) res[i+1]=res[i];
                res[0]=0;
                for(int j=0;j<(int)Md.size();++ j)
                    res[Md[j]]=(res[Md[j]]-res[k]*md[Md[j]])%mod;
            }
        }

        for(int i=0;i<k;i++) ans=(ans+res[i]*b[i])%mod;
        if (ans<0) ans+=mod;
        return ans;
    }

    VI BM(VI s)
    {
        VI C(1,1),B(1,1);
        int L=0,m=1,b=1;
        for(int n=0;n<(int)s.size();++n)
        {
            ll d=0;
            for(int i=0;i<L+1;++ i)
                d=(d+(ll)C[i]*s[n-i])%mod;
            if(d==0) ++m;
            else if(2*L<=n)
            {
                VI T=C;
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m)C.push_back(0);
                for(int i=0;i<(int)B.size();++i) C[i+m]=(C[i+m]+c*B[i])%mod;
                L=n+1-L; B=T; b=d; m=1;
            }
            else
            {
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m)C.push_back(0);
                for(int i=0;i<(int)B.size();++ i) C[i+m]=(C[i+m]+c*B[i])%mod;
                ++m;
            }
        }
        return C;
    }

    ll gao(VI a,ll n)
    {
        VI c=BM(a);
        c.erase(c.begin());
        for(int i=0;i<(int)c.size();++i) c[i]=(mod-c[i])%mod;
        return (ll)solve(n,c,VI(a.begin(),a.begin()+SZ(c)));
    }
};

int main()
{
    int n,k;
    scanf("%d%d",&k,&n);//给定前k项求第n项

    vector<int>vc;//项数越多越好
    for(int i=1;i<=k;i++)
    {
        int x;
        scanf("%d",&x);
        vc.pb(x);
    }

    printf("%lld\n",linear_seq::gao(vc,n-1));
    return 0 ;
}

二、P5487 【模板】Berlekamp-Massey算法：

注意：这里的Pm其实是第m+1项。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
#include<vector>
#include<map>
#include<set>
#include<cmath>
#include<queue>
#include<bitset>
#define ll long long
#define llu unsigned ll
#define pr make_pair
#define pb push_back
//#define lc (cnt<<1)
//#define rc (cnt<<1|1)
using namespace std;
const int inf=0x3f3f3f3f;
const ll lnf=0x3f3f3f3f3f3f3f3f;
const double dnf=1e15;
const int mod=998244353;
const double eps=1e-8;
const int maxn=550;

namespace linear_seq
{
    #define SZ(x) ((int)(x).size())
    typedef vector<int> VI;
    typedef pair<int,int> PII;
    const int N=10010;
    ll res[N],base[N],c[N],md[N];
    vector<int> Md;
    const ll mod=998244353;

    ll powmod(ll a,ll b)
    {
        ll res=1;
        a%=mod;
        for(;b;b>>=1)
        {
            if(b&1)
                res=res*a%mod;
            a=a*a%mod;
        }
        return res;
    }

    void mul(ll *a,ll *b,int k)
    {
        for(int i=0;i<k+k;++i) c[i]=0;
        for(int i=0;i<k;++i)
            if(a[i])for(int j=0;j<k ;++j)
                    c[i+j]=(c[i+j]+a[i]*b[j])%mod;

        for (int i=k+k-1;i>=k;i--)
            if(c[i])for(int j=0;j<(int)Md.size();++ j)
                    c[i-k+Md[j]]=(c[i-k+Md[j]]-c[i]*md[Md[j]])%mod;

        for(int i=0;i<k;++i) a[i]=c[i];
    }

    int solve(ll n,VI a,VI b)
    {
        ll ans=0,pnt=0;
        int k=SZ(a);
        for(int i=0;i<k;++i) md[k-1-i]=-a[i];
        md[k]=1;Md.clear();
        for(int i=0;i<k;++i)
            if(md[i]!=0)Md.push_back(i);

        for(int i=0;i<k;++i) res[i]=base[i]=0;
        res[0]=1;

        while ((1ll<<pnt)<=n)
            pnt++;

        for (int p=pnt;p>=0;p--)
        {
            mul(res,res,k);
            if ((n>>p)&1)
            {
                for (int i=k-1;i>=0;i--) res[i+1]=res[i];
                res[0]=0;
                for(int j=0;j<(int)Md.size();++ j)
                    res[Md[j]]=(res[Md[j]]-res[k]*md[Md[j]])%mod;
            }
        }

        for(int i=0;i<k;i++) ans=(ans+res[i]*b[i])%mod;
        if (ans<0) ans+=mod;
        return ans;
    }

    VI BM(VI s)
    {
        VI C(1,1),B(1,1);
        int L=0,m=1,b=1;
        for(int n=0;n<(int)s.size();++n)
        {
            ll d=0;
            for(int i=0;i<L+1;++ i)
                d=(d+(ll)C[i]*s[n-i])%mod;
            if(d==0) ++m;
            else if(2*L<=n)
            {
                VI T=C;
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m)C.push_back(0);
                for(int i=0;i<(int)B.size();++i) C[i+m]=(C[i+m]+c*B[i])%mod;
                L=n+1-L; B=T; b=d; m=1;
            }
            else
            {
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m)C.push_back(0);
                for(int i=0;i<(int)B.size();++ i) C[i+m]=(C[i+m]+c*B[i])%mod;
                ++m;
            }
        }
        return C;
    }

    ll gao(VI a,ll n)
    {
        VI c=BM(a);
        c.erase(c.begin());
        for(int i=0;i<(int)c.size();++i) c[i]=(mod-c[i])%mod;
		
		//vector c 保存的是最短递推式的系数
        for(int i=0;i<(int)c.size();++i)
        {
            if(i!=0) putchar(' ');
            printf("%d",c[i]);
        }
        putchar('\n');

        return (ll)solve(n,c,VI(a.begin(),a.begin()+SZ(c)));
    }
};

int main()
{
    int n,m;
    scanf("%d%d",&n,&m);
    vector<int>vc;
    for(int i=1;i<=n;i++)
    {
        int x;
        scanf("%d",&x);
        vc.pb(x);
    }

    printf("%lld\n",linear_seq::gao(vc,m));
    return 0 ;
}

三、一个板子搞定模数为质数的情况与模数不为质数的情况。
（1）模数任意：
牛客：The power of Fibonacci

注意：序列是从第0项开始的。
斐波那契数列前n项m次幂和是一个m+2阶的线性递推方程。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
#include<vector>
#include<map>
#include<set>
#include<cmath>
#include<queue>
#include<bitset>
#define ll long long
#define llu unsigned ll
#define pr make_pair
#define pb push_back
//#define lc (cnt<<1)
//#define rc (cnt<<1|1)
using namespace std;
const int inf=0x3f3f3f3f;
const ll lnf=0x3f3f3f3f3f3f3f3f;
const double dnf=1e15;
const int mod=1000000000;
const double eps=1e-8;
const int maxn=550;


namespace linear_seq
{
    #define SZ(x) ((ll)(x).size())
    #define FOR(i,a,b) for(int i(a);i<=(b);++i)
    #define FOL(i,a,b) for(int i(a);i>=(b);--i)
    typedef vector<ll> VI;
    using vec = std::vector<ll>;
    const int N=10010;
    ll res[N],base[N],c[N],md[N];
    const ll mod = 1000000000;
    vector<int> Md;

    ll powmod(ll a,ll b)
    {
        ll res=1LL;
        while(b)
        {
            if(b&1) res=res*a%mod;
            a=a*a%mod;
            b>>=1;
        }
        return res;
    }

    void mul(ll *a,ll *b,int k)
    {
        for(int i=0;i<k+k;++i) c[i]=0;
        for(int i=0;i<k;++i)
            if(a[i])for(int j=0;j<k ;++j)
                    c[i+j]=(c[i+j]+a[i]*b[j])%mod;

        for (int i=k+k-1;i>=k;i--)
            if(c[i])for(int j=0;j<(int)Md.size();++ j)
                    c[i-k+Md[j]]=(c[i-k+Md[j]]-c[i]*md[Md[j]])%mod;

        for(int i=0;i<k;++i) a[i]=c[i];
    }
    int solve(ll n,VI a,VI b)
    {   // a 系数 b 初值 b[n+1]=a[0]*b[n]+...
        ll ans=0,pnt=0;
        int k=SZ(a);
        for(int i=0;i<k;++i) md[k-1-i]=-a[i];
        md[k]=1;Md.clear();
        for(int i=0;i<k;++i) if(md[i]!=0)Md.push_back(i);

        for(int i=0;i<k;++i) res[i]=base[i]=0;
        res[0]=1;

        while ((1ll<<pnt)<=n) pnt++;

        for (int p=pnt;p>=0;p--)
        {
            mul(res,res,k);
            if ((n>>p)&1)
            {
                for (int i=k-1;i>=0;i--) res[i+1]=res[i];
                res[0]=0;
                for(int j=0;j<(int)Md.size();++j)
                    res[Md[j]]=(res[Md[j]]-res[k]*md[Md[j]])%mod;
            }
        }
        for(int i=0;i<k;i++) ans=(ans+res[i]*b[i])%mod;
        if (ans<0) ans+=mod;
        return ans;
    }

    VI BM(VI s)
    {
        VI C(1,1),B(1,1);
        int L=0,m=1,b=1;
        for(int n=0;n<(int)s.size();++n)
        {
            ll d=0;
            for(int i=0;i<L+1;++ i) d=(d+(ll)C[i]*s[n-i])%mod;
            if(d==0) ++m;
            else if(2*L<=n)
            {
                VI T=C;
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m)C.push_back(0);
                for(int i=0;i<(int)B.size();++i) C[i+m]=(C[i+m]+c*B[i])%mod;
                L=n+1-L; B=T; b=d; m=1;
            }
            else
            {
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m)C.push_back(0);
                for(int i=0;i<(int)B.size();++ i) C[i+m]=(C[i+m]+c*B[i])%mod;
                ++m;
            }
        }
        return C;
    }

    static void extand(vec &a, size_t d, ll value=0)
    {
        if (d<=a.size()) return;
        a.resize(d, value);
    }

    static void exgcd(ll a, ll b, ll &g, ll &x, ll &y)
    {
        if (!b) x=1,y=0,g=a;
        else
        {
            exgcd(b,a%b,g,y,x);
            y-=x*(a/b);
        }
    }

    static ll crt(const vec &c, const vec &m)
    {
        int n=c.size();
        ll M=1,ans=0;
        for(int i=0;i<n;++i) M*=m[i];
        for(int i=0;i<n;++i)
        {
            ll x,y,g,tm=M/m[i];
            exgcd(tm,m[i],g,x,y);
            ans=(ans+tm*x*c[i]%M)%M;
        }
        return (ans+M)%M;
    }

    static vec ReedsSloane(const vec &s, ll mod)
    {
        auto inverse=[](ll a, ll m)
        {
            ll d,x,y;
            exgcd(a,m,d,x,y);
            return d==1?(x%m+m)%m:-1;
        };

        auto L=[](const vec &a, const vec &b)
        {
            int da=(a.size()>1||(a.size()==1&&a[0]))?a.size()-1:-1000;
            int db=(b.size()>1||(b.size()==1&&b[0]))?b.size()-1:-1000;
            return std::max(da,db+1);
        };

        auto prime_power = [&](const vec &s, ll mod, ll p, ll e)
        {
            // linear feedback shift register mod p^e, p is prime
            std::vector<vec> a(e),b(e),an(e),bn(e),ao(e),bo(e);
            vec t(e),u(e),r(e),to(e,1),uo(e),pw(e + 1);
            pw[0]=1;
            for(int i=pw[0]=1;i<=e;++i) pw[i]=pw[i-1]*p;
            for(ll i=0;i<e;++i)
            {
                a[i]={pw[i]},an[i]={pw[i]};
                b[i]={0},bn[i]={s[0]*pw[i]%mod};
                t[i]=s[0]*pw[i]%mod;
                if (t[i]==0)
                {
                    t[i]=1,u[i]=e;
                }
                else
                {
                    for(u[i]=0;t[i]%p==0;t[i]/=p,++u[i]);
                }
            }

            for(size_t k=1;k<s.size();++k)
            {
                for(int g=0;g<e;++g)
                {
                    if (L(an[g],bn[g])>L(a[g], b[g]))
                    {
                        ao[g]=a[e-1-u[g]];
                        bo[g]=b[e-1-u[g]];
                        to[g]=t[e-1-u[g]];
                        uo[g]=u[e-1-u[g]];
                        r[g]=k-1;
                    }
                }
                a=an,b=bn;
                for(int o=0;o<e;++o)
                {
                    ll d=0;
                    for(size_t i=0;i<a[o].size()&&i<=k;++i)
                    {
                        d=(d+a[o][i]*s[k-i])% mod;
                    }
                    if(d==0)
                    {
                        t[o]=1,u[o]=e;
                    }
                    else
                    {
                        for(u[o]=0,t[o]=d;t[o]%p==0;t[o]/=p,++u[o]);
                        int g=e-1-u[o];
                        if(L(a[g],b[g])==0)
                        {
                            extand(bn[o],k+1);
                            bn[o][k]=(bn[o][k]+d)%mod;
                        }
                        else
                        {
                            ll coef=t[o]*inverse(to[g],mod)%mod*pw[u[o]-uo[g]]%mod;
                            int m=k-r[g];
                            extand(an[o],ao[g].size()+m);
                            extand(bn[o],bo[g].size()+m);
                            for(size_t i=0;i<ao[g].size();++i)
                            {
                                an[o][i+m]-=coef*ao[g][i]%mod;
                                if(an[o][i+m]<0) an[o][i+m]+=mod;
                            }
                            while(an[o].size()&&an[o].back()==0) an[o].pop_back();
                            for(size_t i=0;i<bo[g].size();++i)
                            {
                                bn[o][i+m]-=coef*bo[g][i]%mod;
                                if(bn[o][i + m]<0) bn[o][i+m]-=mod;
                            }
                            while(bn[o].size()&&bn[o].back()==0) bn[o].pop_back();
                        }
                    }
                }
            }
            return std::make_pair(an[0],bn[0]);
        };
        std::vector<std::tuple<ll,ll,int> >fac;
        for(ll i=2;i*i<=mod;++i)
            if(mod%i==0)
            {
                ll cnt=0,pw=1;
                while(mod%i==0) mod/=i,++cnt,pw*=i;
                fac.emplace_back(pw,i,cnt);
            }
        if(mod>1) fac.emplace_back(mod,mod,1);
        std::vector<vec> as;
        size_t n=0;
        for (auto &&x:fac)
        {
            ll mod,p,e;
            vec a,b;
            std::tie(mod,p,e)=x;
            auto ss=s;
            for(auto &&x:ss) x%=mod;
            std::tie(a,b)=prime_power(ss,mod,p,e);
            as.emplace_back(a);
            n=std::max(n,a.size());
        }
        vec a(n),c(as.size()),m(as.size());
        for(size_t i=0;i<n;++i)
        {
            for(size_t j=0;j<as.size();++j)
            {
                m[j]=std::get<0>(fac[j]);
                c[j]=i<as[j].size()?as[j][i]:0;
            }
            a[i]=crt(c,m);
        }
        return a;
    }
    ll gao(VI a,ll n,ll mod,bool prime=true)
    {
        VI c;
        if(prime) c=BM(a);
        else c=ReedsSloane(a,mod);
        c.erase(c.begin());
        for(int i=0;i<SZ(c);++i) c[i]=(mod-c[i])%mod;
        return solve(n,c,VI(a.begin(),a.begin()+SZ(c)));
    }
};


int main()
{
    ll n, m;
    scanf("%lld%lld",&n,&m);

    vector<ll>f;
    f.pb(0),f.pb(1);
    for(int i=2;i<=2000;i++) f.pb((f[i-1]+f[i-2])%mod);
    for(int i=1;i<=2000;i++) f[i]=linear_seq::powmod(f[i],m);
    for(int i=1;i<=2000;i++) f[i]=(f[i-1]+f[i])%mod;
    printf("%lld\n",linear_seq::gao(f,n,mod,false));//false表示mod为非质数。
    return 0;
}

（2）模数为质数：
牛客： 斐波那契和

注意：从第1项开始。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
#include<vector>
#include<map>
#include<set>
#include<cmath>
#include<queue>
#include<bitset>
#define ll long long
#define llu unsigned ll
#define pr make_pair
#define pb push_back
//#define lc (cnt<<1)
//#define rc (cnt<<1|1)
using namespace std;
const int inf=0x3f3f3f3f;
const ll lnf=0x3f3f3f3f3f3f3f3f;
const double dnf=1e15;
const int mod=998244353;
const double eps=1e-8;
const int maxn=550;


namespace linear_seq
{
    #define SZ(x) ((ll)(x).size())
    #define FOR(i,a,b) for(int i(a);i<=(b);++i)
    #define FOL(i,a,b) for(int i(a);i>=(b);--i)
    typedef vector<ll> VI;
    using vec = std::vector<ll>;
    const int N=10010;
    ll res[N],base[N],c[N],md[N];
    const ll mod = ::mod;
    vector<int> Md;

    ll powmod(ll a,ll b)
    {
        ll res=1LL;
        while(b)
        {
            if(b&1) res=res*a%mod;
            a=a*a%mod;
            b>>=1;
        }
        return res;
    }

    void mul(ll *a,ll *b,int k)
    {
        for(int i=0;i<k+k;++i) c[i]=0;
        for(int i=0;i<k;++i)
            if(a[i])for(int j=0;j<k ;++j)
                    c[i+j]=(c[i+j]+a[i]*b[j])%mod;

        for (int i=k+k-1;i>=k;i--)
            if(c[i])for(int j=0;j<(int)Md.size();++ j)
                    c[i-k+Md[j]]=(c[i-k+Md[j]]-c[i]*md[Md[j]])%mod;

        for(int i=0;i<k;++i) a[i]=c[i];
    }
    int solve(ll n,VI a,VI b)
    {   // a 系数 b 初值 b[n+1]=a[0]*b[n]+...
        ll ans=0,pnt=0;
        int k=SZ(a);
        for(int i=0;i<k;++i) md[k-1-i]=-a[i];
        md[k]=1;Md.clear();
        for(int i=0;i<k;++i) if(md[i]!=0)Md.push_back(i);

        for(int i=0;i<k;++i) res[i]=base[i]=0;
        res[0]=1;

        while ((1ll<<pnt)<=n) pnt++;

        for (int p=pnt;p>=0;p--)
        {
            mul(res,res,k);
            if ((n>>p)&1)
            {
                for (int i=k-1;i>=0;i--) res[i+1]=res[i];
                res[0]=0;
                for(int j=0;j<(int)Md.size();++j)
                    res[Md[j]]=(res[Md[j]]-res[k]*md[Md[j]])%mod;
            }
        }
        for(int i=0;i<k;i++) ans=(ans+res[i]*b[i])%mod;
        if (ans<0) ans+=mod;
        return ans;
    }

    VI BM(VI s)
    {
        VI C(1,1),B(1,1);
        int L=0,m=1,b=1;
        for(int n=0;n<(int)s.size();++n)
        {
            ll d=0;
            for(int i=0;i<L+1;++ i) d=(d+(ll)C[i]*s[n-i])%mod;
            if(d==0) ++m;
            else if(2*L<=n)
            {
                VI T=C;
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m)C.push_back(0);
                for(int i=0;i<(int)B.size();++i) C[i+m]=(C[i+m]+c*B[i])%mod;
                L=n+1-L; B=T; b=d; m=1;
            }
            else
            {
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m)C.push_back(0);
                for(int i=0;i<(int)B.size();++ i) C[i+m]=(C[i+m]+c*B[i])%mod;
                ++m;
            }
        }
        return C;
    }

    static void extand(vec &a, size_t d, ll value=0)
    {
        if (d<=a.size()) return;
        a.resize(d, value);
    }

    static void exgcd(ll a, ll b, ll &g, ll &x, ll &y)
    {
        if (!b) x=1,y=0,g=a;
        else
        {
            exgcd(b,a%b,g,y,x);
            y-=x*(a/b);
        }
    }

    static ll crt(const vec &c, const vec &m)
    {
        int n=c.size();
        ll M=1,ans=0;
        for(int i=0;i<n;++i) M*=m[i];
        for(int i=0;i<n;++i)
        {
            ll x,y,g,tm=M/m[i];
            exgcd(tm,m[i],g,x,y);
            ans=(ans+tm*x*c[i]%M)%M;
        }
        return (ans+M)%M;
    }

    static vec ReedsSloane(const vec &s, ll mod)
    {
        auto inverse=[](ll a, ll m)
        {
            ll d,x,y;
            exgcd(a,m,d,x,y);
            return d==1?(x%m+m)%m:-1;
        };

        auto L=[](const vec &a, const vec &b)
        {
            int da=(a.size()>1||(a.size()==1&&a[0]))?a.size()-1:-1000;
            int db=(b.size()>1||(b.size()==1&&b[0]))?b.size()-1:-1000;
            return std::max(da,db+1);
        };

        auto prime_power = [&](const vec &s, ll mod, ll p, ll e)
        {
            // linear feedback shift register mod p^e, p is prime
            std::vector<vec> a(e),b(e),an(e),bn(e),ao(e),bo(e);
            vec t(e),u(e),r(e),to(e,1),uo(e),pw(e + 1);
            pw[0]=1;
            for(int i=pw[0]=1;i<=e;++i) pw[i]=pw[i-1]*p;
            for(ll i=0;i<e;++i)
            {
                a[i]={pw[i]},an[i]={pw[i]};
                b[i]={0},bn[i]={s[0]*pw[i]%mod};
                t[i]=s[0]*pw[i]%mod;
                if (t[i]==0)
                {
                    t[i]=1,u[i]=e;
                }
                else
                {
                    for(u[i]=0;t[i]%p==0;t[i]/=p,++u[i]);
                }
            }

            for(size_t k=1;k<s.size();++k)
            {
                for(int g=0;g<e;++g)
                {
                    if (L(an[g],bn[g])>L(a[g], b[g]))
                    {
                        ao[g]=a[e-1-u[g]];
                        bo[g]=b[e-1-u[g]];
                        to[g]=t[e-1-u[g]];
                        uo[g]=u[e-1-u[g]];
                        r[g]=k-1;
                    }
                }
                a=an,b=bn;
                for(int o=0;o<e;++o)
                {
                    ll d=0;
                    for(size_t i=0;i<a[o].size()&&i<=k;++i)
                    {
                        d=(d+a[o][i]*s[k-i])% mod;
                    }
                    if(d==0)
                    {
                        t[o]=1,u[o]=e;
                    }
                    else
                    {
                        for(u[o]=0,t[o]=d;t[o]%p==0;t[o]/=p,++u[o]);
                        int g=e-1-u[o];
                        if(L(a[g],b[g])==0)
                        {
                            extand(bn[o],k+1);
                            bn[o][k]=(bn[o][k]+d)%mod;
                        }
                        else
                        {
                            ll coef=t[o]*inverse(to[g],mod)%mod*pw[u[o]-uo[g]]%mod;
                            int m=k-r[g];
                            extand(an[o],ao[g].size()+m);
                            extand(bn[o],bo[g].size()+m);
                            for(size_t i=0;i<ao[g].size();++i)
                            {
                                an[o][i+m]-=coef*ao[g][i]%mod;
                                if(an[o][i+m]<0) an[o][i+m]+=mod;
                            }
                            while(an[o].size()&&an[o].back()==0) an[o].pop_back();
                            for(size_t i=0;i<bo[g].size();++i)
                            {
                                bn[o][i+m]-=coef*bo[g][i]%mod;
                                if(bn[o][i + m]<0) bn[o][i+m]-=mod;
                            }
                            while(bn[o].size()&&bn[o].back()==0) bn[o].pop_back();
                        }
                    }
                }
            }
            return std::make_pair(an[0],bn[0]);
        };
        std::vector<std::tuple<ll,ll,int> >fac;
        for(ll i=2;i*i<=mod;++i)
            if(mod%i==0)
            {
                ll cnt=0,pw=1;
                while(mod%i==0) mod/=i,++cnt,pw*=i;
                fac.emplace_back(pw,i,cnt);
            }
        if(mod>1) fac.emplace_back(mod,mod,1);
        std::vector<vec> as;
        size_t n=0;
        for (auto &&x:fac)
        {
            ll mod,p,e;
            vec a,b;
            std::tie(mod,p,e)=x;
            auto ss=s;
            for(auto &&x:ss) x%=mod;
            std::tie(a,b)=prime_power(ss,mod,p,e);
            as.emplace_back(a);
            n=std::max(n,a.size());
        }
        vec a(n),c(as.size()),m(as.size());
        for(size_t i=0;i<n;++i)
        {
            for(size_t j=0;j<as.size();++j)
            {
                m[j]=std::get<0>(fac[j]);
                c[j]=i<as[j].size()?as[j][i]:0;
            }
            a[i]=crt(c,m);
        }
        return a;
    }
    ll gao(VI a,ll n,ll mod,bool prime=true)
    {
        VI c;
        if(prime) c=BM(a);
        else c=ReedsSloane(a,mod);
        c.erase(c.begin());
        for(int i=0;i<SZ(c);++i) c[i]=(mod-c[i])%mod;
        return solve(n,c,VI(a.begin(),a.begin()+SZ(c)));
    }
};

ll f[maxn];
ll mypow(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}

int main()
{
    ll n,k;
    scanf("%lld%lld",&n,&k);
    f[1]=f[2]=1;
    for(int i=3;i<maxn;i++)
        f[i]=(f[i-1]+f[i-2])%mod;
    ll ans=0;
    vector<ll>vc;
    for(int i=1;i<=500;i++)
    {
        ans=(ans+mypow(i,k)*f[i])%mod;
        vc.pb(ans);
    }
    //质数的时候false true均可
    //因为false可以处理更一般的情况

    //printf("%lld\n",linear_seq::gao(vc,n-1,mod,false));
    printf("%lld\n",linear_seq::gao(vc,n-1,mod,true));
    return 0;
}