数论在算法方面的基础应用

数论基础

模运算&算数基本定理

模运算常用恒等式&符号

注：除法取模要使用逆元

代码实现

1
2
3
4
5
int a,b,c;
a %= c,b %= c;
int ans1 = (a + b) % c;
int ans2 = ((a - b) % c + c) % c; //三步走：模 加 模
int ans3 = ((long long) a * b) % c; //一定要记得开成longlong

算数基本定理

任何一个大于$1$的自然数$N$，如果$N$不是质数，那么它一定可以唯一分解成有限个质数，即

补充

  • 约数倍数的特征：如果一个数字$x=P_1^{b_1}P_2^{b_2}…P_n^{b_n}$是N的约数，那么一定有$\forall b_i\leq a_1$，同理如果x是N的倍数，那么一定有$\forall b_i\geq a_i$

  • 素因数求法：模拟短除法

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    inline prime_factor(int n){
        for(int i=2;i*i<=n;i++){
            if(n%i==0){
                ++cnt;
                prime_fac[cnt]=i;
            }
            while(n%i==0) n/=i;
        }
        if(n!=1){
            ++cnt;
            prime_fac[++cnt]=n;
        }
    }

关于快速分解质因数，见欧拉函数部分

同余

定义

如果对于两个数字$a,b$，对于一个数字$m$的余数相等，则这两个数对于$m$同余

在数学上的标准定义：设$m$施给定的一个正整数，$a,b$是整数，若满足$m|(a-b)$，则称$a$与$b$对模$m$同余，记作$a\equiv b(mod\ m)$

逆元

  • 对于$\frac{7}{3} \% 2$的结果怎么求呢？（注意：这里不是整除）

定义

乘法逆元：对于一个数字$a$，在膜$p$意义下如果存在一个数字$b$，满足$a*b\equiv 1(mod\ p)$，我们就称$b$是模$p$意义下$a$的乘法逆元，记作$a^{-1}$。容易发现

于是，我们就可以把除法取模转换成乘法取模啦！

接下来就是各种求逆元的方法

欧拉定理

前置知识：

欧拉函数（在下方素数筛法里面有讲）

定理：

如果$gcd(a,p)=1$，即$a,p$互质，则有$a^{\varphi (p)}\equiv 1(mod\ p)$

求逆元：

对于数字$a$，我所需要寻找的$a^{-1}$需要满足$a*a^{-1}\equiv 1(mod\ p)$，而$a^{\varphi (p)}\equiv 1(mod\ p)$，因此有$a* a^{\varphi (p)-1}\equiv 1(mod\ p)$，因此，$a^{-1}=a^{\varphi (p)-1}$，因此只要用一次快速幂求出逆元就可以了

费马小定理：

如果$p$是一个质数，则一定有$a^{p-1}\equiv 1(mod\ p)$

其他用法：

在$a,p$互质的前提下，计算$a^b(mod\ p)$，其中$b$可以大到$10^{10^7}$级别

方法：令$b=k*\varphi (p)+m$，原式转化为$a^b(mod\ p)=(a^{\varphi (p)})^k*a^m$，由欧拉定理得$a^{\varphi (p)}\equiv 1(mod\ p)$，因此原式转化为$a^b(mod\ p)=a^m$，而$m=b\% \varphi (p)$，所以最后答案就是$a^b\equiv a^{b\% \varphi (p)}(mod\ p)$

UPD 2020.6.12

扩展欧拉定理

洛谷P5091 【模板】扩展欧拉定理

题意：若$a,p$不互质，计算$a^b(mod\ p)$，其中$b$可以大到$10^{10^7}$级别

此时，欧拉定理就不适用了，我们需要找到一个更为广泛的解法，即 对欧拉定理进行扩展，先给出结论：

证明

image-20200612112612672

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
ll a,m,b,phim;
bool flag=false;

