基础数论算法汇总

乘法逆元

给定 n 个正整数 ai,求它们在模 p 意义下的乘法逆元。

逆元是模意义下的倒数,能够将模意义下无法直接计算的除法转化为乘法。先来总结一下常用的求单个逆元的方法:

扩展欧几里得

O(logn) 地求一个数的逆元,要求 a,p 互质即可(p 为模数),原理为解线性同余方程 ax1(modp)

void exgcd(int a, int b, int& x, int& y) {
    if (b == 0) return x = 1, y = 0, void();
    exgcd(b, a % b, y, x);
    y = y - a / b * x;
}
exgcd(a, p, inv, temp);

快速幂

根据费马小定理,当 p 为质数时,a 在模 p 意义下的逆元是 ap2,可直接用快速幂 O(logn) 求解。注意: 扩展欧几里得的条件是 a,p 互质即可,但快速幂法要求 p 必须是质数!

下面两种方法用于求多个数的逆元:

线性求逆元

显然 111(modp),我们从 1 开始向后迭代求逆元。设 k=pij=pmodi,有 p=ki+j。所以 ki+j0(modp)。两边同乘 i1,j1,得到 kj1+i10(modp),所以要求的 i1=kj1。因为 j=pmodi<i,而在迭代到 i 时,[1,i) 的逆元都已经求出,则可以直接根据 j1 推出 i1 的值。

代码实现比较简单:

inv[1] = 1
for i in range(2, n + 1):
    k = p // i
    j = p % i
    inv[i] = (p - k) * inv[j] % p # p-k 避免出现负数

线性求任意数逆元

当要求逆元的数不连续时(如本题中的情况),可以通过前缀积处理。

//求前缀积 s[i]
s[0] = 1;
for (int i = 1; i <= n; i++) s[i] = (long long)s[i - 1] * i % p;
// 费马小定理求 n 个数积的逆元(也就是逆元的积)
sv[n] = qpow(s[n], p - 2, p);
// 逐个消去 a[i] 的逆元,得到逆元前缀积
for (int i = n; i > 1; i--) sv[i - 1] = (long long)sv[i] * a[i] % p;
// 对于每个逆元前缀积,消去 a[1...i-1] 的逆元积,得到 a[i] 的逆元
for(int i = 1; i <= n; i++) inv[i] = (long long)sv[i] * s[i - 1] % p;

这个题也可以直接通分:设 x=ai,原式 i=1nkiai(modp)=i=1nkixaix,其中 xai 可以通过预处理前后缀积求出,求一次 x 的逆元即可。

裴蜀定理

a,b 为不全为 0 的整数,ax+by=c 有整数解当且仅当 gcd(a,b)|c。证明如下:设 d=gcd(a,b),显然 d|a,d|b,取 x,yZ,有 d|ax,d|by,所以 d|(ax+by),得证。容易推广到多个整数的情况。

此题中,由裴蜀定理的推广得,gcd(A1,A2An)|S,取 S 为最小公约数即可。

def gcd(a, b):
    if b == 0:
        return a
    return gcd(b, a % b)

n = int(input())
a = list(map(abs, map(int, input().split())))
result = a[0]

for i in range(1, len(a)):
    result = gcd(result, a[i])

print(result)

扩展欧几里得

【模板】同余方程

求关于 x 的同余方程 ax1(modb) 的最小正整数解。根据取模的性质,这个方程相当于 ax+by=1,其中 y 为负数,形式类似于扩展欧几里得的经典形式 ax+by=gcd(a,b)

方程 ax+by=m 有整数解的必要条件是 gcd(a,b)|m,此处 m=1,所以有 gcd(a,b)=1=m。也就是说方程 ax+by=1ax+by=gcd(a,b),我们可以直接使用扩展欧几里得求解。

假设我们已知一组整数 x1,y1 使得 bx1+(amodb)y1=gcd(a,b) 成立,则有 ax+by=bx1+(amodb)y1 成立。由取模的性质得,ax+by=bx1+(abab)y1=ay1+b(x1+aby1),得到一组解 x=y1,y=x1+aby1。这样,只要知道了一组 x1,y1,就能够得到一组 x,y

那么怎么求 x1,y1 呢?我们发现,定义 x1,y1 的方程 bx1+(amodb)y1=gcd(a,b)ax+by=gcd(a,b) 形式相同,可以用同样的方法,令 ax+by=gcd(a,b) 中的 a=b,b=amodb,再去找一组 x2,y2 满足 bx2+(amodb)y2=gcd(a,b)。知道了 x2,y2,就知道了 x1,y1……依此类推。

