参考资料:
拉格朗日插值法(图文详解),
a
t
t
a
c
k
attack
attack,
c
o
m
m
a
n
d
_
b
l
o
c
k
command\_block
command_block
一些定义
拉格朗日插值法 是已知点值 求多项式系数的方法:
已知
n
+
1
n+1
n+1 个点值
(
x
i
,
y
i
)
(
i
∈
[
0
,
n
]
∩
N
)
(x_i,y_i)(i\in [0,n]\cap \N)
(xi,yi)(i∈[0,n]∩N), 求一个
n
n
n 次多项式
L
(
x
)
L(x)
L(x),使得
L
(
x
i
)
=
y
i
L(x_i)=y_i
L(xi)=yi.
求法:
L
(
x
)
=
∑
i
=
0
n
ℓ
j
(
x
)
y
j
,
ℓ
j
(
x
)
=
∏
j
≠
i
x
−
x
i
x
j
−
x
i
L(x)=\sum_{i=0}^n \ell_j(x)y_j,~~~~~\ell_j(x)=\prod_{j\ne i} \dfrac{x-x_i}{x_j-x_i}
L(x)=∑i=0nℓj(x)yj, ℓj(x)=∏j=ixj−xix−xi.
容易发现
ℓ
j
(
x
j
)
=
1
,
ℓ
j
(
x
i
)
=
0
(
i
≠
j
)
\ell_j(x_j)=1,\ell_j(x_i)=0(i\ne j)
ℓj(xj)=1,ℓj(xi)=0(i=j) .
所以
L
(
x
i
)
=
(
∑
j
≠
i
ℓ
j
(
x
i
)
)
+
ℓ
i
(
x
i
)
=
0
+
y
i
=
y
i
L(x_i)=\left( \sum_{j\ne i}\ell_j(x_i) \right) + \ell_i(x_i)=0+y_i=y_i
L(xi)=(∑j=iℓj(xi))+ℓi(xi)=0+yi=yi.
我们称 L ( x ) L(x) L(x) 为拉格朗日插值多项式, ℓ j ( x ) \ell_j(x) ℓj(x) 为拉格朗日基本多项式或插值基函数.
性质:
- 拉格朗日插值多项式的唯一性.假设存在两个 n n n次多项式 A , B A,B A,B,那么 A − B A-B A−B在 x i x_i xi的取值都为 0 0 0.所以 A − B = λ ( ∏ i x − x i ) A-B=\lambda (\prod_i x-x_i) A−B=λ(∏ix−xi),这显然是一个 n + 1 n+1 n+1次多项式,所以和 A , B A,B A,B均为 n n n次多项式矛盾, A − B = 0 A-B=0 A−B=0.
- 插值基函数 组成 一个 n n n次多项式的线性空间. 假设存在系数 λ 0 , λ 1 , λ 2 . . . , λ n \lambda_0,\lambda_1,\lambda_2...,\lambda_n λ0,λ1,λ2...,λn使得 多项式 P = ∑ i = 0 n λ i ℓ i = 0 P=\sum_{i=0}^n \lambda_i \ell_i =0 P=∑i=0nλiℓi=0,则 P P P为 ( x i , λ i ) (x_i,\lambda_i) (xi,λi)的拉格朗日插值多项式,因为 P = 0 P=0 P=0,所以 λ i = 0 \lambda_i=0 λi=0,所以 ℓ i \ell_i ℓi线性无关.
模板
需要注意的是我们并不需要求出拉格朗日插值多项式的具体每一项的系数,只需要记录其点值
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi)即可.
void solve(int n,ll k,int *x,int *y) {
ll ans=0;
for(int i=0;i<n;i++) {
ll s1=1,s2=1;
for(int j=0;j<n;j++) if(i^j) {
s1=s1*(k-x[j]+mod)%mod;
s2=s2*(x[i]-x[j]+mod)%mod;
}
ans+=s1*inv(s2)%mod*y[i]%mod;
}
pr2(ans%mod);
}
重心拉格朗日插值法
可以发现
ℓ
i
(
x
)
\ell_i(x)
ℓi(x) 的分子部分有大部分重复部分,所以设
l
(
x
)
=
∏
x
−
x
i
l(x)=\prod x-x_i
l(x)=∏x−xi.
同时设
w
i
=
∏
j
≠
i
1
x
i
−
x
j
w_i=\prod_{j\ne i} \dfrac 1{x_i-x_j}
wi=∏j=ixi−xj1.
则有 :
L
(
x
)
=
l
(
x
)
∑
i
=
0
n
w
i
⋅
y
i
x
−
x
i
(1)
L(x)=l(x)\sum_{i=0}^n \dfrac {w_i\cdot y_i}{x-x_i}\tag 1
L(x)=l(x)i=0∑nx−xiwi⋅yi(1).
这样对于带加入操作的情况,我们就
O
(
n
)
O(n)
O(n)修改
w
w
w即可.
单次代入也是
O
(
n
)
O(n)
O(n)的了.,
我们设
g
(
x
)
=
1
g(x)=1
g(x)=1,用插值法得到:
g
(
x
)
=
l
(
x
)
∑
i
=
0
n
w
i
x
−
x
i
g(x)=l(x)\sum_{i=0}^n \dfrac{w_i}{x-x_i}
g(x)=l(x)i=0∑nx−xiwi
L
(
x
)
=
L
(
x
)
g
(
x
)
=
∑
i
=
0
n
w
i
⋅
y
i
x
−
x
i
∑
i
=
0
n
w
i
x
−
x
i
(2)
L(x)=\dfrac {L(x)}{g(x)}=\dfrac{\sum_{i=0}^n \dfrac{w_i\cdot y_i}{x-x_i} } {\sum_{i=0}^n \dfrac {w_i}{x-x_i} } \tag 2
L(x)=g(x)L(x)=∑i=0nx−xiwi∑i=0nx−xiwi⋅yi(2)
实际效果同上.
x x x 连续的问题
此时
x
i
=
i
,
i
∈
[
0
,
n
]
x_i=i,i\in[0,n]
xi=i,i∈[0,n].
则有
L
(
x
)
=
∑
i
=
0
n
y
i
∏
j
≠
i
x
−
x
j
x
i
−
x
j
=
∑
i
=
0
n
y
i
∏
j
≠
i
x
−
j
i
−
j
L(x)=\sum_{i=0}^n y_i \prod_{j\ne i}\dfrac {x-x_j}{x_i-x_j}=\sum_{i=0}^n y_i \prod_{j\ne i}\dfrac {x-j}{i-j}
L(x)=∑i=0nyi∏j=ixi−xjx−xj=∑i=0nyi∏j=ii−jx−j
这种题目我们可以
O
(
n
)
O(n)
O(n)处理.
定义
p
r
e
[
i
]
=
∏
j
=
0
i
(
x
−
j
)
,
s
u
f
[
i
]
=
∏
j
=
i
n
(
x
−
j
)
,
i
n
v
[
i
]
=
∏
j
=
1
i
1
j
pre[i]=\prod_{j=0}^i (x-j),suf[i]=\prod_{j=i}^n (x-j),inv[i]=\prod_{j=1}^i \dfrac 1j
pre[i]=∏j=0i(x−j),suf[i]=∏j=in(x−j),inv[i]=∏j=1ij1.
则有
L
(
x
)
=
∑
y
i
⋅
p
r
e
[
i
−
1
]
⋅
s
u
f
[
i
+
1
]
i
n
v
[
i
]
∗
i
n
v
[
n
−
i
]
∗
(
−
1
)
n
−
i
L(x)=\sum y_i \cdot \dfrac{pre[i-1]\cdot suf[i+1]}{inv[i]*inv[n-i]*(-1)^{n-i}}
L(x)=∑yi⋅inv[i]∗inv[n−i]∗(−1)n−ipre[i−1]⋅suf[i+1].
自然数幂之和及其差分思想
S
(
x
,
k
)
=
f
(
x
)
=
∑
i
=
1
x
i
k
S(x,k)=f(x) = \sum_{i=1}^x i ^k
S(x,k)=f(x)=∑i=1xik.
这是一个
k
+
1
k+1
k+1 次的多项式, 可以运用上面一个部分的方法实现
O
(
k
log
k
)
O(k\log k)
O(klogk) 求解.
差分法证明:
- 设
n
n
n 次多项式
g
(
x
)
=
∑
i
=
0
n
a
i
x
i
g(x) = \sum_{i=0}^n a_i x^i
g(x)=∑i=0naixi.
对于常数 D ≠ 0 D\ne 0 D=0,设 g ( x ) − g ( x − D ) = ∑ i = 0 n b i x i g(x)-g(x-D)= \sum_{i=0}^n b_i x^i g(x)−g(x−D)=∑i=0nbixi.
我们利用二项式定理,对 ( x − D ) n (x-D)^n (x−D)n进行展开,可以发现 b n = 0 b_n=0 bn=0.
所以, g ( x ) g(x) g(x)的一阶差分是 n − 1 n-1 n−1 次多项式. - 运用数学归纳法推理:
当 k = 0 k=0 k=0时, S ( n , k ) = n S(n,k)=n S(n,k)=n,成立.
否则:
( n + 1 ) k + 1 − 1 = ∑ i = 1 n ( i + 1 ) k + 1 − i k + 1 = ∑ i = 1 n ∑ j = 0 k ( k + 1 j ) i j = ∑ j = 0 k ( k + 1 j ) ∑ i = 1 n i j = ∑ j = 0 k ( k + 1 j ) S ( n , j ) (n+1)^{k+1}-1=\sum_{i=1}^n (i+1)^{k+1}-i^{k+1}=\sum_{i=1}^n\sum_{j=0}^k \dbinom {k+1}{j}i^j=\sum_{j=0}^k\dbinom {k+1}j \sum_{i=1}^n i^j=\sum_{j=0}^k\dbinom {k+1}j S(n,j) (n+1)k+1−1=∑i=1n(i+1)k+1−ik+1=∑i=1n∑j=0k(jk+1)ij=∑j=0k(jk+1)∑i=1nij=∑j=0k(jk+1)S(n,j).
故 S ( n , k ) = ( n + 1 ) k + 1 − 1 − ∑ j = 0 k − 1 ( k + 1 j ) S ( n , j ) k + 1 S(n,k) = \dfrac{(n+1)^{k+1} -1 -\sum_{j=0}^{k-1}\dbinom {k+1}j S(n,j)}{k+1} S(n,k)=k+1(n+1)k+1−1−∑j=0k−1(jk+1)S(n,j).
因为 S ( n , j ) S(n,j) S(n,j)为 j + 1 j+1 j+1次多项式,所以次数为 max ( k + 1 , k ) = k + 1 \max(k+1,k)=k+1 max(k+1,k)=k+1.
推论: 任意一阶差分为 n n n次多项式,则函数为 n + 1 n+1 n+1次多项式.
例题
P4593 [TJOI2018]教科书般的亵渎
显然法术使用次数为
m
+
1
m+1
m+1,所以我们需要
m
+
2
m+2
m+2个点值.
可以暴力拉格朗日插值,复杂度为
O
(
T
m
3
)
O(Tm^3)
O(Tm3).
ll f(ll n,ll m,ll *x,ll *y) {
ll res=0; n=n%mod+mod;
for(int i=0;i<=m;i++) {
ll a=1,b=1;
for(int j=0;j<=m;j++) if(i^j)
a=a*(n-x[j])%mod,b=b*(x[i]-x[j]+mod)%mod;
res += y[i]*a%mod*power(b)%mod;
}
return res%mod;
}
int main() {
qr(T); while(T--) {
qr(n); qr(m);
for(int i=1;i<=m;i++) qr(a[i]);
sort(a+1,a+m+1); a[++m]=++n;
y[0]=0; for(int i=1;i<=m+1;i++) x[i]=i,y[i]=(y[i-1]+power(i,m))%mod;
ll ans=0;
for(int i=0;i<=m;i++) {
ans += f(n-a[i],m+1,x,y);
for(int j=i+1;j<=m;j++) ans -= power((a[j]-a[i])%mod,m);
}
pr2((ans%mod+mod)%mod);
}
}
也可以应用上面所述的 x x x连续的方法计算,复杂度为 O ( T m 2 ) O(Tm^2) O(Tm2).
int T;
ll n,m,y[N],pre[N],suf[N],inv[N],a[N];
ll power(ll a,ll b=mod-2) {
ll c=1;
while(b) {
if(b&1) c=c*a%mod;
b /= 2; a=a*a%mod;
} return c;
}
ll f(ll n,int m) {
int t=m+1; n=n%mod+mod; pre[0]=n; suf[t+1]=1; y[0]=0;
for(int i=1;i<=t;i++) pre[i]=pre[i-1]*(n-i)%mod;
for(int i=t;i>=1;i--) suf[i]=suf[i+1]*(n-i)%mod;
for(int i=1;i<=t;i++) y[i]=(y[i-1]+power(i,m))%mod;
ll res=0;
for(int i=1;i<=t;i++)
res += y[i]*pre[i-1]%mod*suf[i+1]%mod*inv[i]%mod*((t-i)&1?mod-inv[t-i]:inv[t-i])%mod;
return res%mod;
}
int main() {
m=53; inv[0]=inv[1]=1; for(int i=2;i<=m;i++) inv[i]=inv[mod%i]*(mod-mod/i)%mod;
for(int i=2;i<=m;i++) inv[i]=inv[i-1]*inv[i]%mod;
qr(T); while(T--) {
qr(n); qr(m);
for(int i=1;i<=m;i++) qr(a[i]);
sort(a+1,a+m+1); a[++m]=++n; ll ans=0;
for(int i=0;i<=m;i++) {
ans += f(n-a[i],m);
for(int j=i+1;j<=m;j++) ans -= power((a[j]-a[i])%mod,m);
}
pr2((ans%mod+mod)%mod);
} return 0;
}
bzoj #3453. tyvj 1858 XLkxc
给 定 k , a , n , d , p = 1234567891. 给定 k,a,n,d,p=1234567891. 给定k,a,n,d,p=1234567891.
f ( i ) = 1 k + 2 k + 3 k + . . . . . . + i k f(i)=1^k+2^k+3^k+......+i^k f(i)=1k+2k+3k+......+ik
g ( x ) = f ( 1 ) + f ( 2 ) + f ( 3 ) + . . . . + f ( x ) g(x)=f(1)+f(2)+f(3)+....+f(x) g(x)=f(1)+f(2)+f(3)+....+f(x)
求 ( g ( a ) + g ( a + d ) + g ( a + 2 d ) + . . . . . . + g ( a + n d ) ) m o d p 求(g(a)+g(a+d)+g(a+2d)+......+g(a+nd))\mod p 求(g(a)+g(a+d)+g(a+2d)+......+g(a+nd))modp
1 ≤ k ≤ 1230 ≤ a , n , d ≤ 123456789 1\le k\le 123 0\le a,n,d\le 123456789 1≤k≤1230≤a,n,d≤123456789
f f f 为 k + 1 k+1 k+1 次, g g g为 k + 2 k+2 k+2次, h ( x ) = ∑ i = 0 x g ( a + i ∗ d ) h(x)=\sum_{i=0}^x g(a+i*d) h(x)=∑i=0xg(a+i∗d)为 k + 3 k+3 k+3次.
注意 2 p 2p 2p会爆 i n t int int.
ll power(ll a, ll b = mod - 2) {
ll c = 1;
while(b) {
if (b & 1) c = c * a % mod;
b /= 2; a = a * a % mod;
}
return c;
}
struct LA {
int n;
ll y[N];
ll& operator [](int i) {return y[i];}
int F(ll k) {
k %= mod;
ll ans=0;
rep(i, 0, n) {
ll a = y[i], b = 1;
rep(j, 0, n) if(j ^ i)
a = a * (k - j) % mod,
b = b * (i - j) % mod;
ans = (ans + a * power(b)) % mod;
}
return (ans + mod) % mod;
}
} f, g;
void solve() {
ll k, a, n, d;
qr(k); qr(a); qr(n); qr(d);
f.n = k + 3;
FOR(i, f.n) f[i] = (f[i - 1] + power(i, k)) % mod;
FOR(i, f.n) f[i] = (f[i - 1] + f[i]) % mod;
g.n = f.n;
rep(i, 0, g.n) g[i] = ((i ? g[i - 1] : 0) + f.F(a + i * d)) % mod;
pr2(g.F(n));
}
P4463 [集训队互测2012] calc
选出一个大小为 n n n的序列 a a a,使得每个数都在 [ 1 , k ] [1,k] [1,k]内,且两两不相等.
定义一个方案的权值为 ∏ i a i \prod_i a_i ∏iai. 求所有方案的总和 m o d p \mod p modp.
n ≤ 500 , k ≤ 1 0 9 , p ∈ P r i m e n\le 500, k \le 10^9,p\in Prime n≤500,k≤109,p∈Prime.
首先把序列转化为集合.
然后定义状态
f
[
n
]
[
k
]
f[n][k]
f[n][k],有:
f
[
n
]
[
k
]
=
f
[
n
−
1
]
[
k
−
1
]
∗
k
+
f
[
n
]
[
k
−
1
]
f[n][k]=f[n-1][k-1]*k + f[n][k-1]
f[n][k]=f[n−1][k−1]∗k+f[n][k−1].
然后,由于
k
k
k很大,我们猜想
f
[
n
]
f[n]
f[n]是一个
g
(
n
)
g(n)
g(n)次多项式.
f
[
n
]
[
k
]
−
f
[
n
]
[
k
−
1
]
=
f
[
n
−
1
]
[
k
−
1
]
∗
k
f[n][k]-f[n][k-1]= f[n-1][k-1]*k
f[n][k]−f[n][k−1]=f[n−1][k−1]∗k.
右边是
g
(
n
−
1
)
+
1
g(n-1)+1
g(n−1)+1次,故
f
[
n
]
f[n]
f[n]是
g
(
n
)
=
g
(
n
−
1
)
+
1
+
1
=
g
(
n
−
1
)
+
2
g(n)=g(n-1)+1+1=g(n-1)+2
g(n)=g(n−1)+1+1=g(n−1)+2.
显然
g
(
0
)
=
0
g(0)=0
g(0)=0,所以
g
(
n
)
=
2
n
g(n)=2n
g(n)=2n.
我们只用跑
n
∗
3
n
n*3n
n∗3n的
d
p
dp
dp,就可以插值啦.
这道题是用 插值快速维护
d
p
dp
dp 的方法.
总复杂度为 O ( n 2 ) O(n^2) O(n2).
struct LA {
int n, x[N], y[N];
int& operator [](int i) {return y[i];}
int operator () (int k) {
ll ans=0;
rep(i, 0, n) {
ll a = y[i], b = 1;
rep(j, 0, n) if(j ^ i)
a = a * (k - x[j]) % mod,
b = b * (x[i] - x[j]) % mod;
ans = (ans + a * power(b)) % mod;
}
return ans;
}
} f;
ll v[N];
void solve() {
qr(k); qr(n); qr(mod);
f.n = 2 * n;
fill(v, v + n + f.n, 1);
ll s = 1;
FOR(i, n) {
s = s * i % mod;
ll t = v[i - 1], tmp; v[i - 1] = 0;
rep(j, i, n + f.n) {
tmp = v[j];
v[j] = (t * j + v[j - 1]) % mod;
t = tmp;
}
}
for(int i = 0, j = n; i <= f.n; i++, j++)
f.x[i] = j, f[i] = v[j] * s % mod;
pr2(f(k));
}
bzoj #4559. [JLoi2016]成绩比较
容斥+ d p dp dp.
对于
S
(
n
,
i
)
,
i
∈
[
0
,
k
)
S(n,i),i\in [0,k)
S(n,i),i∈[0,k)
我学会了复杂的优化复杂度的做法: 即用二项式定理划开式子进行
d
p
dp
dp.
可以由
O
(
k
2
log
k
)
O(k^2\log k)
O(k2logk)变成
O
(
k
2
)
O(k^2)
O(k2).
但是实际上,我们如果幂用递推求,然后用拉格朗日插值也是
O
(
k
2
)
O(k^2)
O(k2)的.
void init() {
jc[0] = inv[0] = 1;
FOR(i, n + 5) jc[i] = jc[i - 1] * i % mod, inv[i] = power(jc[i]);
}
ll C(int x,int y) {return jc[x] * inv[y] % mod * inv[x - y] % mod;}
ll g[N];
void solve() {
qr(n); qr(m); qr(k); init();
ll ans = 0;
int mn = n - 1;
FOR(i, m) qr(U[i]);
FOR(i, m) qr(R[i]), cmin(mn, n - R[i]);
rep(i, k, mn) {
ll s = C(n - 1, i);
FOR(j, m) s = s * C(n - i - 1, R[j] - 1) % mod;
ans = (ans + s * ((i - k) & 1? mod - 1 : 1) % mod * C(i, k)) % mod;
}
FOR(i, m) {
int u = U[i] + 1;
ll p = u;
FOR(j, n) {
g[j - 1] = p - 1;
p = p * u % mod;
rep(t, 0, j - 2)
g[j - 1] = (g[j - 1] - C(j, t) * g[t]) % mod;
g[j - 1] = (g[j - 1] + mod) * power(j) % mod;
}
ll sum = 0; p = 1;
rep(j, 0, R[i] - 1) {
sum = (sum + p * ((R[i] - j - 1) & 1 ? mod - 1 : 1) % mod * C(R[i] - 1, j) % mod * g[n - j - 1]) % mod;
p = p * U[i] % mod;
}
ans = ans * sum % mod;
}
pr2(ans);
}
P6271 [湖北省队互测2014]一个人的数论
给定 n n n的因式分解: n = ∏ i = 1 w p i a i n=\prod_{i=1}^w p_i^{a_i} n=∏i=1wpiai.
求 ∑ i = 1 n i k [ g c d ( i , n ) = 1 ] \sum_{i=1}^n i^k[gcd(i,n) = 1] ∑i=1nik[gcd(i,n)=1].
k ≤ 100 , w ≤ 1000 , a i , p i ≤ 1 0 9 k\le 100,w\le 1000,a_i,p_i\le 10^9 k≤100,w≤1000,ai,pi≤109.
构造
F
(
n
)
=
∑
i
=
1
n
i
k
,
G
(
n
)
=
∑
i
=
1
n
i
k
[
g
c
d
(
i
,
n
)
=
1
]
F(n)=\sum _{i=1}^n i^k,G(n) =\sum_{i=1}^ni^k[gcd(i,n)=1]
F(n)=∑i=1nik,G(n)=∑i=1nik[gcd(i,n)=1].
有
F
(
n
)
=
∑
d
∣
n
d
k
G
(
n
d
)
F(n)=\sum_{d|n} d^k G (\dfrac n d)
F(n)=∑d∣ndkG(dn).
运用反演:
G
(
n
)
=
∑
d
∣
n
d
k
μ
(
d
)
F
(
n
d
)
G(n) = \sum_{d|n} d^k\mu(d) F(\dfrac n d)
G(n)=∑d∣ndkμ(d)F(dn).
由于
n
n
n较大,我们必须利用积性函数的性质来优化.
而
F
F
F不是积性函数,而单项式是积性函数.
所以,我们可以把
F
F
F按照定义式展开:
L
(
x
)
=
∑
i
=
1
k
+
2
ℓ
(
x
)
y
i
∏
j
≠
i
1
x
i
−
x
j
(
x
−
x
i
)
L(x)=\sum _{i=1}^{k+2} \ell(x)\dfrac {y_i\prod_{j\ne i}\dfrac 1{x_i-x_j}}{(x-x_i)}
L(x)=∑i=1k+2ℓ(x)(x−xi)yi∏j=ixi−xj1.
我们可以预处理
ℓ
(
x
)
\ell(x)
ℓ(x),然后用多项式除法
O
(
k
2
)
O(k^2)
O(k2)可以求出每一项的系数(
L
(
x
)
=
∑
i
=
0
k
+
1
a
i
x
i
L(x)=\sum_{i=0}^{k+1} a_ix^i
L(x)=∑i=0k+1aixi中的
a
i
a_i
ai.)
如果懒得写多项式除法,我们可以直接乘
O
(
k
3
)
O(k^3)
O(k3).
代入: G ( n ) = ∑ d ∣ n d k μ ( d ) F ( n d ) = ∑ d ∣ n d k μ ( d ) ∑ i = 0 k + 1 a i ( n / d ) i = ∑ i = 0 k + 1 a i n i ∑ d ∣ n d k − i μ ( d ) G(n) = \sum_{d|n} d^k\mu(d) F(\dfrac n d)=\sum_{d|n} d^k\mu(d) \sum_{i=0}^{k+1}a_i(n/d)^i=\sum_{i=0}^{k+1}a_in^i\sum_{d|n}d^{k-i}\mu(d) G(n)=∑d∣ndkμ(d)F(dn)=∑d∣ndkμ(d)∑i=0k+1ai(n/d)i=∑i=0k+1aini∑d∣ndk−iμ(d).
后面部分显然是积性函数,且 μ ( p 2 ) \mu(p^2) μ(p2)开始就=0,所以每一个质因数只有俩项.
总复杂度为 O ( k 2 + w k ) O(k^2+wk) O(k2+wk).
TP void ad(o& x,o y) {x += y; if(x >= mod) x -= mod;}
TP void dl(o& x,o y) {x -= y; if(x < 0) x += mod;}
ll k, n;
int w, p[1010];
ll power(ll a,ll b=mod-2) {
ll c=1;
while(b) {
if(b&1) c=c*a%mod;
b /= 2; a=a*a%mod;
}
return c;
}
ll f[N], L[N], a[N], y[N];
int len;
void add_poly() {
rep(i, 0, k + 2)
ad(f[i], a[i]), a[i] = 0;
}
void mul_poly(ll a, ll b) {
REP(i, 1, k + 2)
L[i] = (L[i - 1] * a + L[i] * b) % mod;
L[0] = L[0] * b % mod;
}
void div_poly(ll b) {// /(x+b)
ll u = L[k + 2];
REP(i, 0, k + 1) {
a[i] = u;
u = ((L[i] - u * b) % mod + mod) % mod;
}
}
void mult(ll *a,ll b) {
rep(i, 0, k + 2)
a[i] = a[i] * b % mod;
}
void solve() {
qr(k);
n = 1;
qr(w);
FOR(i, w) {
int x; qr(p[i]); qr(x);
n = n * power(p[i], x) % mod;
}
L[0] = 1;
FOR(i, k + 2) y[i] = (y[i - 1] + power(i, k)) % mod, mul_poly(1, mod - i);
FOR(i, k + 2) {
div_poly(mod - i);
ll s = 1;
FOR(j, k + 2) if(i ^ j) s = s * (i - j + mod) % mod;
s = y[i] * power(s) % mod;
mult(a, s); add_poly();
}
ll ans = 0, pw = 1;
rep(i, 0, k + 1) {
ll s = f[i] * pw % mod;
pw = pw * n % mod;
FOR(j, w) s = s * (1 - power(p[j], k - i + mod - 1) + mod) % mod;
ad(ans, s);
}
pr2(ans);
}
CF995F Cowmpany Cowmpensation
树形结构，给每个节点分配工资 ( [ 1 , d ] ) ([1, d]) ([1,d])，子节点不能超过父亲节点的工资，问有多少种分配方案。
定义状态
f
[
i
]
[
j
]
f[i][j]
f[i][j]表示以
i
i
i为节点的子树中
i
i
i的权值为
j
j
j的合法方案数.
则有:
f
′
[
i
]
[
j
]
=
f
[
i
]
[
j
]
∗
∑
k
=
1
j
f
[
z
]
[
k
]
(
z
∈
s
o
n
i
)
f'[i][j] = f[i][j]*\sum_{k=1}^jf[z][k](z\in son_i)
f′[i][j]=f[i][j]∗∑k=1jf[z][k](z∈soni).
由此, 得出
f
[
i
]
f[i]
f[i]的次数是子树大小.
我们只需要求一个 O ( n 2 ) O(n^2) O(n2)的 d p dp dp,然后前缀和后插值即可.
struct LA {
int n, y[N];
int& operator [] (int i) {return y[i];}
ll operator () (ll k) {
ll ans = 0;
FOR(i, n) {
ll a = y[i], b = 1;
FOR(j, n) if(i ^ j)
a = a * (k - j) % mod,
b = b * (i - j) % mod;
ans = (ans + a * power(b)) % mod;
}
return (ans + mod) % mod;
}
} g;
int fa[N];
void solve() {
qr(n); qr(D);
rep(i, 2, n) {
qr(fa[i]);
}
g.n = mx = n + 5;
FOR(i, n) fill(f[i] + 1, f[i] + mx + 1, 1);
REP(i, 1, n) {
ll b = 0;
FOR(j, mx) {
b = (b + f[i][j]) % mod;
f[fa[i]][j] = b * f[fa[i]][j] % mod;
}
}
FOR(i, mx) g[i] = f[1][i] = (f[1][i] + f[1][i - 1]) % mod;
pr2(g(D));
}
P7116 [NOIP2020] 微信步数
首先,有一个80分做法和拉格朗日无关.
我们需要确定哪一维度先碰到边界, 为此对于每一个维度我们设
f
i
r
[
i
]
fir[i]
fir[i]表示第一个到达
i
i
i的时间,容易由周期性
O
(
w
)
O(w)
O(w)预处理(
f
i
r
[
i
+
Δ
]
=
f
i
r
[
i
]
+
n
fir[i+\Delta]=fir[i]+n
fir[i+Δ]=fir[i]+n).
然后,对于一个该维度的坐标
x
x
x,其结束时间
e
d
[
x
]
min
(
f
i
r
[
−
x
]
,
w
−
x
+
1
)
ed[x]\min(fir[-x],w-x+1)
ed[x]min(fir[−x],w−x+1).
然后我们把每一个维度的 e d ed ed数组排序.
现在我们考虑哪一个维度先碰到边界.
假设我们钦定在第
x
x
x维度于第
y
y
y小的结束时间(
e
d
x
[
y
]
ed_x[y]
edx[y])弹出.
那么合法的方案,即其他维度的所有弹出时间比其大的位置的组合.
即
∏
i
≠
x
∑
j
e
d
i
[
j
]
>
e
d
x
[
y
]
\prod_{i\ne x}\sum_j ed_i[j]>ed_x[y]
∏i=x∑jedi[j]>edx[y].
如果采用二分的话,复杂度为
O
(
∑
w
log
w
)
O(\sum w \log w)
O(∑wlogw).
我们可以用
k
k
k指针法优化到
O
(
∑
w
)
O(\sum w)
O(∑w).
上面没有考虑排序的复杂度,比赛的时候直接用sort被卡到65嘤嘤嘤.赛后用个基数排序就80了.
正解部分:
设
r
i
,
j
,
l
i
,
j
r_{i,j},l_{i,j}
ri,j,li,j表示第
i
i
i维度经过
j
j
j步后最大/最小的位置变化量,
d
i
d_i
di为经过
n
n
n步后的位移量.
可以发现,经过一回合操作(
n
n
n步)后,该维度合法的起始位置为
w
i
−
(
r
i
,
n
−
l
i
,
n
)
w_i-(r_{i,n}-l_{i,n})
wi−(ri,n−li,n).
之后,在经过
t
t
t个回合,该维度合法的起始位置为
w
i
−
(
r
i
,
n
−
l
i
,
n
)
−
∣
d
i
∣
t
w_i-(r_{i,n}-l_{i,n})-|d_i|t
wi−(ri,n−li,n)−∣di∣t.
如果在走
z
z
z步的话,该维度合法起始位置为
f
i
(
t
)
=
w
i
−
(
max
(
r
i
,
n
,
r
i
,
z
+
d
)
−
min
(
l
i
,
n
,
l
i
,
z
+
d
)
−
∣
d
i
∣
t
f_i(t)=w_i-(\max(r_{i,n},r_{i,z}+d)-\min(l_{i,n},l_{i,z}+d)-|d_i|t
fi(t)=wi−(max(ri,n,ri,z+d)−min(li,n,li,z+d)−∣di∣t.
那么总体的方案即为
∏
i
=
1
k
f
i
(
t
)
\prod_{i=1}^k f_i(t)
∏i=1kfi(t).
这是一个关于
t
t
t的
k
k
k次多项式,
t
t
t的取值,从0开始,到出边界时中间的回合数
T
T
T.
这显然就是一个关于
T
T
T的
k
+
2
k+2
k+2次多项式.
插值一下就好了.
总体复杂度为
O
(
n
k
2
)
O(nk^2)
O(nk2).
int n, k;
struct rec {
int l[N], r[N], d, w, A;
inline void upd(int now) {
cmin(l[now] = l[now - 1], d);
cmax(r[now] = r[now - 1], d);
}
bool not_out() {
return !d && r[n] - l[n] + 1 <= w;
}
int maxt(int now) {
return !d ? 1e9: floor(1.0 * get(now, 0) / abs(d));
}
ll get(int i,int t) {// 中间有t个整周期, 再加i步.
return w - (max(r[n], r[i] + d) - min(l[n], l[i] + d) + t * abs(d));
}
} A[11];
ll power(ll a, ll b = mod - 2) {
ll c = 1;
while(b) {
if(b & 1) c = c * a % mod;
b /= 2; a = a * a % mod;
}
return c;
}
void ad(ll &x,ll y) {x += y; if(x >= mod) x -= mod;}
void dl(ll &x,ll y) {x -= y; if(x < 0) x += mod;}
struct Poly {
ll f[15]; int len;
Poly(int x = 0) {memset(f, 0, sizeof f); f[0] = x; len = 0;}
ll& operator [](int i) {return f[i];}
ll operator () (ll k) {
ll ans = 0;
REP(i, 0, len) ans = (ans * k + f[i]) % mod;
return ans;
}
Poly mul(ll a,ll b) {// * (ax + b)
Poly res; res.len = len + 1;
rep(i, 0, len) ad(res[i], f[i] * b % mod), ad(res[i + 1], f[i] * a % mod);
return res;
}
Poly div(ll a,ll b) {// / (ax + b)
Poly res; res.len = len - 1;
ll u = f[len]; a = power(a);
REP(i, 0, len - 1) {
res[i] = u * a % mod;
dl(u = f[i], res[i] * b % mod);
}
return res;
}
} ans(1);
ll jc[15], inv[15];
void init() {
jc[0] = inv[0] = 1;
FOR(i, 14) jc[i] = jc[i - 1] * i % mod, inv[i] = power(jc[i]);
}
struct Lagrange {
int n;
ll y[15], pre[15], suf[15];
ll& operator [](int i) {return y[i];}
ll operator() (ll k) {
if(k <= n) return y[k];
ll ans = 0;
pre[0] = k; suf[n + 1] = 1;
FOR(i, n) pre[i] = pre[i - 1] * (k - i) % mod;
REP(i, 0, n) suf[i] = suf[i + 1] * (k - i) % mod;
rep(i, 0, n)
ans = (ans + y[i] * (i ? pre[i - 1] : 1) % mod * suf[i + 1] % mod
* inv[i] % mod * inv[n - i] % mod * ((n - i) & 1 ? -1 : 1)) % mod;
return (ans + mod) % mod;
}
} L;
int c[N];
void solve() {
qr(n); qr(k);
FOR(i, k) qr(A[i].w);
FOR(i, n) {
qr(c[i]); int d; qr(d);
A[c[i]].d += d;
FOR(j, k) A[j].upd(i);
}
bool flag = 1;
FOR(i, k)
flag &= A[i].not_out(),
ans = ans.mul(mod - abs(A[i].d), A[i].get(0, 0));
if(flag) exit(puts("-1") & 0);
init();
ll sum = 0; L.n = k + 2;
rep(i, 0, n - 1) {
int o = c[i];
if(A[o].r[i] - A[o].l[i] > A[o].w) break;
ll max_t = INF;
FOR(j, k) cmin(max_t, (ll)A[j].maxt(i));
//t + 2轮, t >= 0 的情况.
if(max_t >= 0) {
L[0] = ans(0);
FOR(j, L.n)
ad(L[j] = L[j - 1], ans(j));
ad(sum, L(max_t));
}
//只有一回合的情况.
ll s = 1;
FOR(j, k)
s = s * (A[j].w - (A[j].r[i] - A[j].l[i])) % mod;
ad(sum, s);
o = c[i + 1];
if(A[o].get(i,0) != A[o].get(i + 1, 0))
ans = ans.div(mod - abs(A[o].d), A[o].get(i, 0)),
ans = ans.mul(mod - abs(A[o].d), A[o].get(i + 1, 0));
}
pr2(sum);
}
1221
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