// a^b mod m = a^(b mod phi[m] + phi[m]) mod m {if b>=phi[m]}
template<class T>
inline void _read(T &x) {
	x=0;int f=0;char ch=getchar();
	while(ch<48||ch>57){f|=(ch=='-');ch=getchar();}
	while(ch>=48&&ch<=57){
		x=(x<<1)+(x<<3)+(ch^48);
		if(x>=phim) flag=true,x%=phim;// 边读入边取模，并标记是否有x>=phi[m]
		ch=getchar();
	}
	x=f?-x:x;
	return;
}

// get phi[n]
inline ll getphi(ll n){
	ll res=n;
	for(ll i=2;i*i<=n;i++){
		if(n%i==0) res=res/i*(i-1);
		while(n%i==0) n/=i;
	}
	if(n>1) res=res/n*(n-1);
	return res;
}

int main(){
	read(a),read(m);
	phim=getphi(m);// 求m欧拉函数
	_read(b);
	// get a^b mod m
	ll ans=qpow(a,b+(flag?phim:0),m);
	printf("%lld\n",ans);
	return 0;
}

扩展欧几里得算法

前置知识：

裴蜀定理：对于方程$ax+by=c$有整数解当且仅当$c$是$gcd(a,b)$的倍数，即$(gcd(a,b)|c)$

扩欧算法：

通常写作$exgcd$算法，它是用来求下面这个同余方程的一组解的

于是我们从另一个角度看逆元方程

因此，现在只需要找出一组$x,y$使得上式成立即可（为简单起见，可以把上式的减号换成加号），这个式子同样也印证了欧拉定理中如果想要存在$a$对$p$的逆元，必须要$gcd(a,p)=1$

如何找出一组解？

代码

1
2
3
4
5
6
7
8
9
10
11
12
ll x,y;
void exgcd(ll a,ll b){
    if(b==0){
        x=1,y=0;
        return;
    }
    exgcd(b,a%b);
    ll tx=x;
    ll ty=y;
    x=ty;
    y=tx-a/b*ty;
}

因此，求一个数$a$的逆元也可以转化为求解上述$x$的过程

那么，如果对于一个一般形式的二元一次不定方程（如下）

或者一个一般的同余方程（如下）

该如何求解呢？

代码如下（exgcd的过程中顺便求出gcd(a,b)）

1
2
3
4
5
6
7
8
9
10
11
12
13
ll x,y;
ll exgcd(ll a,ll b){
    if(b==0){
        x=1,y=0;
        return a;
    }
    ll t=exgcd(b,a%b);
    ll tx=x;
    ll ty=y;
    x=ty;
    y=tx-a/b*ty;
    return t
}//于是答案就是x*c/exgcd(a,b)

除此之外，如果我们又想要所有的解，或者$x$最小或$y$最小的正整数解，这就又涉及到另外一个问题，它所有解的通式要怎么写呢？

这里直接给出答案（证明太复杂，略）

如果再将其转化为一般式，则有

如果要求一个最小的$x$，那么一定要记得把答案$mod\frac{b}{gcd(a,b)}$

二元一次不定方程

因此这里有一道模板题【模板】二元一次不定方程，要求求出一个一般的二元一次不定方程的根，还添加了若干要求：

  • 如果无整数解，则输出-1
  • 否则如果没有正整数解，则输出最小的正整数x和最小的正整数y
  • 否则输出正整数解的数量、最小的x和最小的y、最大的x和最大的y

题解：

无整数解很简单，就是不满足裴蜀定理嘛，这里不多说

如何判断没有正整数解：由上面推出的公式可以发现，解x和y是此消彼长的关系，因此当x取最小正整数时，如果y非正了，那么就没有正整数解，否则就有，同理y也一样

如何求正整数解的个数：由上面推出的公式可以很容易得出最小的正整数x（即$把exgcd的x解mod\frac{b}{gcd(a,b)}$，但是记得如果得0了要加一个模数把它调正了），而最大正整数x就是y取最小正整数的时候的x，如法炮制一遍就行啦，最后的个数就是$\frac{x{min}-x\{max}}{\frac{b}{gcd(a,b)}}+1$（纸上写写就动了）

