【学习笔记 / 长期更新】OI 中的数论

目录

  • Preface

    • 0.1 前言
    • 0.2 常用符号
    • 0.3 快速幂
  • 数论

    • 1 数论基础

      • 1.1 整除
      • 1.2 带余除法
      • 1.3 同余
      • 1.4 数论函数
    • 2 素数

      • 2.1 唯一分解定理
      • 2.2 判定素数
      • 2.3 Miller-Rabin 判定素数
        • 2.3.1 费马小定理
        • 2.3.2 二次探测定理
        • 2.3.3 Miller-Rabin
      • 2.4 素数筛
        • 2.4.1 埃氏筛
        • 2.4.2 欧拉筛 / 线性筛
        • 2.4.3 筛法的一些优化
      • 2.5 分解质因数
        • 2.5.1 朴素算法
        • 2.5.2 素数筛优化
        • 2.5.3 Pollard-Rho
      • 2.6 欧拉函数
        • 2.6.1 性质
        • 2.6.2 欧拉函数
        • 2.6.3 引理
        • 2.6.4 求一个数的欧拉函数
          • 2.6.4.1 朴素求法
          • 2.6.4.2 Pollard-Rho
        • 2.6.5 筛法求欧拉函数
    • 3 最大公约数

      • 3.1 最大公约数
      • 3.2 最小公倍数
      • 3.3 扩展欧几里得算法
        • 3.3.1 用 exgcd 求解的一个问题
        • 3.3.2 有理数取余:用 exgcd 求解的另一个问题
        • 3.3.3 exgcd 例题
    • 4 数论分块

      • 4.1 算法
      • 4.2 例题
      • 4.3 用数论分块解决的一个问题
  • Reference

- Preface

0.1 前言

本文意为作者从 0 开始学习数论,同时也对 OI Wiki 的某些内容做补充说明。

如果你看到有一些小标题没有内容,很正常,作者会填坑的。

同时本文用了大量 LATEX,加载公式较慢属正常现象。

都看到这了不点个赞吗 /kel

Page Views Count

0.2 常用符号

  1. 整除符号:xy,表示 x 整除 y

  2. 取模符号:xmody,表示 x 除以 y 的余数。

  3. 互质符号:xy,表示 xy 互质。

  4. 最大公约数:gcd(x,y)(x,y)

  5. 最小公倍数:lcm(x,y)[x,y]

  6. 阶乘符号:n!,表示 1×2×3××n,特别的,0!=1

  7. 求和符号:,表示特定条件下几个数的和,例如:

    • i=1ni 表示 1n 的和。
    • xn,xn1 表示 1n 里有几个数和它互质,即 φ(n)
  8. 求积符号:,表示特定条件下几个数的积,例如:

    • i=1ni 表示 n 的阶乘,即 n!
    • xdx 表示 d 的所有因数的乘积。
  9. 向上取整符号:x,表示大于等于 x 的最大的整数。

  10. 向下取整符号:x,表示小于等于 x 的最大的整数。

  11. 组合数:(xy)

  12. 整数集:Z,表示所有的整数。

  13. 自然数集:N,表示所有的自然数。

  14. 实数集:R,表示所有的实数。

  15. 存在:x:P(x),表示至少存在一个 x 使得 P(x) 的值为真。

0.3 快速幂 Θ(logn)

从熟悉的快速幂开始学起。

题意:给定 a,b,求 ab

朴素算法是 Θ(b) 的,较慢。

b 的二进制为 (ntnt1n2n1)2,考虑把 b 二进制分解,即 nt2p+nt12p1++n121+n020,其中 n{0,1}p 为非负整数,易得 p=log2n

又由幂的性质可得,ax+y=ax×ay

于是 ab 就可写为如下形式:

ab=ant2p+nt12p1++n121+n020=ant2p×ant12p1×an020

又有 x2y=xyxy=(xy)2,我们就可以在 log2n 的时间算完上式了。

  • 例如:

13=(1101)2=1×23+1×22+0×21+1×20

513=523×522×520=58×54×51=390625×625×5=1220703625

