Math
expectation and posibility
一些前置芝士
概率
P ( A ∪ B ) = P ( A ) + P ( B ) − P ( A ∩ B ) P ( B − A ) = P ( B ) − P ( A ) ( P ( A ) ≤ P ( B ) ) P ( B − A ) = P ( B ) − P ( A ∩ B ) ( ∀ A , B ) P ( A ∣ B ) = P ( A ∩ B ) P ( B ) P ( A B ) = P ( A ) × P ( B ∣ A ) = P ( B ) × P ( A ∣ B ) ( A , B 有 关 联 ) \begin{aligned} & P(A\cup B) = P(A) + P(B) - P(A \cap B) \\ & P(B-A) = P(B) - P(A) \qquad (P(A) \leq P(B)) \\ & P(B-A) = P(B) - P(A \cap B) \qquad (\forall A, B) \\ & P(A\mid B) = \frac{P(A \cap B)}{P(B)} \\ & P(AB) = P(A) \times P(B\mid A) = P(B) \times P(A\mid B) \quad (A,B 有关联) \\ \end{aligned} P(A∪B)=P(A)+P(B)−P(A∩B)P(B−A)=P(B)−P(A)(P(A)≤P(B))P(B−A)=P(B)−P(A∩B)(∀A,B)P(A∣B)=P(B)P(A∩B)P(AB)=P(A)×P(B∣A)=P(B)×P(A∣B)(A,B有关联)
全概率公式
P ( A ) = ∑ i = 1 n P ( A ∣ H i ) × P ( H i ) P(A) = \sum_{i=1}^n P(A\mid H_i) \times P(H_i) P(A)=i=1∑nP(A∣Hi)×P(Hi)
贝叶斯定理
P ( A ∣ B ) = P ( A ) P ( B ∣ A ) P ( B ) P(A\mid B) = \frac{P(A)P(B\mid A)}{P(B)} P(A∣B)=P(B)P(A)P(B∣A)
期望
E ( x ) = ∑ i p i x i E ( C ) = C E ( C x ) = C E ( x ) E ( x + y ) = E ( x ) + E ( y ) E ( x y ) = E ( x ) E ( y ) ( x , y 相 互 独 立 ) \begin{aligned} & E(x) = \sum_i p_ix_i \\ & E(C) = C \quad E(Cx) = CE(x) \quad E(x+y) = E(x) + E(y) \\ & E(xy) = E(x)E(y) \qquad (x,y相互独立) \end{aligned} E(x)=i∑pixiE(C)=CE(Cx)=CE(x)E(x+y)=E(x)+E(y)E(xy)=E(x)E(y)(x,y相互独立)
全期望公式
E ( Y ) = E ( E ( Y ∣ X ) ) = ∑ i P ( x = x i ) E ( Y ∣ X = x i ) E(Y) = E(E(Y \mid X)) = \sum_i P(x=x_i)E(Y \mid X = x_i) E(Y)=E(E(Y∣X))=i∑P(x=xi)E(Y∣X=xi)
luogu1291 SHOI2002百事世界杯之旅
(2021.11.3)
“……在 2002 2002 2002 年 6 6 6 月之前购买的百事任何饮料的瓶盖上都会有一个百事球星的名字。只要凑齐所有百事球星的名字,就可参加百事世界杯之旅的抽奖活动,获得球星背包,随声听,更可赴日韩观看世界杯。还不赶快行动!”
你关上电视,心想:假设有 n n n 个不同的球星名字,每个名字出现的概率相同,平均需要买几瓶饮料才能凑齐所有的名字呢?
设
E
i
,
j
E_{i, j}
Ei,j 表示总共有 i 个,还有 j 个没有抽到的期望。,我们考虑下一次抽取:有
k
n
\frac k n
nk 的概率抽到我们还没有的,有
n
−
k
n
\frac{n-k}n
nn−k 的概率抽到我们已经抽到的,于是我们可以得到以下的递推式:
E
n
,
k
=
k
n
E
n
,
k
−
1
+
n
−
k
n
E
n
,
k
+
1
E_{n, k} = \frac k n E_{n,k-1} + \frac{n-k}n E_{n, k} + 1
En,k=nkEn,k−1+nn−kEn,k+1
然后稍微移项化简一下,就会变成这样:
E
n
,
k
=
k
n
E
n
,
k
−
1
+
n
−
k
n
E
n
,
k
+
1
E
n
,
k
−
n
−
k
n
E
n
,
k
=
k
n
E
n
,
k
=
k
n
E
n
,
k
−
1
+
1
E
n
,
k
=
E
n
,
k
−
1
+
n
k
\begin{aligned} & E_{n, k} = \frac k n E_{n,k-1} + \frac{n-k}n E_{n, k} + 1 \\ & E_{n, k} - \frac{n-k}n E_{n, k} = \frac k n E_{n,k} = \frac k n E_{n,k-1} + 1 \\ & E_{n, k} = E_{n, k-1} + \frac n k \end{aligned}
En,k=nkEn,k−1+nn−kEn,k+1En,k−nn−kEn,k=nkEn,k=nkEn,k−1+1En,k=En,k−1+kn
又因为我们知道
E
1
,
1
=
1
E_{1,1} = 1
E1,1=1,,所以我们得到通项公式:
E
(
x
)
=
n
∑
i
=
1
n
1
i
E(x) = n\sum_{i=1}^n\frac 1 i
E(x)=ni=1∑ni1
骰子
(2021.11.3)
一个骰子,扔 n 次的点数之和大于 x 的概率。
首先我们考虑算出扔 n 次点数之和大于 x 的方案数,最后再除以
6
n
6^n
6n 即可。我们设
f
i
,
j
f_{i, j}
fi,j 表示扔了 i 次,点数之和为 j 的方案数。我们就有:
f
i
,
j
=
∑
k
=
1
6
f
i
−
1
,
j
−
k
f_{i, j} = \sum_{k=1}^6 f_{i-1,j-k}
fi,j=k=1∑6fi−1,j−k
所以整个算法就是先 DP 然后最后的答案就是:
a n s = ∑ i = x 6 n f n , i 6 n ans = \frac{\sum\limits_{i = x}^{6n}f_{n,i}}{6^n} ans=6ni=x∑6nfn,i
OX 序列
(2022.2.6)
给定序列,每个位置是 o o o 的概率为 p i p_i pi,为 x x x 的概率是 1 − p i 1 - p_i 1−pi。对于一个 o x ox ox 序列,连续 a a a 长度的 o o o 会得到 a 3 a^3 a3 的收益,问最终得到的 o x ox ox 序列的期望收益是多少。
我们设 f i f_i fi 为某一种具体的放置方法中已经放到第 i i i 位的收益。 x i x_i xi 表示某一种具体的放置方法中从最后一位开始连续的 o o o 的长度。那么我么需要求解的就是 E ( f n ) E(f_n) E(fn)。首先我们很容易知道一个关于 f f f 的递推式:
f i = f i − 1 + ( x i − 1 + 1 ) 3 − x i − 1 3 = f i − 1 + 3 x i − 1 2 + 3 x i − 1 + 1 f_i = f_{i-1} + (x_{i-1} + 1)^3 - x_{i-1}^3 = f_{i-1} + 3x_{i-1}^2 + 3x_{i-1} + 1 fi=fi−1+(xi−1+1)3−xi−13=fi−1+3xi−12+3xi−1+1
那么:
E ( f i ) = E ( f i − 1 ) + 3 E ( x i − 1 2 ) + 3 E ( x i − 1 ) + 1 E(f_i) = E(f_{i-1}) + 3E(x_{i-1}^2) + 3E(x_{i-1}) + 1 E(fi)=E(fi−1)+3E(xi−12)+3E(xi−1)+1
化成这样之后,我们就可以发先如果我们求出 3 E ( x i − 1 2 ) + 3 E ( x i − 1 ) 3E(x_{i-1}^2) + 3E(x_{i-1}) 3E(xi−12)+3E(xi−1),这就是一个关于 E ( f ) E(f) E(f) 的递推式。就可以直接求出来。所以我们考虑怎样求出 E ( x i ) E(x_{i}) E(xi) 和 E ( x i 2 ) E(x_{i}^2) E(xi2)。首先处理一个简单一点的:
E ( x i ) = ? \begin{aligned} E(x_{i}) = ? \end{aligned} E(xi)=?
要知道这个我们首先需要知道 x i x_i xi 的递推式:
{ 第 i 位 o , 概 率 为 p i , x i = x i − 1 + 1 第 i 位 x , 概 率 为 1 − p i , x i = 0 \begin{cases} 第i位o, \; 概率为p_i, \; x_i = x_{i-1}+1 \\ 第i位x, \; 概率为1-p_i, \; x_i = 0 \end{cases} {第i位o,概率为pi,xi=xi−1+1第i位x,概率为1−pi,xi=0
那么:
E ( x i ) = p i × E ( x i − 1 + 1 ) + E ( 1 − p i ) × 0 = p i E ( x i − 1 ) + p i \begin{aligned} E(x_i) = p_i \times E(x_{i-1} + 1 ) + E(1-p_i) \times 0 = p_iE(x_{i-1}) + p_i \end{aligned} E(xi)=pi×E(xi−1+1)+E(1−pi)×0=piE(xi−1)+pi
这样之后我们也可以用递推解决 E ( x i ) E(x_i) E(xi) 的值。然后考虑 E ( x i 2 ) E(x_i^2) E(xi2) 的值。首先:
{ 第 i 位 o , 概 率 为 p i , x i 2 = ( x i − 1 + 1 ) 2 = x i − 1 2 + 2 x i − 1 + 1 第 i 位 x , 概 率 为 1 − p i , x i 2 = 0 \begin{cases} 第i位o, \; 概率为p_i, \; x_i^2 = (x_{i-1}+1)^2 = x_{i-1}^2 + 2x_{i-1} + 1 \\ 第i位x, \; 概率为1-p_i, \; x_i^2 = 0 \end{cases} {第i位o,概率为pi,xi2=(xi−1+1)2=xi−12+2xi−1+1第i位x,概率为1−pi,xi2=0
那么:
E ( x i 2 ) = p i E ( x i − 1 2 + 2 x i − 1 + 1 ) + ( 1 − p i ) × 0 = p i E ( x i − 1 2 ) + 2 p i E ( x i − 1 ) + p i \begin{aligned} E(x_i^2) = & p_iE(x_{i-1}^2 + 2x_{i-1} + 1) + (1-p_i) \times 0 \\ = & p_i E(x_{i-1}^2) + 2p_iE(x_{i-1} ) + p_i \end{aligned} E(xi2)==piE(xi−12+2xi−1+1)+(1−pi)×0piE(xi−12)+2piE(xi−1)+pi
显然这也是一个递推式,那我们就可以 O ( n ) O(n) O(n) 求出最终答案了。
luogu4316 绿豆蛙的归宿
(2022.4.8)
给定一张 m m m 条边 n n n 个点的,起点为 1 1 1 终点为 n n n 的有向无环图。现在在这张图上走,到达一个顶点的时候(假设这个点的度为 k k k)可以选择一条路径离开这个点,并且走向每条边的概率是 1 k \frac 1k k1。问你从 1 1 1 到 n n n 的期望路径总长。
显然现在给的是一张 D A G DAG DAG。然后考虑期望概率 d p dp dp。我们设 f i f_i fi 表示 i i i 点到终点的期望路径总长。显然我们要求的答案就是 f 1 f_1 f1,且有 f n = 0 f_n = 0 fn=0。然后就是 d p dp dp 方程了
因为我们知道了终点的状态的值,所以我们考虑逆推,我们就能得到:
f x = ∑ ( x , y ) ∈ G x f y + w ( x , y ) ∣ G x ∣ f_x = \sum_{(x, y)\in G_x}\frac{f_y + w_{(x, y)}}{|G_x|} fx=(x,y)∈Gx∑∣Gx∣fy+w(x,y)
转移的过程就可以通过连反向边然后沿着拓扑序走一遍转移就行了。时间复杂度就是 O ( n + m ) O(n + m) O(n+m)。
number theory
一些前置芝士
欧拉定理:
∀
g
c
d
(
a
,
n
)
=
1
→
a
φ
(
n
)
≡
1
(
m
o
d
n
)
\forall gcd(a, n) = 1 \rightarrow a^{\varphi(n)} \equiv 1 \pmod n
∀gcd(a,n)=1→aφ(n)≡1(modn)
欧拉定理的推论:
∀ g c d ( a , n ) = 1 → a b m o d φ ( n ) ≡ 1 ( m o d n ) \forall gcd(a, n) = 1 \rightarrow a^{b \mod \varphi(n)} \equiv 1 \pmod n ∀gcd(a,n)=1→abmodφ(n)≡1(modn)
费马小定理:
∀ p 是 质 数 → a p − 1 ≡ ( m o d p ) \forall p 是质数 \rightarrow a^{p - 1} \equiv \pmod p ∀p是质数→ap−1≡(modp)
exgcd():
int exgcd(int a, int b, int &x, int &y){
if(b == 0) { x = 1; y = 0; return a; }
int d = exgcd(b, a % b, x, y);
int z = x; x = y; y = z - y * (a / b);
return d;
}
中国剩余定理:
( i : 1 → n ) { x ≡ a i ( m o d m i ) m = ∏ i = 1 n m i M i = m m i M i t i ≡ 1 ( m o d m i ) → x = ∑ i = 1 n a i M i t i (i:1 \to n) \;\;\; \begin{cases} x \equiv a_i \pmod {m_i} \\ m = \prod\limits_{i = 1}^nm_i \\ M_i = \frac{m}{m_i} \\ M_it_i \equiv 1 \pmod {m_i} \end{cases} \rightarrow x = \sum_{i=1}^na_iM_it_i (i:1→n)⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧x≡ai(modmi)m=i=1∏nmiMi=mimMiti≡1(modmi)→x=i=1∑naiMiti
L u c a s Lucas Lucas 定理(把 n n n 和 m m m 表示成 p p p 进制数,然后每一位计算组合数再乘起来):
p 是 质 数 → C n m ≡ C n m o d p m m o d p × C n p m p ( m o d p ) p 是质数 \rightarrow C_n^m \equiv C_{n\mod p}^{m\mod p} \times C_{\frac np}^{\frac mp} \pmod p p是质数→Cnm≡Cnmodpmmodp×Cpnpm(modp)
埃氏筛
埃氏筛的思想很简单,就是从前往后,遇到一个质数就把所有后面小于 n 的这个质数的倍数全都标记为合数。后面遇到他们就直接跳过。
筛质数
按照上面的思想直接模拟就好了,时间复杂度 O ( n log log n ) O(n \log \log n) O(nloglogn):
#include<bits/stdc++.h>
using namespace std;
#define MAXN 100100
int v[MAXN] = { 0 }; // 合数标记
void primes(int n){
memset(v, 0, sizeof(v));
for(int i = 2; i <= n; i++){
if(v[i]) continue;
cout << i << endl;
for(int j = 1; j <= n / i; j++) v[i * j] = 1;
}
}
int main(){
primes(100);
return 0;
}
筛欧拉函数
利用欧拉函数的计算式和埃氏筛,能在 O ( n log n ) O(n \log n) O(nlogn) 的时间内筛出欧拉函数的值。就是首先初始化所有的 phi[i] = 1,然后再枚举从 2 到 n 的每个数,用埃氏筛找出质数,然后把后面质数的倍数根据计算式计算就好了:
φ ( n ) = n ∏ p ∣ n p − 1 p \varphi(n) = n \prod_{p | n} \frac {p-1}p φ(n)=np∣n∏pp−1
#include<bits/stdc++.h>
using namespace std;
#define MAXN 100100
int phi[MAXN] = { 0 };
void euler(int n){
for(int i = 1; i <= n; i++) phi[i] = i;
for(int i = 2; i <= n; i++){
if(phi[i] == i)
for(int j = i; j <= n; j += i)
phi[j] = phi[j] / i * (i - 1);
}
}
int main(){
euler(20);
for(int i = 1; i <= 20; i++){
cout << "phi[" << i << "] = " << phi[i] << endl;
}
return 0;
}
筛莫比乌斯函数
首先把所有的 mu 初始化为 1,接下来对于每一个质数 p 把 mu[p] 变成 -1。并扫描 p 的倍数们看他们能否被 p 2 p^2 p2 整除,如果可以 mu[x] = 0,如果不行 mu[x] = -mu[x]。
#include<bits/stdc++.h>
using namespace std;
#define MAXN 100100
int v[MAXN] = { 0 };
int mu[MAXN] = { 0 };
void mobius(int n){
memset(v, 0, sizeof(v));
for(int i = 2; i <= n; i++) mu[i] = 1;
for(int i = 2; i <= n; i++){
if(v[i]) continue;
mu[i] = -1;
for(int j = 2 * i; j <= n; j += i){
v[j] = 1;
if((j / i) % i == 0) mu[j] = 0;
else mu[j] *= -1;
}
}
}
int main(){
mobius(20);
for(int i = 1; i <= 20; i++){
cout << "mu[" << i << "] = " << mu[i] << endl;
}
return 0;
}
线性筛
顾名思义,这玩意儿能在线性的时间内筛出质数。具体的做法是优化一下埃氏筛,在埃氏筛里面,一个合数可能被很多个质数标记过,比如 12 就被 2 和 3 同时标记过。这样就会降低我们的效率。线性筛的想法就是让所有的合数只有一种被筛出来的方式。也就是被它的最小质因子筛出来。具体做法如下:
- 枚举 2 ∼ n 2 \sim n 2∼n 的所有数,考虑用一下方式维护一个数组 v(表示这个数的最小质因子)。
- 如果 v[i] = i,那么这个是是质数,把它存下来。
- 扫描不大于 v[i] 的所有质数 p,令 v[ip] = p。
筛质数
#include<bits/stdc++.h>
using namespace std;
#define MAXN 100100
int v[MAXN] = { 0 };
int prime[MAXN] = { 0 };
int primes(int n){
memset(v, 0, sizeof(v));
int cnt = 0; // 质数的个数
for(int i = 2; i <= n; i++){
if(!v[i]){
v[i] = i; prime[++cnt] = i;
}
for(int j = 1; j <= cnt; j++){ // 枚举所有质数
if(prime[j] > v[i] or prime[j] > n / i) break; // 质数比它大或者超过范围
v[i * prime[j]] = prime[j]; // 更新最小质因子
}
}
return cnt;
}
int main(){
int num = primes(100);
for(int i = 1; i <= num; i++){
cout << prime[i] << endl;
}
return 0;
}
筛欧拉函数
这就要用到几个欧拉函数的性质了:
- 如果 p ∣ n p \mid n p∣n 且 p 2 ∣ n p^2 \mid n p2∣n 则有 φ ( n ) = φ ( n p ) × p \varphi(n) = {\varphi(\frac np)} \times p φ(n)=φ(pn)×p
- 如果 p ∣ n p \mid n p∣n 且 p 2 ∤ n p^2 \nmid n p2∤n 则有 φ ( n ) = φ ( n p ) × ( p − 1 ) \varphi(n) = \varphi(\frac np) \times (p-1) φ(n)=φ(pn)×(p−1)
我们发现这两个判断和线性筛的两个判断很类似,所以我们可以利用线性筛,从 φ ( n p ) \varphi(\frac np) φ(pn) 递推到 φ ( n ) \varphi(n) φ(n)
#include<bits/stdc++.h>
using namespace std;
#define MAXN 100100
int v[MAXN] = { 0 };
int phi[MAXN] = { 0 };
void euler(int n){
memset(v, 0, sizeof(v));
int cnt = 0;
for(int i = 2; i <= n; i++){
if(!v[i]){
v[i] = i; prime[++cnt] = i;
phi[i] = i - 1;
}
for(int j = 1; j <= cnt; j++){
if(prime[j] > v[i] or prime[j] > n / i) break;
v[i *prime[j]] = prime[j];
if(i % prime[j]) phi[i * prime[j]] = phi[i] * (prime[j] - 1);
else phi[i * prime[j]] = phi[i] * prime[j];
}
}
}
int main(){
euler(20);
for(int i = 1; i <= 20; i++){
cout << "phi[" << i << "] = " << phi[i] << endl;
}
return 0;
}
luogu2185 SDOI2008 仪仗队
(2021.11.7)
有一个点阵(方阵),你站在左下角,问你能总共看到几个点。然后我们可以发现,如果我们将左下角的点放在坐标轴 原点,能被看见的点的横纵坐标一定是互质的(除了两个和它相邻的点)。所以我们要求解的问题的答案就是:
a
n
s
=
2
+
∑
i
=
1
n
∑
j
=
1
n
[
g
c
d
(
i
,
j
)
=
1
]
ans = 2 + \sum_{i=1}^n \sum_{j=1}^n [gcd(i, j) = 1]
ans=2+i=1∑nj=1∑n[gcd(i,j)=1]
然后可以对这个式子进行化简:
a
n
s
−
2
=
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
[
g
c
d
(
i
,
j
)
=
1
]
=
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
ε
(
g
c
d
(
i
,
j
)
)
=
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
∑
d
∣
g
c
d
(
i
,
j
)
μ
(
d
)
=
∑
d
=
1
n
−
1
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
[
d
∣
g
c
d
(
i
,
j
)
]
μ
(
d
)
=
∑
d
=
1
n
−
1
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
[
d
∣
i
∧
d
∣
j
]
μ
(
d
)
=
∑
d
=
1
n
−
1
μ
(
d
)
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
[
d
∣
i
∧
d
∣
j
]
=
∑
d
=
1
n
−
1
μ
(
d
)
⌊
n
−
1
d
⌋
2
\begin{aligned} ans-2 = & \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} [gcd(i, j) = 1] \\ = & \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} \varepsilon(gcd(i, j)) = \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} \sum_{d \mid gcd(i, j)}\mu(d) \\ = & \sum_{d = 1}^{n-1} \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} [d \mid gcd(i, j)]\mu(d) = \sum_{d = 1}^{n-1} \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} [d \mid i \wedge d \mid j]\mu(d) \\ = & \sum_{d=1}^{n-1} \mu(d) \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} [d \mid i \wedge d \mid j] \\ = & \sum_{d=1}^{n-1} \mu(d) \left\lfloor \frac {n-1}d \right\rfloor^2 \end{aligned}
ans−2=====i=1∑n−1j=1∑n−1[gcd(i,j)=1]i=1∑n−1j=1∑n−1ε(gcd(i,j))=i=1∑n−1j=1∑n−1d∣gcd(i,j)∑μ(d)d=1∑n−1i=1∑n−1j=1∑n−1[d∣gcd(i,j)]μ(d)=d=1∑n−1i=1∑n−1j=1∑n−1[d∣i∧d∣j]μ(d)d=1∑n−1μ(d)i=1∑n−1j=1∑n−1[d∣i∧d∣j]d=1∑n−1μ(d)⌊dn−1⌋2
然后 μ ( d ) \mu(d) μ(d) 可以用埃氏筛预处理然后 O ( 1 ) O(1) O(1) 查询。后面的 ⌊ n d ⌋ 2 \left\lfloor \frac nd \right\rfloor^2 ⌊dn⌋2 显然也是可以 O ( 1 ) O(1) O(1) 查询的,那么总的复杂度就是 O ( n log n + n ) O(n\log n + n) O(nlogn+n)
P2522 [HAOI2011]Problem b
(2021.11.7)
对于给出的 n 个询问,每次求出有多少个数对 ( x , y ) (x, y) (x,y) 满足 a ≤ x ≤ b , c ≤ y ≤ d a \leq x \leq b, c \leq y \leq d a≤x≤b,c≤y≤d 且 g c d ( x , y ) = k gcd(x, y) = k gcd(x,y)=k。
所以这个题的答案就是:
a
n
s
=
∑
i
=
a
b
∑
j
=
c
d
[
g
c
d
(
i
,
j
)
=
k
]
ans = \sum_{i = a}^b \sum_{j = c}^d [gcd(i, j) = k]
ans=i=a∑bj=c∑d[gcd(i,j)=k]
这个式子就和上面那个题的试着长得很像了,但是这个上下界看着还是有点不舒服,因为我们想把它转化成下界为 1 的求和,于是我们利用容斥原理,转化一下问题:
a
n
s
=
∑
i
=
a
b
∑
j
=
c
d
[
g
c
d
(
i
,
j
)
=
k
]
=
∑
i
=
1
b
∑
j
=
1
d
[
g
c
d
(
i
,
j
)
=
k
]
−
∑
i
=
1
a
−
1
∑
j
=
1
d
[
g
c
d
(
i
,
j
)
=
k
]
−
∑
i
=
1
b
∑
j
=
1
c
−
1
[
g
c
d
(
i
,
j
)
=
k
]
+
∑
i
=
1
a
−
1
∑
j
=
1
c
−
1
[
g
c
d
(
i
,
j
)
=
k
]
\begin{aligned} ans = & \sum_{i = a}^b \sum_{j = c}^d [gcd(i, j) = k] \\ = & \sum_{i = 1}^b \sum_{j = 1}^d [gcd(i, j) = k] - \sum_{i = 1}^{a-1} \sum_{j = 1}^d [gcd(i, j) = k] - \sum_{i = 1}^b \sum_{j = 1}^{c-1} [gcd(i, j) = k] + \sum_{i = 1}^{a-1} \sum_{j = 1}^{c-1} [gcd(i, j) = k] \end{aligned}
ans==i=a∑bj=c∑d[gcd(i,j)=k]i=1∑bj=1∑d[gcd(i,j)=k]−i=1∑a−1j=1∑d[gcd(i,j)=k]−i=1∑bj=1∑c−1[gcd(i,j)=k]+i=1∑a−1j=1∑c−1[gcd(i,j)=k]
于是问题就被转化成求 4 个和上道题差不多的式子。由上一道题我们知道:
∑
i
=
1
n
∑
j
=
1
m
[
g
c
d
(
i
,
j
)
=
k
]
=
∑
i
=
1
n
∑
j
=
1
m
ε
(
g
c
d
(
i
,
j
)
k
)
=
∑
i
=
1
n
∑
j
=
1
m
∑
p
∣
g
c
d
(
i
,
j
)
/
k
μ
(
p
)
=
∑
p
=
1
m
i
n
(
n
,
m
)
∑
i
=
1
n
∑
j
=
1
m
[
p
∣
g
c
d
(
i
,
j
)
k
]
μ
(
p
)
=
∑
p
=
1
m
i
n
(
n
,
m
)
μ
(
p
)
∑
i
=
1
n
∑
j
=
1
m
[
p
∣
g
c
d
(
i
,
j
)
k
]
=
∑
p
=
1
m
i
n
(
n
,
m
)
μ
(
p
)
∑
i
=
1
n
∑
j
=
1
m
[
p
∣
i
k
∧
p
∣
j
k
]
=
∑
p
=
1
m
i
n
(
n
,
m
)
μ
(
p
)
∑
i
=
1
n
∑
j
=
1
m
[
k
p
∣
i
∧
k
p
∣
j
]
=
∑
p
=
1
m
i
n
(
n
,
m
)
μ
(
p
)
⌊
n
k
p
⌋
⌊
m
k
p
⌋
\begin{aligned} & \sum_{i = 1}^n \sum_{j = 1}^m [gcd(i, j) = k]\\ = & \sum_{i = 1}^n \sum_{j = 1}^m \varepsilon\left( \frac{gcd(i, j)}{k} \right)= \sum_{i = 1}^n \sum_{j = 1}^m \sum_{p | gcd(i, j) / k} \mu(p) = \sum_{p = 1}^{min(n, m)}\sum_{i = 1}^n \sum_{j = 1}^m \left[ p \Big| \frac{gcd(i, j)}{k} \right]\mu(p) \\ = &\sum_{p = 1}^{min(n, m)}\mu(p) \sum_{i = 1}^n \sum_{j = 1}^m \left[ p \Big| \frac{gcd(i, j)}{k} \right] = \sum_{p = 1}^{min(n, m)}\mu(p) \sum_{i = 1}^n \sum_{j = 1}^m \left[ p \Big| \frac{i}{k} \wedge p \Big| \frac{j}{k} \right] = \sum_{p = 1}^{min(n, m)}\mu(p) \sum_{i = 1}^n \sum_{j = 1}^m \left[ kp \mid i \wedge kp \mid j \right] \\ = & \sum_{p = 1}^{min(n, m)}\mu(p) \left\lfloor \frac{n}{kp} \right\rfloor \left\lfloor \frac{m}{kp} \right\rfloor \end{aligned}
===i=1∑nj=1∑m[gcd(i,j)=k]i=1∑nj=1∑mε(kgcd(i,j))=i=1∑nj=1∑mp∣gcd(i,j)/k∑μ(p)=p=1∑min(n,m)i=1∑nj=1∑m[p∣∣∣kgcd(i,j)]μ(p)p=1∑min(n,m)μ(p)i=1∑nj=1∑m[p∣∣∣kgcd(i,j)]=p=1∑min(n,m)μ(p)i=1∑nj=1∑m[p∣∣∣ki∧p∣∣∣kj]=p=1∑min(n,m)μ(p)i=1∑nj=1∑m[kp∣i∧kp∣j]p=1∑min(n,m)μ(p)⌊kpn⌋⌊kpm⌋
然后我们只需要写一个函数来处理这个,然后每次询问调用 4 次就好了。计算这个式子的时间复杂度是 O ( m i n ( n , m ) log m i n ( n , m ) + m i n ( n , m ) ) O(min(n, m) \log min(n, m) + min(n, m)) O(min(n,m)logmin(n,m)+min(n,m)),所以如果我还要处理 5 × 1 0 4 5 \times 10^4 5×104 组这样的数据,就会超时,直接计算能拿 30 分。
下面我们考虑怎么样优化处理这个式子的时间。首先看看这个式子的样子,他长得是不是就很像整除分块,那就直接整除分块计算这个式子,复杂度就被降到了 O ( m i n ( n , m ) ) O(\sqrt{min(n, m)}) O(min(n,m)),·这样一来就能过了。
luogu2480 古代猪文
(2022.4.2)
给你 q , n ∈ [ 1 , 1 0 9 ] q, n \in [1, 10^9] q,n∈[1,109],算 q ∑ d ∣ n C n d m o d p q^{\sum\limits_{d | n}C_n^d} \mod p qd∣n∑Cndmodp, 其中 p = 999911659 p = 999911659 p=999911659。
如果 q = p q = p q=p,那么上式结果为 0 0 0。否则,因为 p p p 是质数,所以 g c d ( q , n ) gcd(q, n) gcd(q,n) 互质。所以根据欧拉定理的推论,我们就能得到:
q ∑ d ∣ n C n d ≡ q ∑ d ∣ n C n d m o d ( p − 1 ) ( m o d p ) q^{\sum\limits_{d | n}C_n^d} \equiv q^{\sum\limits_{d | n}C_n^d \mod (p-1)} \pmod p qd∣n∑Cnd≡qd∣n∑Cndmod(p−1)(modp)
所以我们只需要知道怎么计算 ∑ d ∣ n C n d m o d ( p − 1 ) \sum\limits_{d|n}C_n^d \mod (p-1) d∣n∑Cndmod(p−1) 就可以了。经过观察我们发现 p − 1 = 2 × 3 × 4679 × 35617 p - 1 = 2 \times 3 \times 4679 \times 35617 p−1=2×3×4679×35617,所以它是一个 s q u a r e − f r e e − n u m b e r square-free-number square−free−number(就是所有质因子的指数都是 1 1 1 的数 qwq)。所以在这里我们可以枚举 n n n 的因数 d d d 然后用 L u c a s Lucas Lucas 定理来计算 C n d m o d ( p − 1 ) C_n^d \mod (p-1) Cndmod(p−1),然后分别计算出 ∑ d ∣ n C n d \sum\limits_{d|n}C_n^d d∣n∑Cnd 对 2 , 3 , 4679 , 35617 2, 3, 4679, 35617 2,3,4679,35617 四个质数取模的结果。我们把他们记为 a 1 , a 2 , a 3 , a 4 a_1, a_2, a_3, a_4 a1,a2,a3,a4。
然后我们解出下面这个方程组的根 x x x 就是 ∑ d ∣ n C n d m o d ( p − 1 ) \sum\limits_{d|n}C_n^d \mod (p-1) d∣n∑Cndmod(p−1) 的最小非负整数解。
{ x ≡ 2 ( m o d a 1 ) x ≡ 3 ( m o d a 2 ) x ≡ 4679 ( m o d a 3 ) x ≡ 35617 ( m o d a 4 ) \begin{cases} x \equiv 2 \pmod {a_1} \\ x \equiv 3 \pmod {a_2} \\ x \equiv 4679 \pmod {a_3} \\ x \equiv 35617 \pmod {a_4} \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧x≡2(moda1)x≡3(moda2)x≡4679(moda3)x≡35617(moda4)
解这个方程很显然就是中国剩余定理。然后我们就用快速幂求 q x q^x qx 就是答案了。
linear algebra
Power of Matrix
(2022.4.30)
给一个 n × n n \times n n×n 的矩阵 A A A,给定一个整数 k k k 求:
a n s = ∑ i = 1 k A i ans = \sum_{i = 1}^k A_i ans=i=1∑kAi
我们将这个式子看成递推也就是:
F k = ∑ i = 1 k A i F_k = \sum_{i = 1}^k A^i Fk=i=1∑kAi
其中 F 0 = O = [ 0 0 ⋯ 0 0 0 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ 0 ] F_0 = O = \begin{bmatrix}0 & 0 & \cdots& 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{bmatrix} F0=O=⎣⎢⎢⎢⎡00⋮000⋮0⋯⋯⋱⋯00⋮0⎦⎥⎥⎥⎤, F 1 = A F_1 = A F1=A,然后我们可以得到矩阵版的递推关系:
[ F k E ] = [ F k − 1 E ] × S = [ F 0 E ] × S k \begin{bmatrix} F_k & E \end{bmatrix} = \begin{bmatrix} F_{k - 1} & E \end{bmatrix} \times S = \begin{bmatrix} F_0 & E \end{bmatrix} \times S^{k} [FkE]=[Fk−1E]×S=[F0E]×Sk
上面的式子中的 E E E 表示单位矩阵,也就是:
E = [ 1 0 ⋯ 0 0 1 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ 1 ] E = \begin{bmatrix}1 & 0 & \cdots& 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix} E=⎣⎢⎢⎢⎡10⋮001⋮0⋯⋯⋱⋯00⋮1⎦⎥⎥⎥⎤
使用单位矩阵是为了更方便的写出 S S S 来,我们设:
S = [ a b c d ] S = \begin{bmatrix} a & b \\ c & d \end{bmatrix} S=[acbd]
其中 a , b , c , d a, b, c, d a,b,c,d 都是 n × n n \times n n×n 的矩阵, S S S 是 2 n × 2 n 2n \times 2n 2n×2n 的矩阵,这样才能满足上面的递推关系。我们根据上面的关系又能得出:
[ F k E ] = [ F k − 1 E ] × [ a b c d ] → F k = F k − 1 × a + c , E = F k − 1 × b + d \begin{aligned} & \begin{bmatrix} F_k & E \end{bmatrix} = \begin{bmatrix} F_{k - 1} & E \end{bmatrix} \times \begin{bmatrix} a & b \\ c & d \end{bmatrix} \\ \rightarrow & \quad F_k = F_{k - 1} \times a + c \; , \;E \; = F_{k - 1} \times b + d \end{aligned} →[FkE]=[Fk−1E]×[acbd]Fk=Fk−1×a+c,E=Fk−1×b+d
又因为我们有:
F k = F k − 1 + A k F_k = F_{k - 1} + A^k Fk=Fk−1+Ak
就能得到:
F k − 1 + A k = F k − 1 × a + c , E = F k − 1 × b + d F_{k - 1} + A^k = F_{k - 1} \times a + c \; , \; E = F_{k - 1} \times b + d Fk−1+Ak=Fk−1×a+c,E=Fk−1×b+d
带入 k = 1 k = 1 k=1,得到:
F 0 + A = F 0 × a + c , E = F 0 × b + d F_0 + A = F_0 \times a +c \; , \; E = F_0 \times b + d F0+A=F0×a+c,E=F0×b+d
带入 k = 2 k = 2 k=2 得到:
F 1 + A 2 = F 1 × a + c , E = F 1 × b + d F_1 + A^2 = F_1 \times a + c \; , \; E = F_1 \times b + d F1+A2=F1×a+c,E=F1×b+d
然后就能得到:
a = A , b = O , c = A , d = E a = A \; , \; b = O \; , \; c = A \; , \; d = E a=A,b=O,c=A,d=E
也就是:
S = [ A O A E ] , S k = [ A k O ∑ i = 1 k A i E ] S = \begin{bmatrix} A & O \\ A & E \end{bmatrix}, S^k = \begin{bmatrix} A^k & O \\ \sum\limits_{i = 1}^k A^i & E \end{bmatrix} S=[AAOE],Sk=⎣⎡Aki=1∑kAiOE⎦⎤
也就是说 S k S^k Sk 的左下角就是答案。
球形空间产生器
(2022.4.30)
给一个 n ∈ [ 1 , 10 ] n \in [1, 10] n∈[1,10] 维空间中的 n + 1 n + 1 n+1 个点的坐标,求出过这 n + 1 n +1 n+1 个点的 n n n 维空间上的 “球” 的球心的坐标。
我们发现对于球心的坐标: ( x 1 , x 2 , x 3 , ⋯ , x n ) (x_1, x_2, x_3, \cdots , x_n) (x1,x2,x3,⋯,xn),我们都有:
∑ j = 0 n ( a i , j − x i ) 2 = C o n s t \sum_{j = 0}^n (a_{i, j} - x_i)^2 = Const j=0∑n(ai,j−xi)2=Const
其中 a i , j a_{i, j} ai,j 第 i i i 个点的第 j j j 维坐标。我们发现这个方程是一个 n n n 元 2 2 2 次方程组,总共有 n + 1 n + 1 n+1 个方程,不是一个线性方程,但是我们可以通过一些变换将它变成线性方程然后通过高斯消元来求解。具体方法我们将相邻的两个方程作差:
∑ j = 1 n ( a i , j 2 − a i + 1 , j 2 − 2 x i ( a i , j − a i + 1 , j ) ) = 0 , ( i ∈ [ 1 , n ] , i ∈ Z ) → ∑ j = 1 n 2 ( a i , j − a i + 1 , j ) x j = ∑ j = 1 n ( a i , j 2 − a i + 1 , j 2 ) , ( i ∈ [ 1 , n ] , i ∈ Z ) \begin{aligned} &\sum_{j = 1}^n (a_{i, j}^2 - a_{i +1, j}^2 - 2x_i(a_{i, j} - a_{i +1, j})) = 0, \quad (i \in [1, n], i \in Z) \\ \rightarrow &\sum_{j = 1}^n2(a_{i, j} - a_{i +1, j})x_j = \sum_{j = 1}^n(a_{i, j}^2 - a_{i +1, j}^2), \quad (i \in [1, n], i \in Z) \end{aligned} →j=1∑n(ai,j2−ai+1,j2−2xi(ai,j−ai+1,j))=0,(i∈[1,n],i∈Z)j=1∑n2(ai,j−ai+1,j)xj=j=1∑n(ai,j2−ai+1,j2),(i∈[1,n],i∈Z)
现在我们发现这就是一个现线性的方程组了,也就是我们高斯消元这个增广矩阵:
[ 2 ( a 1 , 1 − a 2 , 1 ) 2 ( a 1 , 2 − a 2 , 2 ) ⋯ 2 ( a 1 , n − a 2 , n ) ∑ j = 1 n ( a 1 , j 2 − a 2 , j 2 ) 2 ( a 2 , 1 − a 3 , 1 ) 2 ( a 2 , 2 − a 3 , 2 ) ⋯ 2 ( a 2 , n − a 3 , n ) ∑ j = 1 n ( a 2 , j 2 − a 3 , j 2 ) ⋮ ⋮ ⋱ ⋮ ⋮ 2 ( a n , 1 − a n + 1 , 1 ) 2 ( a n , 2 − a n + 1 , 2 ) ⋯ 2 ( a n , n − a n + 1 , n ) ∑ j = 1 n ( a n , j 2 − a n + 1 , j 2 ) ] \begin{bmatrix} 2(a_{1, 1} - a_{2, 1}) & 2(a_{1, 2} - a_{2, 2}) & \cdots & 2(a_{1, n} - a_{2, n}) & \sum\limits_{j = 1}^n (a_{1, j}^2 - a_{2, j}^2) \\ 2(a_{2, 1} - a_{3, 1}) & 2(a_{2, 2} - a_{3, 2}) & \cdots & 2(a_{2, n} - a_{3, n}) & \sum\limits_{j = 1}^n (a_{2, j}^2 - a_{3, j}^2) \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 2(a_{n, 1} - a_{n +1, 1}) & 2(a_{n, 2} - a_{n +1, 2}) & \cdots & 2(a_{n, n} - a_{n +1, n}) & \sum\limits_{j = 1}^n (a_{n, j}^2 - a_{n +1, j}^2) \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡2(a1,1−a2,1)2(a2,1−a3,1)⋮2(an,1−an+1,1)2(a1,2−a2,2)2(a2,2−a3,2)⋮2(an,2−an+1,2)⋯⋯⋱⋯2(a1,n−a2,n)2(a2,n−a3,n)⋮2(an,n−an+1,n)j=1∑n(a1,j2−a2,j2)j=1∑n(a2,j2−a3,j2)⋮j=1∑n(an,j2−an+1,j2)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
就能得到答案了。