code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll gcd(ll a,ll b){
	return b==0?a:gcd(b,a%b);
}
ll x,y,a,b,c;
void exgcd(ll a,ll b){ //exgcd
	if(b==0){
		x=1,y=0;
		return;
	}
	exgcd(b,a%b);
	ll tx=x,ty=y;
	x=ty;
	y=tx-a/b*ty;
}
int main(){
	int t;
	scanf("%d",&t);
	while(t--){
		x=0,y=0;
		scanf("%lld%lld%lld",&a,&b,&c);
		ll g=gcd(a,b);
		if(c%g!=0){ //无解
			printf("-1\n");
			continue;
		}else{
		a/=g,b/=g,c/=g;
		exgcd(a,b);
		x*=c,y*=c; //求特解
		ll minx=(x%b+b)%b;
		if(!minx) minx+=b; //调正
		ll miny=(y%a+a)%a;
		if(!miny) miny+=a; //同上
		ll maxy=(c-(a*minx))/b;
		ll maxx=(c-(b*miny))/a;
		ll cnt=(maxx-minx)/b+1;
		if(cnt) printf("%lld %lld %lld %lld %lld\n",cnt,minx,miny,maxx,maxy);
		else printf("%lld %lld\n",minx,miny);
		}
	}
	return 0;
}

线性求逆元

问题

给定$n,p$，求$[1,n]$中所有整数在模$p$意义下的乘法逆元

算法

代码

对应题目【模板】乘法逆元

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int inv[3000005];
int main(){
	int n,p;
	scanf("%d%d",&n,&p);
	inv[0]=0,inv[1]=1;
	printf("1\n");
	for(int i=2;i<=n;i++){
		inv[i]=(ll)(p-p/i)*inv[p%i]%p;
		printf("%d\n",inv[i]);
	}
	return 0;
}

阶乘逆元

问题

如何求出$[1,n]$各数阶乘的逆元呢？

算法

这样一来，就可以倒着递推一遍就可以了

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
inline int get_inv(int x,int y){
    int base=x,ans=1;
    y-=2;
    while(y>0){
        if(y&1) ans=(ll)ans*base%mod;
        base=(ll)base*base%mod;
        y>>=1;
    }
    return ans;
}
int main(){
    int ifac[N],fac=1,mod;
    for(int i=2;i<=n;i++) 
        fac=(ll)fac*i % mod;
    ifac[n]=get_inv(fac, mod);
    for(int i=n-1;i>=0;i--) 
        ifac[i]=(ll)ifac[i+1]*(i+1)%mod;
}

逆元小结

  • 欧拉定理法（核心公式$a^{\varphi(p)}\equiv 1(mod\ p)$）
    • 衍生：$a^b=a^{b\%\varphi(p)}$（$a,p$互质）
    • 扩展欧拉定理：在$b\geq \varphi(p)时,$$a^b\equiv a^{b\%\varphi(p)+\varphi(p)}(mod\ p)$
  • 拓展gcd（核心公式$ax+by=bx’+(a\%b)y’$）
    • 衍生：对任意不定方程$ax+by=c$的解法
  • 线性求逆元（对固定模p）（核心公式p=ki+b）
    • 注意转正
  • 阶乘逆元（核心公式$(n-1)!\cdot n\cdot(n!)^{-1}=1$）

中国剩余定理及其扩展

中国剩余定理（CRT）

洛谷P1495 【模板】中国剩余定理(CRT)/曹冲养猪

定义

中国剩余定理（孙子定理，Chinese remainder theorem，简称CRT）是中国古代用来求解一次同余式方程组的一种算法，一次同余方程组的数学语言表达如下：

中国剩余定理成立的前提是：$b_i$两两互质，在满足此条件的前提下，可以通过构造得出该方程的解：