由于取模并不影响乘法的运算,所以这种方法也可以求 abmodk

Sample Code

// 0.3 qpow
template<class T = long long> T qpow(T a, T b,const T mod = b - 2){
	T ans = 1, base = a;
	while(b){
		if(b & 1) ans = ans * base % mod;
		base = base * base % mod;
		b >>= 1;
	}
	return ans;
}

- 数论

1 数论基础

1.1 整除

a,bZ,a0,且如果 qZ,使得 b=aq,那么我们就说 a 整除 b,或者 b 能被 a 整除,记作 ab,其中 ba 的倍数,ab 的因数。

否则我们就称 b 不能被 a 整除,记作 ab

性质:

  • ab=ab=ab=|a||b|
  • abbcac
  • abacx,yZ,a(xb+yc)
  • abbab=±a
  • m0,则 abmamb
  • b0,那么 ab|a||b|
  • a0,b=qa+c,那么 abac

以上性质的证明

对于整数 b0b 的约束只有有限个,0 是任何不等于 0 的整数的倍数。

对于整数 b0,1b 有四个显然因数为 ±1,±b,特别的,1 只有两个显然因数。

对于整数 b0 的真因数称为 b 非显然因数的因数。

  • b0,那么当 d 遍历 b 的因数的时候,bd 也在遍历 b 的因数。
  • b>0,那么当 d 遍历 b 的正因数的时候,bd 也在遍历 b 的正因数。

1.2 带余除法

a,bZa0,则有 b=qa+r,其中 q=ba,0r|a|

这里的 r 被称为最小非负余数,也可以写成 r=bmoda

其中 r=baq=baba

性质

  • 连续 a 个整数除以 a,余数一定且正好取到 0a1 各一个,其中有且仅有一个正好被 a 整除。
  • 任何一个整数除以 a,余数一定为 0a1 其中一个。

1.3 同余

1.4 数论函数

1.4.1 积性函数

指的是如果 x,yN,gcd(x,y)=1f(1)=1,那么 f(x×y)=f(x)×f(y)

特殊的,如果 x,yNf(1)=1,那么 f(x×y)=f(x)×f(y),则称 f(n)完全积性函数。

  • 如欧拉函数 φ(a×b)=φ(a)×φ(b) {gcd(a,b)=1}

2 素数

定义:若整数 b0,±1 只有四个显然因数的时候,b 是素数,否则 b 是合数。

±b 总是同为素数或同为合数。特别的,若无特殊说明,则素数指正素数

2.1 唯一分解定理

  • 合数总能分解成几个素数的乘积。

形式化的,设 a 为合数,p 为素数集,则

a=p1α1p2α2pnαn

  • 例如

24=2×2×2×3=23×31

2.2 判定素数 Θ(n)

通过素数的性质可知,对于素数 nkZ,k[2,n),有 kn,于是不难想到一个 Θ(n) 的朴素程序,用 2n1 的所有正整数对 n 进行试除,如果 nmodk=0,说明 p 除了显然因数外,还有 knk 的因数,因此判定为合数。

Sample Code

// 1.2.2 prime check O(n)
bool check_linear(int x){
	for(int i = 2; i < x; i++){
		if(x % i == 0) return false;
	}
	return true;
}

显然,当 n109 级别时,此种方法便会超时。

注意到 1.1 的一句话:b>0,那么当 d 遍历 b 的正因数的时候,bd 也在遍历 b 的正因数。

也就是说,当 k 访问到 n 的时候,如果此时还存在 nmodk=0,那么它就会在访问到 nk 的时候被判定为合数,因此,n 之后的判定完全是多余的,因此时间复杂度便可降到 Θ(n) 级别。

Sample Code

// 1.2.2 prime check O(sqrt(n))
bool check_sqrt(int x){
	for(int i = 2; i * i <= x; i++){ //注意这里是 <=x,不然可能把一些完全平方数放过去
		if(x % i == 0) return false;
	}
	return true;
}

2.3 Miller-Rabin 判定素数 Θ(klog3n)