类似欧几里得算法,递归边界是 a=gcd(a,b),b=0,又 axn+byn=gcd(a,b),所以 xn=1,yn 取任意值(为避免溢出,一般取 0),递归回去求出一开始的 x,y 即可。

此时的答案尚不符合题目“最小”的要求,需要进一步处理。设 a=m+bnn 极大,于是原方程为 (m+bn)x+by=1mx+b(nx+y)=1 使得 x 最小,所以直接 xmodb即可,注意调整范围避免出现负数。

void exgcd(int a, int b, int &x, int &y) {
    if (!b) return x = 1, y = 0, void();
    exgcd(b, a % b, y, x);
    y = y - (a / b) * x;
}
exgcd(a, b, x, y);
x = (x % b + b) % b; // 调整答案范围

【模板】二元一次不定方程:没有整数解当且仅当 gcd(a,b)c,直接输出 -1

用 exgcd 解方程 ax+by=gcd(a,b) 得到一组特解 x0,y0。对原方程变形得到 axcgcd(a,b)+bycgcd(a,b)=c,于是有 ax+by=c 的一组特解 x1=xcgcd(a,b),y1=ycgcd(a,b)

接下来,由「同余方程」一题中的技巧,用 x=(x%b+b)%b, y=(y%a+a)%a; 求出 x,y 的最小取值 xmin,ymin,因为 x 增大 y 减小,可以直接分别带入原方程,求出 ymax,xmax。如果 xmax,ymax<0,则没有正整数解,直接输出两个最小值即可;否则有正整数解,分别输出最大值、最小值和解的个数。

扩展欧拉定理

abmodm,b10200000

首先引入三种可以通过取模缩小幂指数的方法。

  1. 费马小定理:当 a,pZ, p 为质数且 pa 时,ap11(mod p),所以有 ababmod(p1)(mod p)
  2. 欧拉定理:当 a,mZ, a,m 互质时,aφ(m)1(mod m),其中欧拉函数 φ(n)[1,n] 中与 n 互质的数的个数。由此, ab=abmodφ(m)(mod m)
  3. 扩展欧拉定理:当 a,mZ 时,
    ab={ab,bφ(m) abmodφ(m)+φ(m),b>φ(m)(mod m)

此题没有给出 a,m 间的互质关系,只能采用扩展欧拉定理解决。只需求出 φ(m) 即可,我们可以采取分解质因数的方法,结合欧拉函数公式
φ(n)={1,n=1 ni=1k(11pi),n1
O(m) 地求解。公式里的 k,pi 来自 n 的唯一分解 n=i=1kpi, pi 为质数。

int phi = m;
for (int p = 2; n > 0; p++) {
    if (n % p != 0) continue;
    phi -= phi / p; // 公式中 (1-1/p) 的变形
    while (n % p == 0) n /= p;
}

下面延伸介绍一下欧拉函数的知识。除了公式,还有五种方法可以计算欧拉函数:

  1. p 为质数,φ(p)=p1
  2. p 为质数,n=pk,有 φ(n)=pkpk1,即 [1,n] 去掉了 p,2p,3ppk1p 共计 pk1 个数;
  3. 积性:若 a,b 互质,φ(ab)=φ(a)φ(b),因为此时 a,b 的质因子 {pi},{qi} 均不相同,所以 φ(ab)=abi=1k(11pi)j=1l(11qi)=ai=1k(11pi)×bj=1l(11qi)=φ(a)φ(b)
  4. 不完全积性:若 p 为质数,np 的倍数,则 nnp 的质因数种类上完全相同,有 φ(np)=npi=1k(11pi)=ni=1k(11pi)×p=φ(n)×p
  5. n 为奇数,则 φ(2n)=φ(n),由方法 1,3 可推得。

结合方法 1,3,4,我们可以用线性筛的方法 O(n) 的计算 φ(i),实现如下:

for (int i = 2; i <= n; i++) {
    if (!isNotPrime[i]) {
        primes.push_back(i);
        phi[i] = i - 1; // 方法 1
    }
    for (int j = 0; j < primes.size() && i * primes[j] > n; j++) {
        isNotPrime[i * primes[j]] = true;
        if (i % primes[j] == 0) {
            phi[i * primes[j]] = phi[i] * primes[j]; // 方法 4
            break; // 线性筛原理:只让合数被「最小」质因数筛掉
        }
        phi[i * primes[j]] = phi[i] * phi[primes[j]]; // 方法 3
    }
}
posted @   XYukari  阅读(28)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
历史上的今天:
2023-11-08 C语言变量分类
点击右上角即可分享
微信分享提示