证明

image-20200612132505697

代码

1
2
3
4
5
6
7
8
9
10
11
inline ll crt(){
    ll M=1, res=0;
    rep(i,1,n) M*=b[i];
    rep(i,1,n){
        ll Mi=M/b[i];
        ll x,y;
        exgcd(Mi,b[i],x,y);
        res=(res+a[i]*Mi*ti)%M;
    }
    return (res%M+M)%M;
}

扩展中国剩余定理（EXCRT）

CRT只能够解模数两两互质的情形，那么对于一般的模数不互质的情况，又如何求解呢？

解法思路

image-20200612133108360

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// O1快速乘
inline ll multi(ll x,ll y,ll Mod){
	ll tmp=(x*y-(ll)((long double)x/Mod*y+1.0e-8)*Mod);
	return tmp<0 ? tmp+Mod : tmp;
}

// x = a[i] (mod b[i]) b[i]不互质
inline ll excrt(){
	ll ans=a[1],M=b[1];// 特判第一个方程
	for(int i=2;i<=n;i++){
		ll aa=M,bb=b[i];
		ll cc=(a[i]-ans%bb+bb)%bb;// aa*x=cc(mod bb)
		ll x,y;
		ll g=exgcd(aa,bb,x,y),bg=bb/g;
		if(cc%g) return -1;// 无解
		x=multi(cc/g,x,bg);// 防爆long long
		ans+=x*M;
		M*=bg;
		ans=(ans%M+M)%M;
	}
	return (ans%M+M)%M;
}

大步小步算法

大步小步算法（Baby Step Giant Step，简称BSGS 拔山盖世）是用来寻找离散对数（即模意义下的对数）的常用算法，具体来说，分为普通的BSGS算法和扩展BSGS，下面一一介绍

BSGS

对于形如$ax\equiv b(\mod p)$的同余方程已经有扩展欧几里得算法的通用解法了，然而，对于形如$a^x\equiv b(\mod p)$的同余方程又该如何求解呢？

模板题链接：洛谷P3846 [TJOI2007] 可爱的质数/【模板】BSGS

题意：求解最小的$x$，满足$a^x\equiv b(\mod p)$，其中，$2\le a,b<p<2^{31}$，且$p$为质数

算法思路：看到这样的指数形式同余式，容易联想到欧拉定理，即$a^{\varphi(p)}\equiv 1(\mod p)$，其中$p$为质数，同时，又有$a^0\equiv 1(\mod p)$，因此容易得到 $a^x$ 有一个长度为$a^{\varphi(p)}$的循环节，当$x>\varphi(x)$后，对应的等号右侧开始循环，因此，要求得最小解，只需要暴力枚举$x$到$\varphi(p)$即可，有因为$p$为质数，因此$\varphi(p)=p-1$，枚举范围为$[0,p-2]$，间复杂度$O(p)$，可以得零分（划掉）

我们发现这个数据范围线性复杂度是不够用的，需要一个根号算法，我们考虑将$x$分解开，令$x=i\cdot m-k$，上式化为$a^{i\cdot m}\equiv ba^k(\mod p)$，这样一来，可以通过预先哈希掉等号右侧的数的可能取值，然后枚举等号左侧的值找哈希表中是否存在过即可，显然这里$m$取$\sqrt p$时能够确保$i\cdot m+k$枚举到所有的$[0,p-2]$的值，因此复杂度为$O(\sqrt p)$

1
2
3
4
5
6
7
8
9
10
11
12
// a^x = b (mod p) <=> a^(im) = b * a^k (mod p)
inline ll BSGS(ll a,ll b,ll p){
	unordered_map<ll,ll> ma;// hash table
	ll m=ceil(sqrt(p));
	for(ll i=0,cur=b%p;i<=m;i++,cur=cur*a%p)
		ma[cur]=i;
	ll am=qpow(a,m,p);
	for(ll i=1,now=am;i<=m;i++,now=now*am%p){
		if(ma.count(now)) return (i*m-ma[now]+p)%p;// x
	}
	return -1;// no solution
}