有时候,要判定的数的级别是 1018 级别的,这时候 Θ(n)check_sqrt 就会 TLE,这时我们可以对其进行 Miller-Rabin 素数测试,来判断其是否为素数。

2.3.1 费马小定理(Fermat's Little Theorem)

有质数 p,则对于正整数 ap,有 ap11(modp)

考虑它的逆否命题,我们就有了一种判断素数的方法:

如果 ap11(modp),则 p 不是素数。

  • 注意 FLT 的逆定理是不成立的,即对于正整数 1a<p,有 ap11(modp),则 p 是质数,一个典型的例子是 Carmichael

2.3.2 二次探测定理

对于上面的 Carmichael 数,FLT 不能判断其素性。此时就可以用到二次探测定理

如果 p 是一个素数,0<x<p,那么 x21(modp) 的解只有两个:x1=1,x2=p1

证明:

x21(modp)x210(modp)(x+1)(x1)0(modp)

于是有 (x+1)(x1)p,故得出 1p1 的正整数只有 1p1 两个解,故原命题得证。

2.3.3 Miller-Rabin

于是我们有了一种判定素数的方法:Miller-Rabin。

具体步骤:

  1. 如果 n=1nmod2=0,直接返回 false;
  2. an1 化成 (an12)2,(an14)22 的形式,直到不能化为止,设最后的底数为 au
  3. 选择素数 a,计算 aumodn
  4. 不断运用二次探测定理,判断当 x21(modp) 成立时,是否 x=1x=n1,如果不是即探测失败,返回 false
  5. 如果满足条件,就判断 an11(modp) 是否成立,如成立,则这一轮 Miller-Rabin 检测结束。

一般来说素数 a 的选择在 long long 范围内选择 2,3,5,7,11,13,17,19,23 这九个数即可 100% 正确。

时间复杂度:Θ(klog3n)k 即为选择了几个素数进行检测。

Sample Code

// 1.2.3.3 Miller-Rabin

typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;

ll mul(ll a,ll b,ll mod){
// 快速乘,原理是用 ull 的溢出解决 ll 的溢出
    ll c=(ld)a/mod*b;
    ll res=(ull)a*b-(ull)c*mod;
    return (res+mod)%mod;
}

ll qpow(ll a,ll b,const ll p){
// 快速幂
	ll ans=1ll,base=a;
	while(b){
		if(b&1) ans=mul(ans,base,p);
		base=mul(base,base,p);
		b>>=1;
	}
	return ans;
}

ll ud[]={2,3,5,7,11,13,17,19,23}; // 素数集

bool MRtest(ll n){
	if(n%2==0 || n<3) return n==2; // 特判
	ll u=n-1,t=0;
	while(u%2==0) u>>=1,++t; // 分解为 a^u
	ll x,pre;
	for(auto p:ud){ // 选择素数
		if(p==n) return true; // 特判
		x=qpow(p,u,n);
		pre=x;
		for(int i=1;i<=t;i++){
			x=mul(x,x,n);
			if(x==1 && pre!=1 && pre!=n-1) return false;
			pre=x; // 二次探测
		}
		if(x!=1) return false; // 费马小定理
	}
	return true;
}

2.4 素数筛

通过 1.2.1 的素数判定,我们知道了如何判定一个数是否是素数,但是,如果要判断 1n 中所有的数是否是素数呢?

一个显然的做法是扫一遍区间,用根号做法或者 Miller-Rabin 判断每一个数是否是素数,但是 Θ(nn)Θ(nklog3n) 的复杂度显然太慢,于是我们引进了素数筛的思想。

2.4.1 埃氏筛 Θ(nloglogn)

唯一分解定理(2.1,每个合数都可以被分解成若干个质数的乘积,换句话说,每个合数都是它的质因数的倍数

于是我们可以得到这样一种筛法(埃拉托斯特尼筛法,简称埃氏筛):建立一个标记数组 vis,每次从小到大选择未被标记的数,把它的倍数都标记为合数,这样在结束后未被标记的数就是素数。

时间复杂度为 Θ(nloglogn),近乎线性。

Sample Code

// 1.2.4.1 Eratosthenes sieve
bool vis[1000005]; //标记此数是不是素数
vector<int> prime; //素数集
void Eratosthenes_sieve(int n){
	vis[1] = 1;
	for(int i = 2; i <= n; i++){
		if(vis[i] == 0){
			prime.push_back(i);
			for(int j = i + i; j <= n; j += i){
				vis[j] = 1;
			}
		}
	}
}

2.4.2 欧拉筛 / 线性筛 Θ(n)

Eratosthenes sieve 已经足够快了,但是碰到一些 n=107 的筛就显然会 TLE。继续改进算法。

我们注意到,在埃氏筛中,42 被筛了三遍(分别为 2,3,7),比较耗费时间,事实上,如果一个数只被被其最小的素因子筛一次,那么复杂度就会降到 Θ(n),此种方法被称作欧拉筛或线性筛,见代码:

Sample Code

// 1.2.4.2 Euler sieve
bool vis[1000005];
vector<int> prime;
void Euler_sieve(int n){
	vis[1] = 1;
	for(int i = 2; i <= n; i++){
		if(!vis[i]) prime.push_back(i);
		for(int j = 0; (size_t)j < prime.size(); j++){
			int it = prime.at(j);
			if(1ll * it * i > n) break;
			vis[1ll * it * i] = 1;
			if(i % it == 0){
				// i % pri[j] == 0
		        // 换言之,i 之前被 pri[j] 筛过了
		        // 由于 pri 里面质数是从小到大的,所以 i 乘上其他的质数的结果一定也是
		        // pri[j] 的倍数 它们都被筛过了,就不需要再筛了,所以这里直接 break
		        // 掉就好了
		        // From OI-Wiki
				break;
			}
		}
	}
}

2.4.3 筛法的一些优化

  • 空间优化:当空间限制很紧的时候,可以用 std::bitset 代替 bool 数组节省空间,但是这会使程序效率降低,可能会 TLE(以空间换时间)。

2.5 分解质因数

给定一个数 nN+,分解出它的质因数。

2.5.1 朴素算法 Θ(n)

根据唯一分解定理,一个合数总能被分解成几个素数的乘积,于是我们就有了一个类似根号时间内判定素数的程序,其时间复杂度也为 Θ(n)

Sample Code

vector<pair<int,int> > factor(ll x){ // 返回一个 pair<int,int> 类型的 vector
	vector<pair<int,int> > v; // first 存储素因子,second 存储素因子的次数
	for(ll i=2;i*i<=x;i++){
		if(x%i==0){
			int t=0;
			while(x%i==0) t++,x/=i;
			v.push_back(make_pair(i,t));
		}
	}
	if(x>1) v.push_back(make_pair(x,1)); // 如果 x 本身就是个素数
	return v;
}

2.5.2 素数筛优化 Θ(nlnn)

我们知道,n 之内的素数个数大约为 nlnn 个,我们如果事先把这些素数筛出来再加上上面朴素算法就能做到以上的复杂度。

代码和上文几乎一样,只是把 i 放到素数集里面循环,代码就不放了。

  • 实际上,筛出素数也需要 Θ(n) 的时间,所以如果 n 很大,还是得选择更优秀的算法,如下文介绍的 Pollard-Rho

2.5.3 Pollard-Rho Θ(n14)

2.6 欧拉函数 φ(n)

2.6.1 性质

首先有欧拉函数是积性函数,也就是说,如果 gcd(a,b)=1,有 φ(a×b)=φ(a)×φ(b)

由素数定义可知,如果 p 为素数,φ(p)=p1

2.6.2 欧拉函数(Euler's totient function)

凭借 2.6.12.1 唯一分解定理
于是我们就有了欧拉函数的计算方法,设合数 n=i=1spiki,则有

φ(n)=n×i=1spi1pi

我们来分析一下这个公式是怎么来的。

2.6.3 引理

p 为素数集,则 φ(pk)=pkpk1

这很好证明,pk 个数中,只有 pk1 个是 p 的倍数,其他数都不含 p 这个因子,故 φ(pk)=pkpk1

于是,根据欧拉函数的积性和唯一分解定理可知

φ(n)=i=1sφ(piki)=i=1spikipiki1=i=1spiki1(pi1)=i=1spikipi1pi=n×i=1spi1pi

2.6.4 求一个数的欧拉函数

2.6.4.1 朴素求法 Θ(n)

2.6.4.2 Pollard-Rho Θ(n14)

2.6.5 筛法求欧拉函数 Θ(n)

我们注意到,欧拉筛每次都会以它最小的素因子筛到这个数。设 n=pq,其中 p 为它最小的素因子,分类讨论:

  1. qmodp=0,这也就是说,q 中有 n 全部的素因子,于是

φ(n)=n×i=1spi1pi=p×q×i=1spi1pi=p×φ(q)

  1. qmodp0,这也就是说 p,q 互质,于是根据积性,有

φ(n)=φ(p)×φ(q)

Sample Code

#include <bits/stdc++.h>
using namespace std;

const int MAXN = 1000000 + 5;

int pri[MAXN],phi[MAXN],cnt;
bool is_prime[MAXN];

void get_euler(const int N){
	for(int i=1;i<=N;i++) is_prime[i]=1;
	is_prime[1]=0;
	phi[1]=1;
	cnt=0;
	for(int i=1;i<=N;i++){
		if(is_prime[i]){
			pri[++cnt]=i;
			phi[i]=i-1; //i为素数,phi[i]=i-1
		}
		for(int j=1;j<=cnt;j++){
			if(i*pri[j]>N) break;
			is_prime[i*pri[j]]=0;
			if(i%pri[j]) // 注意分辨 phi 和 pri
				phi[i*pri[j]]=phi[i]*phi[pri[j]]; // 情况1
			else{
				phi[i*pri[j]]=phi[i]*pri[j]; // 情况2
				break;
			}
		}
	}
}

int n,T;

int main(){
	get_euler(MAXN-5);
	cin>>T;
	while(T--){
		cin>>n;
		cout<<phi[n]<<'\n';
	}
}

3 最大公约数

3.1 最大公约数 gcd(a,b)

定义正整数 a,b 的最大公约数为他们共同的约数中最大的一个(Greatest Common Divisor),简写为 gcd

对于如何快速地求出 gcd(a,b),有欧几里得算法

  • 欧几里得算法:对于正整数 a,b,有 gcd(a,b)=gcd(b,amodb)

下面给出这个算法的证明过程(个人觉得 OI-Wiki 上的有些不严谨):

Proof

  • 等式 1:如果 abba,有 a=±b。(1.1 整除

  • 等式 2abacx,yZ,a(xb+yc)。(1.1 整除

  • 等式 3aZ,nN+,amodn=anan。(1.2 带余除法

  • 推论 1dadbdgcd(a,b)

    • 很显然 a,b 的公约数里有 a,b 共同的因子,然后 gcd(a,b) 显然也有他们共同的因子,根据 唯一分解定理gcd(a,b) 要么比 d 含有的因子数多,要么他们共同含有的因子 gcd(a,b) 的幂次一定大于 d(不然不是最大公约数),于是 dgcd(a,b)

接下来我们来推出 gcd(a,b)gcd(b,amodb)

d=gcd(a,b),r=amodb

等式 3r=abab

于是根据 等式 2 ,由 dadbr=abab,推出 dr,即 damodb

推论 1 ,因为 dbdamodb,于是 dgcd(a,amodb)

gcd(a,b)gcd(b,amodb)


接下来我们来推出 gcd(b,amodb)gcd(a,b)

c=gcd(b,amodb),则有 cbc(amodb)

由带余除法,a=bq+r,其中 q=ab,r=amodb

于是根据 等式 2 ,由 cbc(amodb)a=bab+amodb,得出 ca

于是有 cacb 又因为 推论 1 ,得出 cgcd(a,b)

gcd(b,amodb)gcd(a,b)


于是根据 等式 1gcd(b,amodb)gcd(a,b)gcd(a,b)gcd(b,amodb)gcd(a,b)=gcd(b,amodb)

至此,欧几里得定理(或辗转相除法、gcd 递归定理) 证明完毕。

因为每次 amodb 的值至多变为 a 的一半,所以该算法的时间复杂度为 Θ(logn)

常见的几种写法:

// 以下默认 typedef long long ll 
ll gcd(ll a,ll b){ // 非递归版
    int tmp;
    while(b!=0){
        tmp=b;
        b=a%b;
        a=tmp;
    }
    return a;
}

ll gcd(ll a,ll b){ // 递归版
	if(b==0) return a;
	return gcd(b,a%b);
}

ll gcd(ll a,ll b){ // 递归版简化
	return b==0?a:gcd(b,a%b);
}

ll gcd(ll a,ll b){ // 位运算版
    while(b){
        a%=b;
        b^=a;
        a^=b;
        b^=a;
    }
    return a;
}

ll gcd(ll a,ll b){ // 位运算简化
	while(b^=a^=b^=a%=b);
	return a;
}

3.2 最小公倍数 lcm(a,b)

定义正整数 a,b 的最小公倍数为他们共同的倍数中最小的一个(Least Common Multiple),简写为 lcm

我们接下来来推一推 lcm 要如何快速地求。

唯一分解定理,可得

a=p1ka1p2ka2pskas,b=p1kb1p2kb2pskbs

因为 gcd 为二者的最大公因数,所以有

gcd(a,b)=p1min(ka1,kb1)p2min(ka2,kb2)psmin(kas,kbs)

同理,lcm 即为

lcm(a,b)=p1max(ka1,kb1)p2max(ka2,kb2)psmax(kas,kbs)

因为 ka1+kb1=min(ka1+kb1)+max(ka2+kb2),所以 gcd(a,b)×lcm(a,b)=a×b

也就是说,只要求出了两个数的 gcd,两个数的 lcm 就能 Θ(1) 求出来,容易看出,求解 lcm 的时间复杂度也为 Θ(logn)

ll lcm(ll a,ll b){
	return a/gcd(a,b)*b; // 防止 a*b 爆 ll,所以先 /gcd 再 *b
}

3.3 扩展欧几里得算法 exgcd

扩展欧几里得算法(Extended Euclidean algorithm),简写为 exgcd,是用来求解形如 ax+by=gcd(a,b) 的一组可行解。

接下来我们来分析一下这个算法:

首先,当 b=0 时,显然有 gcd(a,0)=a,此时 {x=1y=0 是方程的一组可行解。

ax1+by1=gcd(a,b)bx2+(amodb)y2=gcd(b,amodb)

于是根据 欧几里得定理,显然有 ax1+by1=bx2+(amodb)y2

于是

ax1+by1=bx2+(amodb)y2ax1+by1=bx2+(abab)y2ax1+by1=bx2+ay2baby2ax1+by1=ay2+b(x2aby2)x1=y2,y1=x2aby2

gcd 递归求解即可。当 b=0 时,带入上文的特解返回。

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

3.3.1 用 exgcd 求解的一个问题

题意:给定 a,b,求出 ax1(modb)x 的最小正整数解,保证有解。

问题转化:ax1(modb)ax+by=1,其中 y 是我们引进的一个负整数

明显地,exgcd 可以用来求解形如 ax+by=gcd(a,b) 的解,现在的问题是如何把 1gcd(a,b) 联系起来。

引理ax+by=m 有整数解的必要条件是 gcd(a,b)m

Proof:因为 gcd(a,b)agcd(a,b)b,所以 gcd(a,b)(ax+by),即 gcd(a,b)m

于是再联系到题目上,就有 gcd(a,b)1,因为 gcd(a,b) 显然是个正整数,所以 gcd(a,b) 只能等于 1 了。

这也就是说,a,b 互质,于是 gcd(a,b)=1,正好符合 exgcd 的原始形式,所以 exgcd 就可以在 Θ(logn) 的时间来求解。

最小正整数解

因为 ax+by=1,所以 ax+by+kabkab=1,所以 a(x±kb)+b(yka)=1

这也就是说,如果 x>0,那么就把 x 一直减 b,直到不能减为止,反之亦然。

这在 C++ 的运算中,用 x = (x % b + b) % b 即可解决,后面的 + b) % b 是为了让 x 从负数变为正数。

Sample Code

#include<bits/stdc++.h>
using namespace std;
using ll = long long;

ll a,b,x,y;

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

int main(){
	cin>>a>>b;
	exgcd(a,b,x,y);
	x=(x%b+b)%b;
	cout<<x<<'\n';
	return 0;
}

3.3.2 有理数取余:用 exgcd 求解的另一个问题

  • 给定 c=ab,求 cmodp 的值。

这个值被定义为 bxa(modp) 的解。

我们来分析一下。

因为同余有性质: 如果 ba(modp),那么两边乘上一个数之后,还是不变: b×da×d(modp),带进原式,就变成了 bxa(modp)

慢着,我们先看上一个问题,他求解的是 bx11(modp),然后两边同乘 a,变成了 b(ax1)a(modp)

于是有 x=ax1,问题解决。

等等,那无解情况呢?我们分情况讨论一下:

  1. bmodp=0,即 bp 的倍数:

    • amodp=0,即 a 也为 p 的倍数,于是原同余式有无数解。
    • amodp0,即 a 不为 p 的倍数,于是原同余式无解(因为等式左边 modp 永远为 0)。
  2. bmodp0,即 b 不为 p 的倍数:

    • b,p 互质,gcd(b,p)=1,于是根据 3.3.1 中的一个结论,原方程必定有解。

总结:如果 bmodp=0,则原方程无解。

Sample Code

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int p = 19260817;

inline ll read(){
	char ch=getchar(); ll x(0),f(0);
	for(; !isdigit(ch); ch=getchar()) f|=(ch=='-');
	for(;  isdigit(ch); ch=getchar()) x=((x<<1)+(x<<3)+(ch^48))%p;
	return f?-x:x;
}

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

ll a,b,x,y;

int main(){
	a=read(); b=read();
	if(b==0)
		return puts("Angry!")&0;
	exgcd(b,p,x,y);
	ll ans=(a*x%p+p)%p;
	printf("%lld\n",ans);
}

3.3.3 exgcd 例题

4 数论分块

4.1 算法

给定 nN+, 求 i=1nni 的值。

朴素算法是 Θ(n) 的,考虑优化。

有结论:对于 ni=nj1ijn,有 j 的最大值为 ni

证明:设 k=ni,显然有 kni

nknni=i=i

j=maxni=nji=nk=nni

容易证得,符合要求的区间约有 n 种。

时间复杂度 Θ(n)

Sample Code

ll H(int n){
	ll res = 0;
	int l = 1, r;
	while(l <= n){
		r = n / (n / l);
		res += 1ll * (r - l + 1) * (n / l);
		l = r + 1;
	}
	return res;
}

4.2 例题

4.3 用数论分块解决的一个问题

P2261 [CQOI2007]余数求和 为例。

给出 n,k,求 i=1nkmodi

题目转化:由带余除法一章,我们知道 kmodi=ki×ki

于是问题就变成了 i=1nki×ki

G(n,k)=i=1nki×ki=n×ki=1ni×ki

其中 i=1ni×ki 部分可以数论分块。因为 ki 相等,由乘法分配律可得 l=1rl×kl=kl×l=1rl。后面的部分可以数论分块顺带解决掉。于是就做完了。

注意这题和上面那题不同。这题 n,k 是分开的,所以当 kl=0 时要特判,不然会 RE。

以及 r 要和 nmin

复杂度:Θ(n)

ll H(ll n, ll k){
	ll ans = k * n;
	int l = 1, r;
	while(l <= n){
		if(k / l != 0) r = k / (k / l);
		else r = n;
		chkmin(r, n);
		ans -= 1ll * (l + r) * (r - l + 1) / 2 * (k / l);
		l = r + 1;
	}
	return ans;
}

- Reference

posted @   TheSky233  阅读(521)  评论(1编辑  收藏  举报
相关博文:
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
点击右上角即可分享
微信分享提示