标签：har 矩阵加速 现在 一个 [1] operator a* bsp 方式
题目看起来非常复杂。
本质不同的细胞这个条件显然太啰嗦，
是否有些可以挖掘的性质？
1.发现，只要第一次分裂不同，那么互相之间一定是不同的（即使总数目相同）。
所以先考虑第一次分裂后，一个固定小球体数量的情况：
2.第一次分裂后，最后的小球体数量固定。想要方案数不同，必须连接方式不同。
可以列出dp式子，f[n]（以n结尾砍一刀）=f[n-2]+f[n-3]+...+f[2]+f[0]，而f[0]=1,f[1]=0
而fibo[n]-1=f[n-2]+...+f[1]
可以猜想和斐波那契数列有关。
实际上，f[n]=fibo[n-1],f[1]=0,f[0]=1
所以，如果我们知道第一次分裂的方案，
即S=a1+...+am（分成m段，ai表示这一段拼成的数字）
那么，f[S]=fibo[S-1]
好了，
所以，题意其实是：
? 给出一个数字序列
? 你要先把序列切成很多段，把每一段看做一个十进制数
??1, ??2, … , ????
? 你要求每一种分割方案的 Fib （??1 + ? + ???? ）的和
? ?? ≤ 1000
现在我们的问题就是怎样快速统计这些所有fibo的和。
设f[i]表示，前i个，i后劈一刀，所有产生的Fib （??1 + ? + ??k ）的和。
由于矩阵有结合律和分配律，
Fibo[i]=Start*(M^j*M^(i-j))（M是转移矩阵）
我可以把f用矩阵表示
经典n^2的dp
然后，我们要快速统计M^(num[j...i])
直接每次高精精快速幂？？（逗笑）
根据套路，可以利用j每次向左走一位，在路上利用之前的基础统计出num[j...i]
那么，我么就要知道所有的K=a*10^p，M^K可以预处理。
然后直接M^num=M^(K+old)=M^k*M^old
即可统计出num[j...i]
转移的时候，f[i]+=f[j]*M^(num[j+1...i])即可，因为乘法分配律，dp式子成立。
注意下，f[S]=Fibo[S-1]
所以，其实f[0]=[0,1,0,0]（写成矩阵）
总之，如果j==0，那么这个段是第一段，直接把这个1减掉。相当于对矩阵乘一个逆
代码：
（注意，逆矩阵：[mod-1,1,1,0]而不是[-1,1,1,0]。不要爆负数）
#include<bits/stdc++.h> using namespace std; typedef long long ll; const int N=1005; const int mod=1e9+7; int n; char s[N]; struct tr{ ll a[3][3]; void pre(){ memset(a,0,sizeof a); } void init(){ a[1][1]=a[2][2]=1; } tr operator *(const tr &b){ tr c;c.pre(); for(int i=1;i<=2;i++) for(int k=1;k<=2;k++) for(int j=1;j<=2;j++) c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j])%mod; return c; } tr operator +(const tr &b){ tr c;c.pre(); for(int i=1;i<=2;i++) for(int j=1;j<=2;j++) c.a[i][j]=(a[i][j]+b.a[i][j])%mod; return c; } void op(){ cout<<endl; for(int i=1;i<=2;i++){ for(int j=1;j<=2;j++){ cout<<a[i][j]<<" "; }cout<<endl; } cout<<endl; } }B[10][1005],f[1005],A,C,D; tr qm(tr c,ll y){ tr ret;ret.pre();ret.init(); while(y){ if(y&1) ret=ret*c; c=c*c; y>>=1; } return ret; } void prewrk(){ for(int i=0;i<=9;i++) B[i][0]=qm(A,i); for(int i=1;i<=n;i++){ for(int j=0;j<=9;j++){ B[j][i]=qm(B[j][i-1],10); } } } int main(){ scanf("%d",&n); A.a[1][1]=0;A.a[1][2]=1;A.a[2][2]=1;A.a[2][1]=1; D.a[1][1]=mod-1;D.a[1][2]=1;D.a[2][1]=1;D.a[2][2]=0; prewrk(); scanf("%s",s+1); //cout<<" hh "<<endl; f[0].a[1][2]=1; for(int i=1;i<=n;i++){ C.pre();C.init(); //cout<<i<<"------------"<<endl; for(int j=i-1;j>=0;j--){ //cout<<s[j+1]-‘0‘<<endl; C=C*B[s[j+1]-‘0‘][i-j-1]; //C.op(); //cout<<">>"<<endl; if(j==0) C=C*D; tr tmp=f[j]*C; f[i]=f[i]+tmp; //cout<<" f[i] "<<endl; //f[i].op(); } } printf("%lld\n",f[n].a[1][1]); return 0; }
标签：har 矩阵加速 现在 一个 [1] operator a* bsp 方式
原文地址：https://www.cnblogs.com/Miracevin/p/9748231.html