EXBSGS

我们可以发现，BSGS是基于欧拉定理的暴力算法的优化，因此只能够处理$\gcd(a,p)=1$的情况，对于更加一般的离散对数求法，我们需要将BSGS进行扩展，得到扩展BSGS算法

模板题链接：洛谷P4195 【模板】扩展BSGS

题意：求满足$a^x\equiv b(\mod p)$的最小自然数$x$，其中，$a,p,b\le 10^9$

算法思路：此时若$\gcd(a,p)\neq 1$，就无法使用BSGS算法了，考虑令$g=\gcd(a,p)$，同时，上式可以转化为 $a\cdot a^{x-1}+p\cdot k\equiv b$，由裴蜀定理可知：这个式子有整数解当且仅当$g|b$，若有解，上式又可以转化为$a^{x-1}\cdot \frac{a}{g}+\frac{p}{g}\cdot k\equiv \frac{b}{g}$，即$a^{x-1}\cdot \frac{a}{g}\equiv \frac{b}{g}(\mod \frac{p}{g})$，令$a’=\frac{a}{g},(a’)^{-1}=(\frac{a}{g})^{-1}$，则得出$a^{x-1}\equiv \frac{b}{g}\cdot (a’)^{-1}(\mod \frac{p}{g})$，此时，可以将这个式子作为新的求解方程，用同样的方法求解，这是一个递归的过程，直到$\gcd(a,p)=1$时，就可以使用BSGS求解了

容易得出，最终的方程是形如下式的同余方程：

在得出可以用BSGS的方程后，先求出$\left(\prod_{i=1}^k{a_i’}\right)^{-1}$，然后BSGS求出对应的解，注意这里求出的解为$x’=x-k$的值，因此，如果有解则最终的解应该为$x’+k$

注意事项：

  • 在算法开始前，需要注意一下特判$b=1$的情况，直接输出$x=1$即可

  • 在EXGSGS中，如果遇到

    则$a^{x-k}=1$，即直接输出$k$即可

  • 好事多模（很重要）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// a^x = b (mod p) <=> a^(x-k) * a/prod(g) = b/prod(g) (mod p/prod(g))
inline ll EXBSGS(ll a,ll b,ll p){
	if(b==1) return 0;// a^x = 1 (mod p)
	ll g,k=0,ag=1;
	while((g=gcd(a,p))>1){
		if(b%g) return -1;// no solution
		k++,b/=g,p/=g,ag=ag*(a/g)%p;
		if(ag==b) return k;// a^(x-k)=1
	}
	ll inv,y;
	exgcd(ag,p,inv,y);// get inv of ag
	inv=(inv+p)%p;
	ll res=BSGS(a,b*inv%p,p);// BSGS
	return res==-1?res:res+k;// solution + k
}

素数筛法

埃氏筛法

算法

  • 从2开始枚举范围内每一个数，碰到是素数的就把他加入到素数表里，然后从这个素数开始把它的所有在范围内的倍数筛去即可
  • 优化：每次筛去时可以直接从$i^2$开始枚举，因为之前的数都已经在以前的$i$的循环里被筛去了
  • 时间复杂度为$O(nlogn)$，不是很给力，但够用了（至少比每次都判一遍素数那个方法靠谱）

代码

1
2
3
4
5
6
7
8
9
10
f[0]=1;
f[1]=1;// 0 and 1 is not prime
for(int i=2;i<=n;i++){
    if(f[i]==0){ //if i is prime number
        prime[++cnt]=i; // add
        for(int j=i*i;j<=n;j+=i){ //filter
            f[j]=1;
        }
    }
}

欧拉筛

