codeforces 1117 D. Magic Gems(矩阵快速幂||杜教BM)

题目：
http://codeforces.com/contest/1117/problem/D

题意：有一种魔法宝石，每一个可以分成m个普通宝石，
问有多少种方法得到总共n个宝石。

思路：
首先动态规划思路可以得到转移方程 dp表示i个宝石的方案数。

dp[n]=dp[n-m]+dp[n-1]
//  分解最后一个宝石+不分解最后一个宝石

数据范围过大不能直接暴力推

方法一：矩阵快速幂优化
将递推转化为矩阵乘法
(这里转化矩阵还是有点懵…)
我的理解找到边界情况 即为 dp[1] …dp[m] 的值。后面所有的dp值均可以由此的推出，那么以上就为 初始矩阵A。

那么相邻一个的dp值关系转化可以用下图表示，那么中间的矩阵即为转化矩阵，记为S；

最终末矩阵记为B
B=S^(n-m)A;
最终答案值即是B矩阵中位于[1][1]的元素；

#include<bits/stdc++.h>
#define fi first
#define se second
#define FOR(a) for(int i=0;i<a;i++)
#define sc(a) scanf("%d",&a)
#define show(a) cout<<a<<endl;
#define show2(a,b) cout<<a<<" "<<b<<endl;
#define show3(a,b,c) cout<<a<<" "<<b<<" "<<c<<endl;
using namespace std;

typedef long long ll;
typedef pair<int, int> P;
typedef pair<P, int> LP;
const ll inf = 1e17 + 10;
const int N = 200000 + 10;
const ll mod = 1e9+7;
const int base=131;
mt19937 rng32(chrono::steady_clock::now().time_since_epoch().count());

map<string, int>ml;
map<int,int> mp;


ll b[N], vis[N], in[N],num[N],a[N],t, n, m, y, k,x;
int  cnt, sum, flag,f[N];

vector<int> v[N],ans;

struct mat{
	ll a[105][105];
	mat(){memset(a,0,sizeof(a));}
};

mat mul(mat a,mat b)
{
	mat ans;
	for(int i = 1; i <= m; i++){
        for(int j = 1; j <= m; j++){
            ans.a[i][j] = 0;
            for(int k = 1; k <= m; k++){
                ans.a[i][j] += a.a[i][k]*b.a[k][j];
			if(ans.a[i][j] > mod)ans.a[i][j] %= mod;
			}
		}
	}
	return ans;
}

mat ksm(mat a,ll b)
{
	mat ans;
	for(int i=1;i<=m;i++) ans.a[i][i]=1;
	while(b)
	{
		if(b&1) ans=mul(ans,a);
		b>>=1;
		a=mul(a,a);
	}
	return ans;
}
int main()
{

	ios::sync_with_stdio(false);
	cin.tie(0);



	cin>>n>>m;

	if(n<m) return cout<<1<<endl,0;
	mat tr;
	for(int i=1;i<m;i++)
		tr.a[i+1][i]=1;
	tr.a[1][1]=tr.a[1][m]=1;
	mat ans;
	for(int i=2;i<=m;i++) ans.a[i][1]=1;//行列关系不要搞错
	ans.a[1][1]=2;//1，1是dp[m]的值为2
	ans=mul(ksm(tr,n-m),ans);//这里是用tr乘ans，顺序不能搞反！！！
	cout<<ans.a[1][1]<<endl;

}

方法二：杜教BM

看到博客上有大神使用这个，第一次见，学到了！！！

这个板子功能超级强大，功能是输入前几项，能够自动推导公式后面的任一项。

前提：公式必须是线性的，前几项越多越好，一般不低于8项，不然推出来不准确。

这里就把前200项丢进去输出就行了。

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <vector>
#include <string>
#include <map>
#include <set>
#include <cassert>
#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
#define pb push_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define fi first
#define se second
#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef long long ll;
typedef pair<int,int> PII;
const ll mod=1000000007;
ll powmod(ll a,ll b) {ll res=1;a%=mod; assert(b>=0); for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}

int _,n;
namespace linear_seq {
    const int N=10010;
    ll res[N],base[N],_c[N],_md[N];

    vector<int> Md;
    void mul(ll *a,ll *b,int k) {
        rep(i,0,k+k) _c[i]=0;
        rep(i,0,k) if (a[i]) rep(j,0,k) _c[i+j]=(_c[i+j]+a[i]*b[j])%mod;
        for (int i=k+k-1;i>=k;i--) if (_c[i])
            rep(j,0,SZ(Md)) _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod;
        rep(i,0,k) a[i]=_c[i];
    }
    int solve(ll n,VI a,VI b) {
        ll ans=0,pnt=0;
        int k=SZ(a);
        assert(SZ(a)==SZ(b));
        rep(i,0,k) _md[k-1-i]=-a[i];_md[k]=1;
        Md.clear();
        rep(i,0,k) if (_md[i]!=0) Md.push_back(i);
        rep(i,0,k) 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;
                rep(j,0,SZ(Md)) res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod;
            }
        }
        rep(i,0,k) 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;
        rep(n,0,SZ(s)) {
            ll d=0;
            rep(i,0,L+1) 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.pb(0);
                rep(i,0,SZ(B)) 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.pb(0);
                rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod;
                ++m;
            }
        }
        return C;
    }
    int gao(VI a,ll n) {
        VI c=BM(a);
        c.erase(c.begin());
        rep(i,0,SZ(c)) c[i]=(mod-c[i])%mod;
        return solve(n,c,VI(a.begin(),a.begin()+SZ(c)));
    }
};
ll f[205];
int main() {
    ll n,m;
    cin>>n>>m;
    for(int i=1;i<m;i++) f[i]=1;
    f[m]=2;
    for(int i=m+1;i<=200;i++)
	f[i]=(f[i-1]+f[i-m])%mod;
    vector<int> v;
    for(int i=1;i<=200;i++)
	v.push_back(f[i]);
   cout<< linear_seq::gao(v,n-1);
}