生成函数 常系数齐次线性递推 P2461题解
Prean
·
·
题解
引用化学老师的一句话：什么矩阵，没有矩阵！
这种板子题怎么能用矩阵呢。
首先设 $ F_n(x)=x^n \bmod {1-P(x)} $，那么我们需要求 $ \sum_{i=1}^n F_i(x) \bmod (1-P(x)) $。然后卷上 $ B(x) $ 就可以得到需要的东西了。
注意到这是等比数列求和，可以使用分治计算等比数列，可以保证复杂度是 $ O(k^2\log n) $ 而不是 $ O(k^2\log^2n) $ 的。
主要到我们在求的实际上是 $ \sum_{i=1}^n (F_i(x) \bmod {P(x)}) $，但是因为这些东西加起来再取模和取模之后再加起来的结果是一样的，所以并无区别。
坑还是比较多的，需要注意一下。
```cpp
#include<cstdio>
typedef unsigned ui;
typedef __uint128_t L;
typedef unsigned long long ull;
const ui M=55;
ui len,P,b[M],p[M];ull n,m;
struct Barrett{
ull b,m;
Barrett(const ull&m=1):m(m),b((L(1)<<64)/m){}
friend inline ull operator%(const ull&a,const Barrett&mod){
ull r=a-mod.m*(L(mod.b)*a>>64);return r>=mod.m?r-mod.m:r;
}
}mod;
inline void add(ui*f,ui*g,const ui&len){
ui i;for(i=0;i^len;++i)f[i]=(f[i]+g[i])%mod;
}
inline void times(ui*f,ui*g,ui*P,const ui&len){
ui i,j,t,x;static ui sav[M];
for(i=0;i^len;++i)if(f[i])for(j=0;j^len;++j)if(g[j])sav[i+j]=(sav[i+j]+1ull*f[i]*g[j])%mod;
for(i=(len<<1)-1;i>=len;--i)if(sav[i])for(t=sav[i],j=len;j<=len;--j)sav[i-j]=(sav[i-j]+1ull*t*P[j])%mod;
for(i=0;i^len;++i)f[i]=sav[i],sav[i]=0;
}
inline ui Solve(ui*b,ui*P,const ui&len,ull n){
if(n>>63)return 0;ui i,ans(0);static ui f[M],g[M],sav[M];sav[0]=g[0]=1;if(len^1)f[1]=1;else f[0]=p[1];
for(;n;n>>=1,++f[0],times(g,f,P,len),--f[0],times(f,f,P,len))if(n&1)times(sav,f,P,len),add(sav,g,len);
for(i=0;i^len;++i)ans=(ans+1ull*sav[i]*b[i+1])%mod,f[i]=g[i]=sav[i]=0;return ans;
}
signed main(){
ui i;scanf("%u",&len);for(i=1;i<=len;++i)scanf("%u",b+i);for(i=1;i<=len;++i)scanf("%u",p+i);
scanf("%llu%llu%u",&n,&m,&P);mod=Barrett(P);p[0]=P-1;for(i=1;i<=len;++i)b[i]=b[i]%mod,p[i]=p[i]%mod;
printf("%u",(Solve(b,p,len,m-1)+P-Solve(b,p,len,n-2))%mod);
}
```
### upd:这道题可以使用新算法。
我们观察得到，答案为 $ [x^n]\frac {B(x)(1-C(x))} {(1-C(x))(1-x)} $。
然后跑一遍老算法，但是需要求逆。
然而注意到分母的零次项一定为 $ 1 $，所以实际上并不需要求逆。
常数比老算法小一点儿，仍然不清楚最优解是什么。。。
```cpp
#include<cstdio>
typedef unsigned ui;
typedef __uint128_t L;
typedef unsigned long long ull;
const ui M=55;
ui len,P,f[M],g[M],b[M],p[M];ull n,m;
struct Barrett{
ull b,m;
Barrett(const ui&m=1):m(m),b((L(1)<<64)/m){}
friend inline ull operator%(const ull&a,const Barrett&mod){
ull r=a-mod.m*(L(mod.b)*a>>64);return r>=mod.m?r-mod.m:r;
}
}mod;
inline void times(ui*f,ui*g,const ui&len){
ui i,j,t;static ui sav[M];
for(i=0;i^len;++i)if(f[i])for(j=0;j^len;++j)if(g[j])sav[i+j]=(sav[i+j]+1ull*f[i]*g[j])%mod;
for(i=0;i<len*2;++i)f[i]=sav[i],sav[i]=0;
}
inline ui Solve(ui*f,ui*g,const ui&len,ull n){
ui i;static ui sav[M];
for(;n;n>>=1){
for(i=0;i<len;++i)sav[i]=i&1?P-g[i]:g[i];times(f,sav,len);times(g,sav,len);
for(i=n&1;i<len*2;i+=2)f[i>>1]=f[i];for(i=0;i<len*2;i+=2)g[i>>1]=g[i];for(i=len;i<len*2;++i)f[i]=g[i]=0;
}
return f[0];
}
signed main(){
ui i,x,y;scanf("%u",&len);++len;for(i=1;i^len;++i)scanf("%u",b+i);for(i=1;i^len;++i)scanf("%u",p+i);
scanf("%llu%llu%u",&n,&m,&P);mod=Barrett(P);p[0]=1;for(i=1;i^len;++i)b[i]=b[i]%mod,p[i]=P-p[i]%mod;
times(b,p,len);b[len++]=0;for(i=len-1;i;--i)p[i]=(p[i]+P-p[i-1])%mod;
for(i=0;i^len;++i)f[i]=b[i],g[i]=p[i];x=Solve(f,g,len,n-1);
for(i=0;i^len;++i)f[i]=b[i],g[i]=p[i];y=Solve(f,g,len,m);
printf("%u",(P+y-x)%mod);
}
```