算法

  • 同样从2开始枚举范围内每一个数，碰到素数就加入素数表里，之后再从当前的$i$开始筛去合数，筛掉的每个合数满足的条件为

    • 最小质因数*最大因数（非自己）=这个合数
    • 保证每一个合数都只会被这个方法筛掉一次

  • 时间复杂度为$O(n)$，非常给力

代码

1
2
3
4
5
6
7
8
9
10
f[0]=f[1]=1;
for(int i=2;i<=n;i++){
    if(f[i]==0)
        prime[++cnt]=i;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
        f[i*prime[j]]=1; //保证每个数被自己的最小质因数筛掉
        if(i%prime[j]==0) //保证线性，保证只被最小质因数筛一次
            break;
    }
}

代码中有一个最难理解的部分if(i%prime[j]==0 break;，这一句的意义在于保证了每个合数只被自己最小的质因数筛掉了唯一一次（保证了线性）

举个例子，设$i=p_1^{k_1}p_2^{k_2}…p_n^{k_n}$，如果上一步中没有及时break，那么j就会在达到$p_1$时继续往后累加，比如说打到了$p_2$，那么此时的合数的最小质因数就是$p_2$了，而$p_2>p_1$，如果在此刻筛掉i*prime[j]，则违背了每个合数只被它的最小质因数筛去一次的原则

欧拉函数

定义

欧拉函数$\varphi (x)$表示小于$x$的整数中与$x$互素的数的个数，特殊的$\varphi (1)=1$

递推式

一些性质

（这里只摘取一些比较简单的性质）

  1. 如果$p$为素数，则$\varphi (p)=p-1$（很容易理解，素数除了它本身之外小于它的数都与它互素）

  2. 欧拉函数在满足一定条件下是积性函数：

    当$m$,$n$为素数时，有$\varphi (m * n)=\varphi (m) * \varphi (n)$

等等。。。

求一个数的欧拉函数值（复杂度$O(\sqrt n)$）

直接根据定义求即可，给出代码如下:

1
2
3
4
5
6
7
8
9
10
// get phi[n]
inline ll getphi(ll n){
	ll res=n;
	for(ll i=2;i*i<=n;i++){
		if(n%i==0) res=res/i*(i-1);
		while(n%i==0) n/=i;
	}
	if(n>1) res=res/n*(n-1);
	return res;
}

利用欧拉筛法筛出区间$[1,n]$所有数的欧拉函数（复杂度$O(n)$）

在筛法中，有三种情况需要讨论：

  1. i是个质数：

    此时，phi[i]=i-1即可

  2. i%prime[j]==0：

    此时，phi[i]已经有包含了prime[j]的所有素因数，即phi[i]和phi[i*prime[j]]的推导式内的质因子种类是一样的，只不过i*prime[j]会比i多出一个prime[j]作为素因数，有下式成立：

    很容易看出

  3. i%prime[j]!=0：

    此时，phi[i]中没有pirme[j]这个质因子，但prime[j]又是满足i*prime[j]这个合数的最小质因子，我们得到下式

由情况3，可以得出欧拉函数是一个不完全积性函数的特性

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
inline void euler(int n){
    phi[1]=1; //欧拉函数特例
    f[0]=f[1]=1; //0和1都不是素数
    for(int i=2;i<=n;i++){
        if(f[i]==0){ //如果该数是素数
            prime[++cnt]=i; //计入素数数组
            phi[i]=i-1; //情况1
        }
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){ //欧拉筛核心代码
            f[i*prime[j]]=1; //筛去
            if(i%prime[j]==0){
                //情况2
                phi[i*prime[j]]=phi[i]*prime[j]; 
                break; //核心
            }else{
                //情况3，欧拉函数是对于素数的积性函数
                phi[i*prime[j]]=phi[i]*phi[prime[j]];
            }
        }
    }
}

快速分解质因数

在欧拉筛内部循环，加入如下代码

1
2
3
4
for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
    from[i*prime[j]]=prime[j];
    ...
}

记录每个数的最小质因数，随后查询一个数x的质因数时，只需要每次求出from[x]，然后继续求from[x/from[x]]，写个循环或者递归即可

