数论学习笔记
BJqxszx_zhuyukun
·
·
算法·理论
前置知识
初等数论，简单数学知识。
群论
很多数论知识都需要群论证明，我们先简单了解一下群论。
设 G 为一个非空集合，定义运算 \cdot，(G,\cdot) 为一个群当且仅当：
- 封闭性：\forall a,b \in G,a \cdot b \in G。
- 结合律：\forall a,b,c \in G,(a \cdot b) \cdot c=a \cdot (b \cdot c)。
- 单位元：\exist e \in G,\forall a \in G,e \cdot a=a \cdot e=a。
- 逆元：\forall a \in G,\exist b \in G,a \cdot b=b \cdot a=e，记作 b=a^{-1}。
如果一个群 G 的元素数量是有限的，那么它就是有限群，它的元素数量叫做它的阶，记作 \lvert G\rvert。
对于有限群 G，如果 H \subseteq G，并且 G,H 的运算相同，则称 H 为 G 的子群。
对于 a \in G，使 a^k=e 成立的最小正整数 k 称作 a 的阶，记作 \operatorname{ord}(a)。
拉格朗日定理
对于有限群 G，若 H 为 G 的子群，则 \lvert H\rvert \mid \lvert G\rvert，即 H 的阶整除 G 的阶。
以下证明中的 \cdot 运算符号省略。
$$
aH=\left\{ah \mid h\in H\right\}
$$
$aH$ 的元素与 $H$ 一一对应，因此 $\lvert aH\rvert=\lvert H\rvert$。
取两个不同的陪集 $aH,bH$，假设它们有公共元素，即 $\exist h_1,h_2 \in H,ah_1=bh_2$，则 $a=b h_2h_1^{-1}$。
又因为 $h_2h_1^{-1}\in H$，所以 $a\in bH$。
$\forall ah \in aH,ah=(bh_2h_1^{-1})h=b(h_2h_1^{-1}h)\in bH$，所以 $aH \subseteq bH$。
同理可证：$bH\subseteq aH$，因此 $aH = bH$。与假设“不同的陪集”矛盾，因此 $aH,bH$ 不存在公共元素。
设不同的陪集有 $k$ 个，可得：$\lvert G\rvert=k\lvert H\rvert$。即 $\lvert H\rvert \mid \lvert G\rvert$。
推论：设 \lvert G \rvert=m，\forall a \in G,a^m=e。
取 a 生成的循环子群 \langle a\rangle=\left\{a^x \mid x\in \mathbb{N}\right\}，容易知道其阶等于 a 的阶，且为 G 的子群。由拉格朗日定理得：\operatorname{ord}(a) \mid m。因此：
a^m=(a^{\operatorname{ord}(a)})^{\frac{m}{\operatorname{ord}(a)}}=e^{\frac{m}{\operatorname{ord}(a)}}=e
费马小定理
若 p 为素数，且 p \nmid a，则：
a^{p-1}\equiv1\pmod p
令 G=(\mathbb{Z}/p\mathbb{Z})^\times（即所有与 p 互质的剩余类），易知 \lvert G\rvert=p-1。
因为 p \nmid a，所以 a \in G，根据拉格朗日定理的推论，a^{\lvert G\rvert}=e，即 a^{p-1}\equiv 1\pmod p。
欧拉定理
若 \gcd(a,n)=1，则：
a^{\varphi(n)}\equiv 1\pmod n
令 G=(\mathbb{Z}/n\mathbb{Z})^\times，易知 \lvert G\rvert=\varphi(n)。
因为 \gcd(a,n)=1，所以 a \in G，a^{\lvert G\rvert}=e，即 a^{\varphi(n)}\equiv 1\pmod n。
扩展欧拉定理
a^b \equiv \begin{cases}
a^{b\ \bmod\ \varphi(p)} & \gcd(a,p)=1 \\
a^b & \gcd(a,p) \neq 1,b<\varphi(p) \\
a^{(b\ \bmod\ \varphi(p))+\varphi(p)} & \gcd(a,p) \neq 1,b \geq \varphi(p)
\end{cases}
\pmod p
> 第一项就是欧拉定理，当 $\gcd(a,p) \neq 1$ 的时候，把 $p$ 质因数分解，对于每一个 $p^k$ 证明原式成立，然后用 CRT 合并（CRT 见后文）即可。懒得证了。
# 裴蜀定理
若 $a,b,c$ 均为整数，$ax+by=c$ 存在整数解 $(x,y)$ 当且仅当 $\gcd(a,b) \mid c$。
> 设 $S=\left\{ax+by \mid x,y \in \mathbb{Z},ax+by>0\right\}$，易知 $S$ 为 $\mathbb{N^+}$ 的非空子集。设 $S$ 中元素的最小值为 $d$。容易知道必然存在整数 $x,y$ 使得 $d=ax+by$。
>
> 可以知道存在整数 $q,r$ 满足 $a=qd+r$，其中 $0 \leq r<d$。
>
> 把 $d$ 带入：$r=a-qd=a-q(ax+by)=a(1-qx)+b(-qy)$。
>
> 如果 $r>0$，那么 $r \in S$ 并且 $r<d$，矛盾，因此 $r=0$，即 $d \mid a$。
>
> 同理，可证：$d \mid b$。因此 $d$ 为 $a,b$ 的公因数。
>
> 设 $k$ 是 $a,b$ 的一个公因数，可得 $k \mid d$，即 $k \leq d$，因此 $d$ 为 $a,b$ 的最大公因数，即 $d=\gcd(a,b)$。
>
> 设 $m=\frac{c}{\gcd(a,b)}$，显然把 $x,y$ 均乘上 $m$ 就可以得到 $ax+by=c$ 的解，所以 $\gcd(a,b) \mid c$ 的充分性得证；
>
> 因为 $d=\gcd(a,b)$，所以 $d \mid (ax+by)$，即 $d \mid c$，$\gcd(a,b) \mid c$ 的必要性得证。
>
> 所以，$\gcd(a,b)\mid c$ 是存在整数解的充要条件。
# 扩展欧几里得算法（exGCD）
求 $ax+by=\gcd(a,b)$ 的一组特解。
由裴蜀定理容易知道一定有解。
如果 $b=0$，$\gcd(a,b)=a$，此时 $x=1,y=0$。
考虑 $a$ 除以 $b$，设 $a=qb+r$，其中 $0 \leq r<b$。可以推出 $\gcd(a,b)=\gcd(b,a\bmod b)=\gcd(b,r)$。
设 $\gcd(b,r)$ 这一步返回了 $x_0,y_0$，得到：
$$
bx_0+ry_0=\gcd(a,b) \\
bx_0+(a-qb)y_0=\gcd(a,b) \\
ay_0+b(x_0-qy_0)=\gcd(a,b)
$$
其中 $q=\lfloor \frac{a}{b} \rfloor$。
显然我们可以令：
$$
x_1=y_0,y_1=x_0-\lfloor \frac{a}{b}\rfloor y_0
$$
显然第一步是正确的，后面每一步是由前一步推导出的，所以也都是正确的。
这就是扩展欧几里得算法的流程。
时间复杂度：$O(\log (a+b))$。
# 扩展中国剩余定理（exCRT）
普通 CRT 感觉不如 exCRT，不讲了。
给定 $a_i,m_i$，解同余方程：
$$
\begin{cases}
x \equiv a_1 \pmod{m_1} \\
\cdots \\
x \equiv a_n \pmod{m_n}
\end{cases}
$$
$m_i$ 不一定两两互质。
考虑两两合并，对于 $x \equiv a_1\pmod{m_1}$ 与 $x\equiv a_2\pmod{m_2}$。
这等价于：$x=a_1+y_1m_1$ 与 $x=a_2+y_2m_2$。带入：$a_1+y_1m_1=a_2+y_2m_2$，即 $m_1y_1-m_2y_2=a_2-a_1$，利用扩展欧几里得算法算出一组 $(y_1,y_2)$，算不出来说明整个同余方程无解，算出来了就可以得到 $x\equiv a_1+y_1m_1\pmod{\operatorname{lcm}(m_1,m_2)}$。
时间复杂度：$O(n \log M)$。
# 大步小步算法（BSGS）
解方程：$a^x \equiv b \pmod p$，$a,p$ 互质。
设 $m=\lceil \sqrt{p}\rceil$，得到 $x=im+j$，其中 $0 \leq i<m,0 \leq j<m$。
带入：
$$
a^{im+j}=b \\
a^j=b(a^{-m})^i
$$
枚举 $j$，计算 $g^j \bmod p$ 记下来，然后枚举 $i$，计算 $t=b(a^{-m})^i \pmod p$，查找 $t$，找到了就说明此时的 $x$ 是解。如果最终都没找到，说明无解。
时间复杂度：$O(\sqrt{p})$。
# 扩展大步小步算法（exBSGS）
解方程：$a^x \equiv b\pmod{p}$，$a,p$ 不一定互质。
首先把 $a,b$ 取模，特判 $p=1$ 或 $b=1$ 的情况。
我们计算 $d=\gcd(a,p)$，如果 $b \bmod d \neq 0$，则无解。否则我们可以得到：$\frac{a}{d}a^{x-1}\equiv\frac{b}{d}\pmod{\frac{p}{d}}$。
令：
$$
b \leftarrow \frac{b}{d} \\
p \leftarrow\frac{p}{d}
$$
不断进行上面的操作，直到 $a,p$ 互质。设我们进行了 $k$ 次操作，$m$ 表示每一轮的 $\frac{a}{d}$ 的乘积，则我们得到了：
$$
ma^{x-k}\equiv b\pmod p \\
a^{x-k}\equiv bm^{-1}\pmod p
$$
此时 $a,p$ 互质，跑 BSGS 即可。
但是还有一种可能，就是答案 $x<k$，所以在每次计算 $d=\gcd(a,p)$ 之前枚举 $x=0,1,\dots,k-1$ 就好。
时间复杂度：$O(\sqrt{p})$。
# Lucas 定理
设 $p$ 为质数，则：
$$
\binom{n}{m}\equiv\binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{m}{p}\rfloor}\binom{n\bmod p}{m\bmod p}\pmod p
$$
> 考虑组合意义，设 $n=qp+r,m=sp+t$，$0\leq r<p,0 \leq t<p$，将 $n$ 个球排成一圈，每 $q$ 个分成一组，每组 $p$ 个球，最终剩下 $r$ 个球没有组。设这 $q$ 组为 $G_1,G_2,\dots,G_q$，没有组的为 $G_0$。
>
> 考虑从 $n$ 个球中选出 $m$ 个，设在 $G_i$ 中选了 $k_i$ 个球（$1 \leq i \leq q,0 \leq k_i \leq p,0 \leq k_0 \leq r$），则 $\sum k_i+k_0=m$，总方案数为：
>
> $$
> \binom{n}{m}\equiv\binom{r}{k_0}\prod_{i=1}^q\binom{p}{k_i}\pmod p
> $$
>
> 假设 $\exist 1 \leq i \leq q,1 \leq k_i \leq p-1$，那么有 $\binom{p}{k_i}\equiv0\pmod p$，因此在总方案数中，它们没有贡献，我们只需考虑 $k_i=0$ 或 $k_i=p$ 的部分产生的贡献。又因为 $m=sp+t$，所以一共有 $s$ 个 $k_i=p$，而 $k_0=t$。所以总方案数为：
> $$
> \binom{q}{s}\binom{r}{t}
> $$
> 因此：
> $$
> \binom{n}{m}\equiv\binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{m}{p}\rfloor}\binom{n\bmod p}{m\bmod p}\pmod p
> $$
# 扩展 Lucas 定理（exLucas）
计算
$$
\binom{n}{m} \bmod p
$$
其中 $p$ 不一定是质数。
我们考虑把 $p$ 质因数分解，对于每一个 $p^k$ 分别求解，最后用 CRT 或 exCRT 合并。
设 $f(n)$ 表示在 $n!$ 中，含有多少个质因子 $p$。容易知道：
$$
f(n)=\sum_{i=1}^\infin \lfloor\frac{n}{p^i}\rfloor=\lfloor\frac{n}{p}\rfloor+f(\lfloor \frac{n}{p}\rfloor)
$$
那么 $\binom{n}{m}$ 中包含的质因子 $p$ 的数量 $c=f(n)-f(m)-f(n-m)$。
设 $g(n)$ 表示在 $n!$ 中，除以所有的 $p$ 之后剩下的数对 $p^k$ 取模结果。即：
$$
g(n)=\frac{n!}{p^{f(n)}} \bmod p^k
$$
推一些式子：
$$
n!=\left(\prod_{1 \leq i \leq n,p \nmid i} i\right)\left(\prod_{i=1}^{\lfloor\frac{n}{p}\rfloor}i\cdot p\right) \\
=\left(\prod_{1\leq i\leq n,p \nmid i}i\right)p^{\lfloor\frac{n}{p}\rfloor}\left(\lfloor\frac{n}{p}\rfloor!\right)
$$
两边同时除以 $p^{f(n)}$：
$$
g(n)=\left(\prod_{1\leq i\leq n,p \nmid i}i\right)g\left(\lfloor\frac{n}{p}\rfloor\right)
$$
其中 $g(0)=1$。
考虑快速计算 $\prod_{1\leq i\leq n,p \nmid i}i$，我们知道，$x\equiv x+p^k \pmod{p^k}$，因此我们可以仅仅计算从 $1$ 到 $\frac{n}{p^k}$ 部分中，不被 $p$ 整除的数的乘积。对其求前缀积 $h(n)$ 可得：
$$
\left(\prod_{1\leq i\leq n,p \nmid i}i\right)\equiv h\left(p^k\right)^{\lfloor\frac{n}{p^k}\rfloor}h\left(n \bmod \left(p^k\right)\right) \pmod{p^k}
$$
综合起来：
$$
\binom{n}{m}\equiv p^cg(n)g(m)^{-1}g(n-m)^{-1} \pmod{p^k}
$$
然后 CRT 或 exCRT 合并，于是就做完了。时间复杂度：$O(p)$。
# 整除分块
求：
$$
\sum_{i=1}^n \lfloor\frac{n}{i}\rfloor
$$
性质：$\lfloor\frac{n}{i}\rfloor$ 取值只有 $O(\sqrt{n})$ 种。
> 当 $i \leq \sqrt{n}$，$i$ 的取值最多 $\sqrt{n}$ 个，因此原式取值最多 $\sqrt{n}$ 个。
>
> 当 $i>\sqrt{n}$，$\lfloor\frac{n}{i}\rfloor<\sqrt{n}$，取值最多 $\sqrt{n}$ 个。
显然每种取值也是连续的。所以我们考虑对于每一段计算。
设当前段左端点为 $l$，右端点为 $r$，那么这一段的值为 $v=\lfloor\frac{n}{l}\rfloor$，$r$ 为满足 $\lfloor\frac{n}{r}\rfloor=v$ 的最大整数。
有如下式子：
$$
v\leq\frac{n}{r}<v+1 \\
r \leq \lfloor\frac{n}{l}\rfloor
$$
因此 $r=\lfloor\frac{n}{\lfloor\frac{n}{l}\rfloor}\rfloor$。
时间复杂度：$O(\sqrt{n})$。
如果所求式子带权：
$$
\sum_{i=1}^n f(i)\lfloor\frac{n}{i}\rfloor
$$
处理 $f$ 的前缀和，每一段乘上 $f$ 的区间和即可。
如果是二维整除分块：
$$
\sum_{i=1}^{\min(n,m)}\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{i}\rfloor
$$
其他部分不变，只考虑 $r$ 怎么求。类似上面的式子推一下就能得到：$r=\min\left(\lfloor\frac{n}{\lfloor\frac{n}{l}\rfloor}\rfloor,\lfloor\frac{m}{\lfloor\frac{m}{l}\rfloor}\rfloor\right)$。
$k$ 维类似。时间复杂度 $O(k\sqrt{n})$。
# 莫比乌斯反演
$$
\mu(n)=\begin{cases}1&n=1\\(-1)^k&n=p_1p_2p_3\dots p_k,\forall 1 \leq i \leq k,p_i\in\mathbb{P},\forall 1 \leq i<j \leq k,p_i \neq p_j \\0 & \exist p \in \mathbb{P},p^2 \mid n\end{cases}
$$
意思就是，如果 $n=1$，那么 $\mu(n)=1$；如果 $n$ 的质因数分解中有 $k$ 个质因子，并且指数都为 $1$，那么 $\mu(n)=(-1)^k$；否则 $\mu(n)=0$。
性质：
$$
\sum_{d \mid n}\mu(d)=\begin{cases}1&n=1\\0&n\neq1\end{cases}
$$
> 首先第一项是显然的，我们证明第二项。
>
> 容易发现形如 $\prod_{i=1}^k p_i$ 形式的 $d$ 才是有贡献的。设 $n$ 有 $k$ 个不同的质因子，我们从中选取若干个组成了 $d$。如果选取了偶数个，则 $\mu(d)=1$，如果选取了奇数个，则 $\mu(d)=-1$。
>
> 设这 $k$ 个不同的质因子中有一个特殊质因子 $p$。我们先在其他 $k-1$ 个质因子中选取，最后决定选不选 $p$。显然，最后选 $p$ 的方案数与不选 $p$ 的方案数相同，但是对于同一种选取其他质因子的方案，选不选 $p$ 之后的集合大小奇偶性不同。因此，选取偶数个的方案数等于选取奇数个的方案数。也就是 $\sum_{d \mid n}\mu(d)=0\ (n\neq1)$。
假设我们已知：
$$
f(n)=\sum_{n \mid d} g(d)
$$
那么根据反演有：
$$
g(n)=\sum_{n\mid d}\mu\left(\frac{d}{n}\right)f(d)
$$
> $$
> \sum_{n\mid d}\mu\left(\frac{d}{n}\right)f(d)=\sum_{n \mid d}\mu\left(\frac{d}{n}\right)\sum_{d \mid k}g(k)
> $$
>
> 令 $t=\frac{d}{n}$，有：
> $$
> \sum_{n \mid d}\mu\left(\frac{d}{n}\right)\sum_{d \mid k}g(k)=\sum_{n \mid k}g(k)\sum_{t \mid \frac{k}{n}} \mu(t)=\sum_{n\mid k}g(k)\left[\frac{k}{n}=1\right]=g(n)
> $$
> 得证。
$\mu(n)$ 一般使用线性筛求出，代码如下：
```cpp
mu[1]=1;
for(int i=2;i<=n;i++){
if(!isp[i]){//是质数
prm[++pcnt]=i;//记录质数
mu[i]=1;
}
for(int j=1;j<=pcnt&&i*prm[j]<=n;j++){
isp[i*prm[j]]=true;//标记合数
if(i%prm[j]==0){//prm[j]是平方因子
mu[i*prm[j]]=0;
break;
}
else mu[i*prm[j]]=-mu[i];
}
}
```
# 例题
### [P5091 【模板】扩展欧拉定理](https://www.luogu.com.cn/problem/P5091)
> 求 $a^b \bmod m$，$1 \leq a\leq 10^9,1 \leq b \leq 10^{20000000},1 \leq m \leq 10^8$。
模板不讲了。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
ll a,m,mod,x,p;
bool flag;
string b;
ll qpow(ll a,ll b){
ll res=1%mod;
while(b){
if(b&1) res=res*a%mod;
a=a*a%mod;
b>>=1;
}
return res;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
cin>>a>>m>>b;
a%=m;
p=mod=m;
for(ll i=2;i*i<=m;i++){
if(m%i==0){
p=p/i*(i-1);
while(m%i==0) m/=i;
}
}
if(m>1) p=p/m*(m-1);
for(auto &i:b){
x=x*10+i-'0';
if(x>=p){
flag=true;
x%=p;
}
}
if(flag) x+=p;
cout<<qpow(a,x);
return 0;
}
```
### [P5656 【模板】二元一次不定方程 (exgcd)](https://www.luogu.com.cn/problem/P5656)
> 给定不定方程 $ax+by=c$，无解输出 $-1$；有正整数解输出正整数解数量，正整数解中 $x$ 的最小值、最大值，正整数解中 $y$ 的最小值、最大值；无正整数解有整数解输出整数解中 $x$ 的最小正整数值，$y$ 的最小正整数值。
判无解用裴蜀定理，求特解 $x_0,y_0$ 使用 exgcd，考虑推出通解形式。有：
$$
a(x_0+db)+b(y_0-da)=c
$$
其中 $db,da$ 为整数。所以，最小的 $d=\frac{1}{\gcd(a,b)}$。令 $d_1=\frac{b}{\gcd(a,b)},d_2=\frac{a}{\gcd(a,b)}$，$s$ 为整数，有通解：
$$
x=x_0+sd_1, y=y_0-sd_2
$$
当 $x$ 为正整数时，有 $x>0$，即：
$$
x_0+sd_1>0 \rightarrow s>-\frac{x_0}{d_1}
$$
当 $y$ 为正整数时，有 $y>0$，即：
$$
y_0-sd_2>0 \rightarrow s<\frac{y_0}{d_2}
$$
分类讨论即可。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define getx(n) (x+(n)*d1)
#define gety(n) (y-(n)*d2)
LL _,a,b,c,x,y,g,d1,d2;
void read(LL &x){
x=0;
char ch=getchar();
int f=1;
while((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
if(ch=='-'){
f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9'){
x=x*10+(ch^48);
ch=getchar();
}
x*=f;
return;
}
void write(const LL &x){
if(x<0){
putchar('-');
write(-x);
return;
}
if(x<=9){
putchar(x^48);
return;
}
write(x/10);
putchar(x%10^48);
return;
}
LL div1(const LL &a,const LL &b){
if(a%b&&(a<0&&b>0||a>0&&b<0)) return a/b-1;
return a/b;
}
LL div2(const LL &a,const LL &b){return -div1(-a,b);}
LL exgcd(const LL &a,const LL &b,LL &x,LL &y){
if(b==0){
x=1,y=0;
return a;
}
LL g=exgcd(b,a%b,y,x);
y-=a/b*x;
return g;
}
void solve(){
read(a);
read(b);
read(c);
g=exgcd(a,b,x,y);
if(c%g){
printf("-1\n");
return;
}
x*=c/g,y*=c/g;
d1=b/g,d2=a/g;
if(div2(y,d2)-1-div1(-x,d1)<=0){
write(getx(div1(-x,d1)+1));
putchar(' ');
write(gety(div2(y,d2)-1));
putchar('\n');;
}
else{
write(div2(y,d2)-1-div1(-x,d1));
putchar(' ');
write(getx(div1(-x,d1)+1));
putchar(' ');
write(gety(div2(y,d2)-1));
putchar(' ');
write(getx(div2(y,d2)-1));
putchar(' ');
write(gety(div1(-x,d1)+1));
putchar('\n');
}
return;
}
int main(){
read(_);
while(_--) solve();
return 0;
}
```
### [P4774 [NOI2018] 屠龙勇士](https://www.luogu.com.cn/problem/P4774)
> 有 $n$ 条龙，$m$ 把剑，你需要按顺序屠龙。每次屠龙，你需要选择攻击力不高于龙血量 $a_i$ 的一把攻击力最高的剑，如果没有选取攻击力最低的，并攻击龙 $x$ 次。设本次攻击力为 $b_i$，则你需要满足 $b_ix \equiv a_i\pmod{p_i}$，并且 $b_ix\geq a_i$。击败这只龙之后你会获得一把新的剑，旧的剑消失。求最小的正整数 $x$ 或者报告无解。
我们先不考虑这个同余式和不等式，思考应该怎么选剑。显然，我们可以用一个 multiset 和二分查找完成。这样我们就确定了每一个 $b_i$。然后我们跑一个 exCRT 求出最小的 $x$，再利用通解公式求出满足不等式的最小解即可。
简单说一下怎么求通解。设 $P=\operatorname{lcm}_{i=1}^n(p_i)$，我们求出来的解为 $x_0$，则有：
$$
x_0\equiv C\pmod P \\
x\equiv C\pmod P
$$
求出 $C$，可以得到：$x=sP+C$，其中 $s$ 为整数。假设不等式的解集为 $x \geq T$，则 $sP+C \geq T$，$s \geq \lceil\frac{T-C}{P}\rceil$，答案即为 $\lceil\frac{T-C}{P}\rceil P+C$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll MAXN=1e5+15;
ll _,n,m,atk,lcm,low,m1,m2,g,c,x,y,ans;
ll a[MAXN],p[MAXN],q[MAXN],b[MAXN];
multiset<ll> mset;
void exgcd(const ll &a,const ll &b,ll &x,ll &y){
if(b==0){
x=1,y=0;
return;
}
exgcd(b,a%b,y,x);
y-=a/b*x;
return;
}
ll inv(const ll &a,const ll &mod){
ll x,y;
exgcd(a,mod,x,y);
return (x%mod+mod)%mod;
}
void solve(){
mset.clear();
cin>>n>>m;
for(ll i=1;i<=n;i++) cin>>a[i];
for(ll i=1;i<=n;i++) cin>>p[i];
for(ll i=1;i<=n;i++) cin>>q[i];
for(ll i=1;i<=m;i++){
cin>>atk;
mset.insert(atk);
}
for(ll i=1;i<=n;i++){
auto it=mset.upper_bound(a[i]);
if(it!=mset.begin()) it--;
b[i]=*it;
mset.erase(it);
mset.insert(q[i]);
}
lcm=1,low=ans=0;
for(ll i=1;i<=n;i++){
low=max((a[i]+b[i]-1)/b[i],low);
g=__gcd(b[i],p[i]);
if(a[i]%g){
cout<<"-1\n";
return;
}
a[i]/=g,b[i]/=g,p[i]/=g;
a[i]=(__int128)a[i]%p[i]*inv(b[i]%p[i],p[i])%p[i];
m1=lcm,m2=p[i];
g=__gcd(m1,m2);
c=a[i]-ans;
if(c%g){
cout<<"-1\n";
return;
}
m1/=g,m2/=g,c/=g;
exgcd(m1,m2,x,y);
lcm*=m2;
x=((__int128)x*c%m2+m2)%m2;
ans=((__int128)x*g%lcm*m1%lcm+ans)%lcm;
}
if(ans<low){
low=(low-ans+lcm-1)/lcm;
ans+=low*lcm;
}
cout<<ans<<'\n';
return;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
cin>>_;
while(_--) solve();
return 0;
}
```
### [P4345 [SHOI2015] 超能粒子炮·改](https://www.luogu.com.cn/problem/P4345)
> 求：
> $$
> \sum_{i=0}^k \binom{n}{i} \bmod 2333
> $$
设 $p=2333$，再设函数 $f$：
$$
f(n,k)=\sum_{i=0}^k \binom{n}{i} \bmod p \\
=\sum_{i=0}^k \binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{i}{p}\rfloor}\binom{n \bmod p}{i \bmod p}\bmod p \\
=\left(\sum_{i=0}^{\left\lfloor\frac{k}{p}\right\rfloor-1}\binom{\lfloor\frac{n}{p}\rfloor}{i}\left(\sum_{j=0}^{p-1}\binom{n \bmod p}{j}\right)+\binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{k}{p}\rfloor}\left(\sum_{j=0}^{k\ \bmod\ p}\binom{n \bmod p}{j}\right)\right) \bmod p \\
=\left(2^{n\ \bmod\ p}\sum_{i=0}^{\left\lfloor\frac{k}{p}\right\rfloor-1}\binom{\lfloor\frac{n}{p}\rfloor}{i}+\binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{k}{p}\rfloor}\left(\sum_{j=0}^{k\ \bmod\ p}\binom{n \bmod p}{j}\right)\right) \bmod p \\
=\left(2^{n\ \bmod\ p}f\left(\left\lfloor\frac{n}{p}\right\rfloor,\left\lfloor\frac{k}{p}\right\rfloor-1\right)+\binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{k}{p}\rfloor}\left(\sum_{j=0}^{k\ \bmod\ p}\binom{n \bmod p}{j}\right)\right) \bmod p
$$
加号后面的可以预处理前缀和，递归计算即可。
时间复杂度：$O(p^2+T\log_p^2n)$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll MAXN=3005,MOD=2333;
ll _,n,k;
ll fact[MAXN],invfac[MAXN],pows[MAXN],pre[MAXN][MAXN];
ll qpow(ll a,ll b){
ll res=1;
while(b){
if(b&1) res=res*a%MOD;
a=a*a%MOD;
b>>=1;
}
return res;
}
ll C(const ll &n,const ll &m){
if(m<0||n<m) return 0;
return fact[n]*invfac[m]%MOD*invfac[n-m]%MOD;
}
void init(){
fact[0]=1;
for(ll i=1;i<MOD;i++) fact[i]=i*fact[i-1]%MOD;
invfac[MOD-1]=qpow(fact[MOD-1],MOD-2);
for(ll i=MOD-1;i;i--) invfac[i-1]=i*invfac[i]%MOD;
pows[0]=1;
for(ll i=1;i<MOD;i++) pows[i]=pows[i-1]*2%MOD;
for(int i=0;i<MOD;i++){
pre[i][0]=1;
for(int j=1;j<MOD;j++) pre[i][j]=(pre[i][j-1]+C(i,j))%MOD;
}
return;
}
ll lucas(const ll &n,const ll &m){
if(m<0||n<m) return 0;
if(n<MOD) return C(n,m);
return lucas(n/MOD,m/MOD)*lucas(n%MOD,m%MOD)%MOD;
}
ll f(const ll &n,const ll &k){
ll res=0;
if(n<MOD){
for(ll i=0;i<=min(k,MOD-1);i++) res=(res+C(n,i))%MOD;
return res;
}
res=pows[n%MOD]*f(n/MOD,k/MOD-1);
res=(res+lucas(n/MOD,k/MOD)*pre[n%MOD][k%MOD]%MOD)%MOD;
return res;
}
void solve(){
cin>>n>>k;
k=min(n,k);
cout<<f(n,k)<<'\n';
return;
}
int main(){
cin>>_;
init();
while(_--) solve();
return 0;
}
```
### [P3773 [CTSC2017] 吉夫特](https://www.luogu.com.cn/problem/P3773)
> 给定长度为 $n$ 的序列 $a$，求有多少个长度至少为 $2$ 的不升子序列 $a'$ 满足：
> $$
> \prod_{i=2}^k \binom{a'_{i-1}}{a'_{i}}\equiv1\pmod 2
> $$
> 求出结果对 $10^9+7$ 取模。
首先根据 Lucas 定理有：
$$
\prod_{i=2}^k \binom{a'_{i-1}}{a'_{i}}\bmod 2 \\
=\prod_{i=2}^k\binom{\left\lfloor\frac{a'_{i-1}}{2}\right\rfloor}{\left\lfloor\frac{a'_{i}}{2}\right\rfloor}\binom{a'_{i-1}\bmod2}{a'_i\bmod2}\bmod2>0
$$
说明 $\binom{a'_{i-1}\ \bmod\ 2}{a'_i\ \bmod\ 2}=1$。
继续展开左边的这一坨分式，我们发现整体式子非 $0$ 当且仅当 $a'_{i-1}$ 与 $a'_i$ 在二进制每一位上 $a'_{i-1,j}\geq a'_{i,j}$ 成立。也就是定义 $\&$ 表示按位与，$a'_{i-1}\&a'_{i}=a'_i$。这也可以看作 $a'_i$ 是 $a'_{i-1}$ 的子集。
注意到 $a_i$ 互不相同，设 $dp_i$ 表示结尾是 $i$ 的序列数量，枚举子集 $j$ 刷表法转移即可。
时间复杂度：$O(3^{\log_2\max a_i})$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll MAXN=1e6+16,MOD=1e9+7;
ll n,ans;
ll a[MAXN],dp[MAXN];
int main(){
cin>>n;
for(ll i=1;i<=n;i++){
cin>>a[i];
for(ll j=a[i]&a[i]-1;j;j=j-1&a[i]) dp[j]=(dp[j]+dp[a[i]]+1)%MOD;
ans=(ans+dp[a[i]])%MOD;
}
cout<<ans;
return 0;
}
```
### [P2260 [清华集训 2012] 模积和](https://www.luogu.com.cn/problem/P2260)
> 求：
> $$
> \sum_{i=1}^n\sum_{j=1}^m(n\bmod i)(m\bmod j),i\neq j
> $$
规定 $n\leq m$，容易想到把取模变成整除形式：
$$
\sum_{i=1}^n\sum_{j=1}^m(n\bmod i)(m\bmod j),i\neq j \\
=\sum_{i=1}^n\sum_{j=1}^m\left(n-i\left\lfloor\frac{n}{i}\right\rfloor\right)\left(m-j\left\lfloor\frac{m}{j}\right\rfloor\right)-\sum_{i=1}^{n}\left(n-i\left\lfloor\frac{n}{i}\right\rfloor\right)\left(m-i\left\lfloor\frac{m}{i}\right\rfloor\right) \\
=\sum_{i=1}^n\sum_{j=1}^m nm-mi\left\lfloor\frac{n}{i}\right\rfloor-nj\left\lfloor\frac{m}{j}\right\rfloor+ij\left\lfloor\frac{n}{i}\right\rfloor \left\lfloor\frac{m}{j}\right\rfloor-\sum_{i=1}^n nm-mi\left\lfloor\frac{n}{i}\right\rfloor-ni\left\lfloor\frac{m}{i}\right\rfloor+i^2\left\lfloor\frac{n}{i}\right\rfloor \left\lfloor\frac{m}{i}\right\rfloor \\
=n^2m^2-n^2m-m^2\sum_{i=1}^n i\left\lfloor\frac{n}{i}\right\rfloor-n^2\sum_{i=1}^m i\left\lfloor\frac{m}{i}\right\rfloor+\left(\sum_{i=1}^n i\left\lfloor\frac{n}{i}\right\rfloor\right)\left(\sum_{i=1}^mi\left\lfloor\frac{m}{i}\right\rfloor\right)+m\sum_{i=1}^ni\left\lfloor\frac{n}{i}\right\rfloor+n\sum_{i=1}^ni\left\lfloor\frac{m}{i}\right\rfloor-\sum_{i=1}^ni^2\left\lfloor\frac{n}{i}\right\rfloor \left\lfloor\frac{m}{i}\right\rfloor
$$
我们知道：
$$
\sum_{i=1}^n i^2=\frac{n(n+1)(2n+1)}{6}
$$
因此整除分块即可。注意一些循环边界。
时间复杂度：$O(\sqrt{n})$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll MOD=19940417;
ll n,m,ans,x,y;
int main(){
cin>>n>>m;
if(n>m) swap(n,m);
ans=n*n%MOD*m%MOD*m%MOD;
ans-=n*n%MOD*m%MOD;
ans=(ans%MOD+MOD)%MOD;
for(ll i=1,j=1;j<=n;i=j+1){
if(n/i==0) break;
j=n/(n/i);
x=(x+n/i*((i+j)*(j-i+1)/2%MOD)%MOD)%MOD;
}
ans-=x*m%MOD*m%MOD;
ans=(ans%MOD+MOD)%MOD;
for(ll i=1,j=1;j<=m;i=j+1){
if(m/i==0) break;
j=m/(m/i);
y=(y+m/i*((i+j)*(j-i+1)/2%MOD)%MOD)%MOD;
}
ans-=y*n%MOD*n%MOD;
ans=(ans%MOD+MOD)%MOD;
ans=(ans+x*y%MOD)%MOD;
ans=(ans+x*m%MOD)%MOD;
x=0;
for(ll i=1,j=1;i<=n;i=j+1){
if(m/i==0) break;
j=min(n,m/(m/i));
x=(x+m/i*((i+j)*(j-i+1)/2%MOD)%MOD)%MOD;
}
ans=(ans+x*n%MOD)%MOD;
x=y=0;
for(ll i=1,j=1;j<=n;i=j+1){
if(n/i==0||m/i==0) break;
j=min(n/(n/i),m/(m/i));
x=(__int128)j*(j+1)*(j*2+1)/6%MOD;
x-=(__int128)(i-1)*i*((i-1)*2+1)/6%MOD;
x=(x%MOD+MOD)%MOD;
y=(y+(n/i)*(m/i)%MOD*x%MOD)%MOD;
}
ans=(ans-y+MOD)%MOD;
cout<<ans;
return 0;
}
```
### [P2522 [HAOI2011] Problem b](https://www.luogu.com.cn/problem/P2522)
> 给出 $n$ 次询问，每次给定 $a,b,c,d,k$，求：
> $$
> \sum_{i=a}^b\sum_{j=c}^d \left[\gcd(i,j)=k\right]
> $$
设：
$$
f(a,b,c,d)=\sum_{i=a}^b\sum_{j=c}^d \left[\gcd(i,j)=k\right]
$$
容斥则有：
$$
ans=f(1,b,1,d)-f(1,a-1,1,d)-f(1,b,1,c-1)+f(1,a-1,1,c-1)
$$
钦定 $n \leq m$，所以只需要求：
$$
T(n,m)=\sum_{i=1}^n\sum_{j=1}^m\left[\gcd(i,j)=k\right] \\
$$
设 $F(d)$ 表示最大公约数为 $d$ 的数对数量，$G(d)$ 表示满足 $d\mid \gcd(i,j)$ 的数对数量。得到：
$$
G(d)=\sum_{d\mid k}F(k) \\
G(d)=\left\lfloor\frac{n}{d}\right\rfloor\left\lfloor\frac{m}{d}\right\rfloor
$$
明显的莫比乌斯反演：
$$
F(k)=\sum_{k\mid d}\mu\left(\frac{d}{k}\right)G(d) \\
=\sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\mu(i)\left\lfloor\frac{n}{ki}\right\rfloor\left\lfloor\frac{m}{ki}\right\rfloor
$$
令 $N=\left\lfloor\frac{n}{k}\right\rfloor,M=\left\lfloor\frac{m}{k}\right\rfloor$：
$$
F(k)=\sum_{i=1}^{N}\mu(i)\left\lfloor\frac{N}{i}\right\rfloor\left\lfloor\frac{M}{i}\right\rfloor
$$
整除分块即可。
时间复杂度：$O(N+n\sqrt{N})$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll MAXN=5e4+54;
ll _,a,b,c,d,k,pcnt;
ll mu[MAXN],pre[MAXN],prm[MAXN];
bool isp[MAXN];
void init(){
mu[1]=1;
for(ll i=2;i<=50000;i++){
if(!isp[i]){
prm[++pcnt]=i;
mu[i]=-1;
}
for(ll j=1;j<=pcnt&&i*prm[j]<=50000;j++){
isp[i*prm[j]]=true;
if(i%prm[j]==0){
mu[i*prm[j]]=0;
break;
}
else mu[i*prm[j]]=-mu[i];
}
}
for(ll i=1;i<=50000;i++) pre[i]=mu[i]+pre[i-1];
return;
}
ll f(const ll &n,const ll &m){
ll tn=n/k,tm=m/k,res=0;
if(tn>tm) swap(tn,tm);
for(ll i=1,j=1;i<=tn;i=j+1){
if(tn/i==0||tm/i==0) break;
j=min(tn/(tn/i),tm/(tm/i));
res+=(pre[j]-pre[i-1])*(tn/i)*(tm/i);
}
return res;
}
void solve(){
cin>>a>>b>>c>>d>>k;
cout<<f(b,d)-f(a-1,d)-f(b,c-1)+f(a-1,c-1)<<'\n';
return;
}
int main(){
init();
cin>>_;
while(_--) solve();
return 0;
}
```
### [P2257 YY的GCD](https://www.luogu.com.cn/problem/P2257)
> 求：
> $$
> \sum_{i=1}^n\sum_{j=1}^m\operatorname{isp}(\gcd(i,j))
> $$
> 其中 $\operatorname{isp}(x)$ 表示 $x$ 是否为质数。
依旧令 $n\leq m$，推式子：
$$
\sum_{i=1}^n\sum_{j=1}^m\operatorname{isp}(\gcd(i,j))\\
=\sum_{p\in\mathbb{P}}\sum_{i=1}^n\sum_{j=1}^m\left[\gcd(i,j)=p\right] \\
$$
显然和上一题一模一样。令 $N=\left\lfloor\frac{n}{p}\right\rfloor,M=\left\lfloor\frac{m}{p}\right\rfloor$：
$$
ans=\sum_{p\in\mathbb{P}}\sum_{i=1}^{N}\mu(i)\left\lfloor\frac{N}{i}\right\rfloor\left\lfloor\frac{M}{i}\right\rfloor\\
=\sum_{i=1}^n\left\lfloor\frac{n}{i}\right\rfloor\left\lfloor\frac{m}{i}\right\rfloor\left(\sum_{p\mid i,p\in\mathbb{P}}\mu\left(\frac{i}{p}\right)\right)\\
$$
右边的一坨令其为 $f(i)$，预处理即可。
时间复杂度：$O(n\log\log n+T\sqrt{n})$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int MAXN=1e7+17;
int _,n,m,pcnt;
long long ans;
long long mu[MAXN],pre[MAXN],f[MAXN],prm[1000005];
bool isp[MAXN];
void init(){
mu[1]=1;
for(int i=2;i<MAXN;i++){
if(!isp[i]){
prm[++pcnt]=i;
mu[i]=-1;
}
for(int j=1;j<=pcnt&&i*prm[j]<MAXN;j++){
isp[i*prm[j]]=true;
if(i%prm[j]==0){
mu[i*prm[j]]=0;
break;
}
else mu[i*prm[j]]=-mu[i];
}
}
for(int i=1;i<=pcnt;i++) for(int j=prm[i];j<MAXN;j+=prm[i]) f[j]+=mu[j/prm[i]];
for(int i=1;i<MAXN;i++) pre[i]=f[i]+pre[i-1];
return;
}
void solve(){
ans=0;
cin>>n>>m;
if(n>m) swap(n,m);
for(int i=1,j=1;i<=n;i=j+1){
if(n/i==0||m/i==0) break;
j=min(n/(n/i),m/(m/i));
ans+=(pre[j]-pre[i-1])*(n/i)*(m/i);
}
cout<<ans<<'\n';
return;
}
int main(){
init();
cin>>_;
while(_--) solve();
return 0;
}
```