[集训队互测 2012] calc【拉格朗日插值+ODE】
[集训队互测 2012] calc
先视序列为无序，最后答案乘上 n ! n! n! 即可。
那么问题变成了，在 [ 1 , A ] [1,A] [1,A] 中选 n n n 个数所有方案的权值之和。
这个东西的生成函数
∏
i
=
1
A
(
1
+
i
x
)
=
exp
(
∑
i
=
1
A
ln
(
1
+
i
x
)
)
\begin{aligned} &\prod_{i=1}^A (1+ix)\\ =&\exp(\sum_{i=1}^A\ln(1+ix)) \end{aligned}
=i=1∏A(1+ix)exp(i=1∑Aln(1+ix))
因为我们知道
−
ln
(
1
−
x
)
=
∑
i
>
1
x
i
/
i
-\ln(1-x)=\sum_{i>1}x^i/i
−ln(1−x)=∑i>1xi/i ，所以有
ln
(
1
+
i
x
i
)
=
−
∑
j
>
1
n
(
−
i
x
)
j
j
/
/
我
们
只
需
要
前
n
项
=
∑
j
>
1
n
(
−
1
)
j
+
1
(
i
x
)
j
j
∑
i
=
1
A
ln
(
1
+
i
x
i
)
=
∑
i
=
1
A
∑
j
>
1
n
(
−
1
)
j
+
1
(
i
x
)
j
j
=
∑
j
>
1
n
(
−
1
)
j
+
1
∑
i
=
1
A
i
j
j
x
j
\begin{aligned} \ln(1+ix^i)&=-\sum_{j>1}^n\frac{(-ix)^j}{j}//我们只需要前n项 \\ &=\sum_{j>1}^n\frac{(-1)^{j+1}(ix)^{j}}{j} \\ \sum_{i=1}^A \ln(1+ix^i)&=\sum_{i=1}^A\sum_{j>1}^n\frac{(-1)^{j+1}(ix)^{j}}{j} \\ &=\sum_{j>1}^n\frac{(-1)^{j+1}\sum_{i=1}^Ai^j}{j}x^{j} \end{aligned}
ln(1+ixi)i=1∑Aln(1+ixi)=−j>1∑nj(−ix)j//我们只需要前n项=j>1∑nj(−1)j+1(ix)j=i=1∑Aj>1∑nj(−1)j+1(ix)j=j>1∑nj(−1)j+1∑i=1Aijxj
然后就成了求经典的求自然数幂和的问题，通常有两种做法
求 ∑ i = 0 A i k \sum_{i=0}^A i^k ∑i=0Aik， A A A 在 int 范围， k ∈ [ 1 , n ] k\in [1,n] k∈[1,n] 。
- 拉格朗日插值 单个
k
k
k 时间复杂度
O
(
k
log
k
)
O(k\log k)
O(klogk) ，求所有
k
k
k 的时间复杂度
O
(
k
2
)
O(k^2)
O(k2) 。
f ( n ) = ∑ i = 0 k + 1 f ( i ) ∏ j = 0 , j ≠ i k + 1 n − j i − j = ∑ i = 0 k + 1 f ( i ) ∏ j = 0 i − 1 ( n − j ) ∑ j = i + 1 k + 1 ( n − j ) i ! ( − 1 ) k − i + 1 ( k + 1 − i ) ! = ∑ i = 0 k + 1 f ( i ) [ ∏ j = 0 k + 1 ( n − j ) ] / ( n − i ) i ! ( − 1 ) k − i + 1 ( k + 1 − i ) ! \begin{aligned} f(n)&=\sum_{i=0}^{k+1}f(i)\prod_{j=0,j\ne i}^{k+1}\frac{n-j}{i-j}\\ &=\sum_{i=0}^{k+1}f(i)\frac{\prod_{j=0}^{i-1}(n-j)\sum_{j=i+1}^{k+1}(n-j)}{i!(-1)^{k-i+1}(k+1-i)!}\\ &=\sum_{i=0}^{k+1}f(i)\frac{[\prod_{j=0}^{k+1}(n-j)]/(n-i)}{i!(-1)^{k-i+1}(k+1-i)!} \end{aligned} f(n)=i=0∑k+1f(i)j=0,j=i∏k+1i−jn−j=i=0∑k+1f(i)i!(−1)k−i+1(k+1−i)!∏j=0i−1(n−j)∑j=i+1k+1(n−j)=i=0∑k+1f(i)i!(−1)k−i+1(k+1−i)![∏j=0k+1(n−j)]/(n−i) - 这个东西的EGF，多项式求逆即可，能求所有的
k
∈
[
1
,
n
]
k\in[1,n]
k∈[1,n] ，
O
(
n
log
n
)
O(n\log n)
O(nlogn) 。
∑ k = 0 ∞ ∑ i = 0 A i k x k k ! = ∑ i = 0 A ∑ k = 0 ∞ ( i x ) k k ! = ∑ i = 0 A e i x = e x ( A + 1 ) − 1 e x − 1 \begin{aligned} &\sum_{k=0}^{\infty}\sum_{i=0}^{A}i^k\frac{x^k}{k!}\\ =&\sum_{i=0}^{A}\sum_{k=0}^{\infty}\frac{{(ix)}^k}{k!}\\ =&\sum_{i=0}^Ae^{ix}\\ =&\frac{e^{x(A+1)-1}}{e^x-1} \end{aligned} ===k=0∑∞i=0∑Aikk!xki=0∑Ak=0∑∞k!(ix)ki=0∑Aeixex−1ex(A+1)−1
下面的多项式常数项为 0 0 0 不能直接求逆，所以要先上下同时除以 x x x 。
这道题的模数不指定，所以拉格朗日插值 O ( n 2 ) O(n^2) O(n2) 算。（苣佬可以用MTT，orz
隔壁加强版就套上多项式板子 O ( n log n ) O(n\log n) O(nlogn) 解决。
最后exp O ( n 2 ) O(n^2) O(n2) 暴力算。具体就是 f = e g , f ′ = f ⋅ g ′ f=e^g,f'=f\cdot g' f=eg,f′=f⋅g′ 。
#include <bits/stdc++.h>
#define N 503
using namespace std;
typedef long long ll;
ll mod,n,A,m;
ll ans[N],f[N],now[N],anss[N];
ll fac[N],inv[N],invv[N],kk[N],bb[N];
ll ksm(ll x,ll y){
ll res=1;
while(y){ if(y&1) res=res*x%mod; x=x*x%mod; y>>=1; }
return res;
}
ll solve(int k){
for(ll i=1;i<=n+1;i++) f[i]=((kk[i]=(kk[i]*i)%mod)+f[i-1])%mod;
if(A<=k+1) return f[A];
ll aa=1,res=0,p;
bb[k+2]=1;
for(ll i=k+1;i>=1;i--) bb[i]=(A-i)*bb[i+1]%mod;
for(ll i=1;i<=k+1;i++){
aa=aa*(A-i+1)%mod;
p=1; if((k+1-i)&1) p=mod-1;
(res+=f[i]*aa%mod*bb[i+1]%mod*inv[i]%mod*p%mod*inv[k+1-i]%mod)%=mod;
}
return res;
}
void init(){
fac[0]=1;
for(ll i=1;i<N;i++) fac[i]=fac[i-1]*i%mod;
inv[N-1]=ksm(fac[N-1],mod-2);
for(ll i=N-1;i>=1;i--) inv[i-1]=inv[i]*i%mod,invv[i]=inv[i]*fac[i-1]%mod;
}
int main(){
cin>>A>>n>>mod;
init();
for(int i=0;i<=n+1;i++) kk[i]=1;
ll p;
for(int i=1;i<=n;i++){
ans[i]=solve(i);
p=1; if((i+1)&1) p=mod-1;
ans[i]=p*ans[i]%mod*invv[i]%mod;
}
for(ll i=1;i<=n;i++)
ans[i-1]=ans[i]*i%mod;
now[0]=1;
for(int i=1;i<=n;i++){
for(int j=0;j<i;j++)
now[i]=(now[i]+now[j]*ans[i-j-1]%mod)%mod;
now[i]=now[i]*invv[i]%mod;
}
printf("%lld",now[n]*fac[n]%mod);
}