练习题

洛谷 P2158 仪仗队

题解

自己在纸上画一画会发现题目意思就是找到一个$n*n$的矩阵里面有多少个点和坐标原点连起来的斜率不同，进而又可以转换成就其中有多少个点的坐标(x,y)满足x与y互素，因为如果不互素的话，就会出现(x/gcd(x,y),y/gcd(x,y))这样的点挡住他们

我们把矩阵沿着对角线切一刀，去掉对角线的一个点（2，2），剩下的就是算一个三角形内部有多少这样的点对，算好之后最终答案就是（2*ans+1），于是可以想到用欧拉函数来求

显然答案就是$\sum_{i=1}^{n-1}(\varphi (i))$

code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#include<bits/stdc++.h>
using namespace std;
int f[40005],prime[40005],ans[40005],cnt;
int main(){
  int n;
  cin>>n;
  if(n==1){
    cout<<0<<endl;
    return 0;
  }
  ans[1]=1;
  f[0]=f[1]=1;
  for(int i=2;i<=n;i++){ //欧拉筛欧拉函数
    if(f[i]==0){
      prime[++cnt]=i;
      ans[i]=i-1;
    }
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      f[i*prime[j]]=1;
      if(i%prime[j]==0){
        ans[i*prime[j]]=ans[i]*prime[j];
        break;
      }else{
        ans[i*prime[j]]=ans[i]*ans[prime[j]];
      }
    }
  }
  for(int i=1;i<n;i++){ //累加
    ans[i]=ans[i]+ans[i-1];
  }
  cout<<ans[n-1]*2+1<<endl;
  return 0;
}

数论分块（整除分块）

算法介绍

整除分块是用来求解形如下式的整除求和式的算法，朴素算法很容易想到，即遍历一遍$[1,n]$求解即可，时间复杂度为$O(n)$，而整除分块的处理能够使得复杂度上界优化到$O(\sqrt n)$

其实数论分块最终就是一个实用结论，其结论为：上式的整除求和式中$\lfloor\frac{k}{i}\rfloor$最多只有$O(\sqrt n)$种取值，且相同的取值总是连续地一块一块出现的，因此称之为分块，那么具体如何分块呢？

经过数学推导，可以知道，对于区间$[1,n]$内的一个左端点$l$，其整除值为$\lfloor\frac{k}{l}\rfloor$，那么就有一段连续的区间整除值与这一块相同，其区间右边界为$\min(n,\lfloor{\frac{k}{\lfloor\frac{k}{l}\rfloor}}\rfloor)$，当然，如果$\lfloor\frac{k}{l}\rfloor=0$，则右边界显然就是$n$（因为此后分母都比分子大了，答案就都是0）

证明

证明选自 题解： P2261 [CQOI2007]余数求和

image-20200612135223242

（但其实这个规律打个表也能发现）

代码

1
2
3
4
5
6
for(int l=,r;l<=n;l=r+1){
    int tp=k/l;
    if(tp) r=min(n,k/tp);
    else t=n;
   	ans+=(r-l+1)*tp;
}

就是这么简单

习题

P2261 [CQOI2007]余数求和（取余式转换+数论分块）

code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
ll n,k,ans;

int main(){
    read(n),read(k);
    for(ll l=1,r;l<=n;l=r+1){
        ll tp=k/l;
        if(tp) r=min(k/tp,n);
        else r=n;
        ans+=(l+r)*(r-l+1)/2*tp;
    }
    ans=n*k-ans;
    printf("%lld\n",ans);
    return 0;
}

P2260 [清华集训2012]模积和（数学推导+前缀和+数论分块）

code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
ll n,m,inv2,ans,inv6,x,y;

inline ll f(ll a,ll b){
	ll res=0;
	for(ll l=1,r;l<=a;l=r+1){
		ll tp=b/l;
		if(tp) r=min(b/tp,a);
		else r=a;
		res=(res%mod+(r-l+1)%mod*(r+l)%mod*inv2%mod*tp%mod+mod)%mod;
	}
	return (res%mod+mod)%mod;
}

inline ll sum(ll x){
	x%=mod;
	return (x%mod*(x+1)%mod*((x+x)%mod+1)%mod*inv6%mod+mod)%mod;
}

ll exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){x=1,y=0;return a;}
	ll res=exgcd(b,a%b,x,y);
	ll tx=x;x=y;y=tx-a/b*y;
	return res;
}

int main
	read(n),read(m);
	if(n>m) swap(n,m);
	exgcd(2,mod,x,y);
	inv2=(x+mod)%mod;
	x=0,y=0;
	exgcd(6,mod,x,y);
	inv6=(x+mod)%mod;
	
	ll fnn=f(n,n)%mod;
	ll fmm=f(m,m)%mod;
	ll fnm=f(n,m)%mod;

	ll aa=(n*n%mod-fnn%mod+mod)%mod;
	ll bb=(m*m%mod-fmm%mod+mod)%mod;
	ans=aa*bb%mod;
	ans=(ans-n*n%mod*m%mod+mod)%mod;
	ans=(ans+m*fnn%mod+mod)%mod;
	ans=(ans+n*fnm%mod+mod)%mod;
	for(ll l=1,r;l<=n;l=r+1){
		r=min(n/(n/l),m/(m/l));
		ll tp=(sum(r)%mod-sum(l-1)%mod+mod)%mod;
		ans=(ans-tp*(n/l)%mod*(m/l)%mod+mod)%mod;
	}
	printf("%lld\n",ans%mod);
    
	return 0;
}// https://www.luogu.com.cn/problem/P2260

P6583 回首过去（数论+数论分块+容斥）

题意： 给出一个正整数$n$，求有序整数对$(x,y)$的个数，满足$1\le x,y\le n$，且$\frac{x}{y}$可以表示为十进制有限小数

题解： 一道数论好题，首先能想到十进制有限小数的分母一定只含有$2,5$两个质因子，因此满足提议的分数一定可以写成$\frac{bc}{2^p5^qc}$的形式，这样一来，可以$[1,n]$枚举分母，然后去掉$2,5$的因子从而获得$c$，最后获得$b$即可算出答案，复杂段接近$O(n)$，只能拿80分

code 
1
2
3
4
5
6
for(int i=1;i<=n;i++){
    int tp=i;
    while(tp%2==0) tp/=2;
    while(tp%5==0) tp/=5;
    ans+=n/tp;
}

进一步观察上方的式子，发现他是一个类似整除分块的式子，即

其中，$f(x)$代表$[1,x]$中质因子只包含$2,5$的数的个数，这样一来，可以先预处理出$f(n)$，然后从小到大排个序，每次整除分块时从后往前寻找$f(\lfloor\frac{n}{i}\rfloor)$，更新答案，但是这里还需要注意一点，就是表达式中的$i$为不包含$2,5$两个因子的数，因此可以通过在整除分块时容斥来解决这一问题（减去$2,5$因子的数的个数，再加上$10$的因子的数的个数）

code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
ll n,cnt,a[maxn],ans;
int main(){
    read(n);
	for(ll i=1;i<=n;i*=2){
		for(ll j=1;i*j<=n;j*=5){
			a[++cnt]=i*j;
		}
	}
	sort(a+1,a+1+cnt);
	ll now=cnt;
	for(ll l=1,r;l<=n;l=r+1){
		ll b=n/l;
		r=n/b;
		ll tp=r-l+1;
		ll tp2=r/2-(l-1)/2;
		ll tp5=r/5-(l-1)/5;
		ll tp10=r/10-(l-1)/10;
		tp=tp-tp2-tp5+tp10;
		while(a[now]>b&&now) now--;
		ans+=tp*b*now;
	}
	printf("%lld\n",ans);
}