Skip to content

扩展欧几里得算法

用于求解 ax+by=gcd(a,b) 的一组可行解。

ax1+by1=gcd(a,b)=gcd(bmoda,a)=(bmoda)x2+ay2

将上式展开:

ax1+by1=(bba×a)x2+ay2

因为只要求找到一组合法解,所以可以构造:

{y1=x2x1=y2ba×x2

所以,在辗转相除法求解 gcd(a,b) 的时候,递归迭代 x,y 即可得到一组构造解。

递归边界是 gcd(0,a),此时 0×x+a×y=gcd(0,a)=a

解得:y=1,xZ,但是根据不同的 x 的取值,迭代结果的值域也是不同的,可以证明当 x=0 时,迭代结果的值的绝对值是不超过 max(a,b) 的。

感性地理解,每次 x1=y2ba×x2 累乘了 ba,那么最终就会不超过 max(a,b) 因为是一路乘上去的,y2ba×x2 不会使绝对值变大。

那么,x=z0,结果值域可以近似认为扩大了 z 倍。

所以取 x=0,y=1 反代即可。

上述过程只是求出了 ax+by=gcd(a,b) 的一组特解。

cpp
void exgcd(int a, int b, int &x, int &y) {
    int t;
    if (!b) {
        x = 1, y = 0;
    } else {
        exgcd(b, a % b, x, y);
        t = x;
        x = y;
        y = t - (a / b) * y;
    }
}

考虑通解:

因为 x 确定时,y 的解是唯一的,所以我们考虑 x 的通解,y 同理。

已知 ax+by=gcd(a,b),令 x=x+δa(x+δ)+by=gcd(a,b)

δ 合法,当且仅当 y 有解,可得:by=gcd(a,b)a(x+δ),那么 bgcd(a,b)a(x+δ)

ax+by 替换 gcd(a,b) 整理式子可得:baδ

所以 δlcm(a,b)aδbgcd(a,b)

综上,原方程关于 x 的通解为:xx(modbgcd(a,b))

y 同理,可得 yy(modagcd(a,b))

求解二元一次方程整数解

裴蜀引理:ax+by=c 有整数解当且仅当 gcd(a,b)c

广义裴蜀引理:i=1naixi=d 有整数解当且仅当 gcd(a1,,an)d

根据「裴蜀引理」,二元一次方程 ax+by=c 若有解,则可以表示成:

ax+by=cgcd(a,b)×gcd(a,b)

使用 exgcd 求出 ax+by=gcd(a,b) 的一组特解,乘以 cgcd(a,b) 得到 ax+by=c 的一组特解:x,y

a(x+d)+b(ye)=cad=beagcd(a,b)d=bgcd(a,b)ed=t×bgcd(a,b),e=t×agcd(a,b)

代入原方程,可得:

原方程的通解为:{x=x+t×bgcd(a,b)y=yt×agcd(a,b)

求解线性同余方程

对于线性同余方程 axb(modc),可以改写成 ax+cy=b 的形式,看作二元一次方程求整数解即可。

cpp
int tyfc(int a, int b, int c) {
    b %= c;
    int x, y;
    if (b % gcd(a, c))
        return -1;
    exgcd(a, c, x, y);
    int g = gcd(a, c);
    x *= b / g;
    c /= g;
    x = (x % c + c) % c;
    return x;
}

欧拉定理

aφ(b)1(modb),gcd(a,b)=1

费马小定理:ap11(modp),pP

费马小定理是欧拉定理在 bP 时的特殊情况。

扩展欧拉定理

ab{abmodφ(m)gcd(a,m)=1abgcd(a,m)1,b<φ(m)a(bmodφ(m))+φ(m)gcd(a,m)1,bφ(m)(modm)

用于大指数取模。

欧拉降幂

aaa...modm,反复使用 exeuler,对指数不断modφ(m)

结论是一个数自取 O(logn)φ(n) 就会降到 1

简单证明:

  • φ(x),x2 是偶数
  • φ(x)x2,xmod2=0

所以至少是两次缩小一半规模。

卢卡斯定理

对于质模数 p,有 (nm)(npmp)×(nmodpmmodp)(modp),规定 n<m,(nm)=0

Lucas 定理揭示了,组合数 (nm) 对质数取模的结果等价于其 n,mp 进制意义下每一位的 (nimi)modp 之积。

即:i=1(nimi)modp

对于 (xy),x,y<p,预处理阶乘及其逆元,即可 O(1) 询问。

所以时间复杂度:O(logpn+p)

cpp
int C(int x, int y, int mod) {
    auto calc = [](int x, int y, int mod) -> int {
        if (x < y)
            return 0;
        if (x == 0 || y == 0)
            return 1;
        return fac[x] * inv[y] % mod * inv[x - y] % mod;
    };
    int res = 1;
    while (x && y) {
        int X = x % mod, Y = y % mod;
        res = res * calc(X, Y, mod) % mod;
        x /= mod;
        y /= mod;
    }
    return res;
}

扩展卢卡斯定理

若模数不为质数,令 m=i=1spiai

根据 crt,等价于求解 {(nk)r1(modp1a1)(nk)r2(modp2a2)(nk)rs(modpsas)

考虑求解 (nk)modpa

(nk)=n!k!(nk)!

vp(n) 表示 n 的质因数分解中质因子 p 的指数,wp(n) 表示 n 除尽 p 后的值。

则原式可以写成 pvp(n!)vp(k!)vp((nk)!)×wp(n!)wp(k!)wp((nk)!)

此时,问题转换成求 vp(n)modpawp(n)modpa

所以 exlucas 本质是解决了阶乘对质数幂取模的问题。

Wilson 定理:(p1)!1(modp),pP

推广:1im,imi±1(modm)

推论:1ipa,ipai{1p=2a21otherwise(modm)

对于 p

wp(n!)(1)np×(nmodp)!×wp(np!)(modp)

可以发现,实际上 lucas 就是模数为 p 的特殊情况。

对于 pa

wp(n!)(1)npa×1i(nmodpa),ipai×wp(np!)(modpa)

预处理 1i(nmodpa),ipai,递归求解,容易得到时间复杂度:预处理 O(pa),单次询问 O(logpn)

cpp
int phi(int x, int y) { return x * (y - 1) / y; }
int nn, a[N], b[N];
int mod[N], mem[N];
int intchina() {
    int i, prod = 1, s = 0;
    for (i = 1; i <= nn; i++)
        prod *= a[i];
    for (i = 1; i <= nn; i++) {
        int t = qz(prod / a[i], phi(a[i], mod[i]) - 1, a[i]);
        s = (s + b[i] * (prod / a[i]) * t) % prod;
    }
    return s;
}
int n, m, p;

void init(int p, int mod) {
    int i;
    mem[0] = 1;
    for (i = 1; i <= mod; i++) {
        mem[i] = mem[i - 1];
        if (i % p)
            mem[i] = mem[i] * i % mod;
    }
}
int l(int x, int p) {
    int res = p, ans = 0;
    while (res <= x) {
        ans += x / res;
        res *= p;
    }
    return ans;
}
int g(int x, int mod) {
    if (x <= mod)
        return mem[x];
    if ((x / mod) % 2 == 0)
        return g(x % mod, mod) % mod;
    return (-g(x % mod, mod) % mod + mod) % mod;
}
int f(int x, int p, int mod) {
    if (x <= 1)
        return 1;
    return f(x / p, p, mod) * g(x, mod) % mod;
}
int C(int n, int m, int p) {
    for (int i = 2; i <= p / i; i++) {
        int res = 1;
        while (p % i == 0) {
            res *= i;
            p /= i;
        }
        if (res > 1) {
            a[++nn] = res;
            mod[nn] = i;
        }
    }
    if (p > 1) {
        a[++nn] = p;
        mod[nn] = p;
    }
    for (int i = 1; i <= nn; i++) {
        int base = l(n, mod[i]) - l(m, mod[i]) - l(n - m, mod[i]);
        int aa = qz(mod[i], base, a[i]);
        int res1, res2;
        init(mod[i], a[i]);
        res1 = f(n, mod[i], a[i]),
        res2 = f(m, mod[i], a[i]) * f(n - m, mod[i], a[i]) % a[i];
        aa = aa * res1 % a[i];
        aa = aa * qz(res2, (a[i] - a[i] / mod[i]) - 1, a[i]);
        b[i] = aa;
    }
    return intchina();
}

BSGS

bsgs 用于求解指数同余方程,即:

axb(modc),gcd(a,c)=1

其基本原理是任意一个 [0,n] 的整数都可以表示成 ak+b,其中 a[0,n1k],b[0,k)

枚举 ab 的待选项只有 k 个,预处理 b,那么求解的时间复杂度即为 O(n1k+k)k=n 时,取到近似最优解,时间复杂度:O(n)

具体到原问题中,首先确定 n,根据鸽巢原理,易得 x<c,axac(modc),所以,若 axb(modc)[0,c) 中无解,那么就无解。

aAc+Bb(modc)aAc×aBb(modc)

gcd(a,c)=1,即 aAcmodc 存在逆元时,aB(aAc)1×b(modc)

预处理 a0,...,an1,查询求解 B 即可。

所以可以发现,bsgs 建立在 gcd(a,m)=1 的前提上。

有求逆元和 map 询问的存在,时间复杂度:O(nlogn)

小优化:

考虑将 aAc+B 写成 aAc+cB 的形式,那么可以将原式写成 aAcBb(modc)aAcb×aB(modc),A[1,n1k+1]

预处理 b×aB,若将 hashmap 的时间复杂度视作 O(1),则时间复杂度为:O(n)

省去了求逆元的时间复杂度。(但是上式仍建立在逆元存在的基础上)

或者实际上本质上是求 O(n) 个元素的逆元,朴素也可以做到 O(n+logn) 求解。

cpp
struct BSGS {
    unordered_map<int, int> mp;
    int bsgs(int a, int b, int p) {
        if (b == 1)
            return 0;
        b %= p;
        int sqt = ceil(sqrt(p));
        int A = b;
        mp.clear();
        for (int i = 0; i < sqt; i++) {
            mp[A] = i;
            A = A * a % p;
        }
        int base = qz(a, sqt, p);
        A = 1;
        for (int i = 1; i <= sqt; i++) {
            A = A * base % p;
            if (mp.count(A))
                return i * sqt - mp[A];
        }
        return -1;
    }
};

扩展 BSGS

用于求解:

axb(modc),gcd(a,c)1

因为 gcd(a,c)1 逆元不存在,不能使用 bsgs 求解。

gcd(a,c)1,根据 exgcd 原同余方程可以写作 agcd(a,c)×ax1bgcd(a,c)(modcgcd(a,c)),其中 gcd(agcd(a,c),cgcd(a,c))=1,所以逆元存在,可以将问题转换成 ax1b(modc) 的子问题,依次类推,直到 gcd(a,c)=1 为止。

此时,原式为:ayA×axybA(modcA)

采用 bsgs 求解后,x+y 即为原同余方程的解,注意此时解的范围是 [y,+),所以要特殊判断 [0,y) 是否存在解。

时间复杂度:O(n+logn) 与 bsgs 相同。

cpp
struct ExBSGS {
    int exbsgs(int a, int b, int c) {
        int sqt, base, i, A = 1;
        hashmap.clear();
        for (int i = 0; i < 32; i++) {
            if (A % c == b)
                return i;
            A = (A * a) % c;
        }
        int d = 1, cnt = 0;
        while (gcd(a, c) != 1) {
            int g = gcd(a, c);
            if (b % g != 0)
                return -1;
            b /= g;
            c /= g;
            d = (d * a / g) % c;
            cnt++;
        }
        sqt = ceil(sqrt(c));
        A = b;
        for (int i = 0; i < sqt; i++) {
            hashmap[A] = i;
            A = A * a % c;
        }
        base = qz(a, sqt, c);
        A = d * base % c;
        for (int i = 1; i <= sqt; i++) {
            if (hashmap.find(A))
                return i * sqt - hashmap[A] + cnt;
            A = A * base % c;
        }
        return -1;
    }
    int calc(int b, int n, int p) {
        b %= p, n %= p;
        if (n == 1 || p == 1)
            return 0;
        if (b == 0 && n != 0)
            return -1;
        return exbsgs(b, n, p);
    }
};

离散对数

对于模数 c,若 c 存在原根,设 g 表示 c 的原根,那么对于 gxb(modc)x[0,p)b[0,p) 不重不漏一一对应,将 gxb(modc) 的解 x 称作 g 关于 bc 的指标,也就是离散对数。

记作 x=indgb,形如对数形式。

所以求解底数是原根时的特殊指数同余方程,可以使用 bsgs 解决。

模数为质数的 N 次剩余

对于 xab(modp),pP

p 存在原根 g,设 x=gc,原式可以写成 (gc)ab(modp)(ga)cb(modp),其中 ga 是常数,那么等价于求解一个离散对数问题,且 gcd(ga,p)=1,使用 bsgs/exbsgs 求解即可。

一般意义下的 N 次剩余

对于模数非质数的 N 次剩余,利用中国剩余定理,将模数写成 m=piai 的形式,对于 pi2piai 存在原根,使用 bsgs/exbsgs 求解后使用中国剩余定理合并即可。

特别地,2x,x>1 不存在原根,不能使用 bsgs/exbsgs 求解 N 次剩余,需要使用其它算法求解。

中国剩余定理

用于求解线性同余方程组 {xr1(modm1)xr2(modm2)    xrk(modmk),其中 mi 两两互质。

M=i=1kmi,Mmi×Mmi11(modmi),ci=Mmi×(Mmi)1modM

则原同余方程组的解为 xi=1kri×ci(modM)

可以证明,线性同余方程组在模意义下解唯一,所以上述解即为原方程组的唯一解。

时间复杂度:O(klogn)

注:若不保证模数为质数,仅保证两两互质,那么可以使用欧拉定理预处理欧拉函数求逆元。或直接使用 exgcd 求解逆元。

cpp
struct CRT {
    int intchina(vector<int> &r, vector<int> &mod, int n) {
        int M = 1;
        for (int i = 1; i <= n; i++)
            M *= mod[i];
        int ans = 0;
        for (int i = 1; i <= n; i++) {
            if (mod[i] == 1)
                continue;
            int now = M / mod[i];
            int inv = qz(now % mod[i], phi[mod[i]] - 1, mod[i]);
            ans = (ans + r[i] * now * inv) % M;
        }
        return ans;
    }
};

garner 算法

a 满足如下线性方程组,且 a<i=1kpi

{aa1(modp1)aa2(modp2)     aak(modpk)

我们可以用以下形式的式子(称作 a 的混合基数表示)表示 a

a=x1+x2p1+x3p1p2+...+xkp1...pk1

garner 算法用来计算系数 x1,...,xk

ri,jpi 在模 pj 意义下的逆:

pi×ri,j1(modpj)

a 代入我们得到的第一个方程:

a1x1(modp1)

代入第二个方程得出:

a2x1+x2p1(modp2)

方程两边减 x1,除 p1 后得:

axx1x2p1(modp2)(a2x1)r1,2x2(modp2)x2(a2x1)r1,2(modp2)

类似地,我们可以得到:

xk=(...((akx1)r1,kx2)r2,k)rk1,kmodpk

时间复杂度:O(k2logp)

三模数 NTT 实现的任意模数 NTT 就是这个原理。

cpp
int Garner(vector<int> &r, vector<int> &mod, int n) {
    vector<int> x(n + 1);
    vector<vector<int>> iv(n + 1, vector<int>(n + 1));
    for (int i = 1; i <= n; i++) {
        for (int j = i + 1; j <= n; j++) {
            iv[i][j] = qz(mod[i], mod[j] - 2, mod[j]);
        }
    }
    for (int i = 1; i <= n; i++) {
        x[i] = r[i];
        for (int j = 1; j < i; j++) {
            x[i] = iv[j][i] * (x[i] - x[j]) % mod[i];
            if (x[i] < 0)
                x[i] += mod[i];
        }
    }
    int res = x[1];
    int now = 1;
    for (int i = 2; i <= n; i++) {
        now = now * mod[i - 1];
        res += now * x[i];
    }
    return res;
}

扩展中国剩余定理

用于求解线性同余方程组 {xr1(modm1)xr2(modm2)    xrk(modmk),其中 mi 不一定满足两两互质。

对于两个同余方程 xr1(modm1),xr2(modm2)

写成二元一次方程的形式:

x+b1m1=r1,x+b2m2=r2r1b1m1=r2b2m2

b1,b2 需要求解,r1,r2,m1,m2 均为常数,所以原式可以写成 m1b1m2b2=r1r2 的二元一次方程。

使用 exgcd 求解即可。

注意会出现 gcd(m1,m2)(r1r2) 时无解的情况。

将解出的 b1b2 代回原式,可得:

x=r1b1m1r1(modm1)r2(modm2)xx(modlcm(m1,m2))

至此,成功将两个同余方程合并成了一个同余方程。

依次类推,合并 k1 次,即可得到最后的同余式。

实际上此处 crt 和 excrt 讨论的都是 x 系数为 1 的同余式,其在模意义下已经是确定的常数了。

对于一般的同余方程 axb(modc),可以使用 exgcd 先求出 x 的同余式通解,再使用上述规则求解即可。

时间复杂度:O(klogn)

cpp
int eyyc(int a, int b, int c) {
    int g = gcd(a, b);
    if (c % g)
        return -1;
    int x, y;
    exgcd(a, b, x, y);
    x *= c / g;
    b /= g;
    x = (x % b + b) % b;
    return x;
}
struct ExCRT {
    int intchina(vector<int> &r, vector<int> &mod, int n) {
        int R = r[1], M = mod[1];
        for (int i = 2; i <= n; i++) {
            int x = eyyc(M, mod[i], R - r[i]); // b 取正值,保证求解得到的一定是非负值
            if (x == -1)
                return -1;
            R -= x * M;
            M = lcm(M, mod[i]);
            R %= M;
        }
        if (R < 0)
            R += M;
        return R;
    }
};

同余最短路:

同余最短路用于求解模意义下从 xy 的最短路。

例如给定 a,求 i=1naixi,xi0 能表示的最小非负整数。

若选中任一 ai,则最终的数都可以通过剩下的 n1 个数和调整 ai 的系数得到。

即可以关于 ai 划分成 ai 的同余类。

在每个同余类中求出最小的可以表示的数,即可求出最小的可以表示的数。

在模意义下,多一个 aj 相当于在 x(x+aj)modai 之间建了一条权值为 aj 的边,从 0 出发到 x[0,ai) 的最短路即为同余类中能表示的最小值。

转圈法求解同余最短路问题:

本质是:xac(modm),gcd(a,m)c

所以对于不同的 x 可以划分成 gcd(a,m) 和同余类,在每个同余类内部正向转移跑完全背包,因为是模意义下,所以要跑两轮。

转移 dp(x)+adp((x+a)modm)。从 x 开始走 k=mgcd(a,m) 步可以更新完 x 能更新的点,那么从 i 开始走 2k 步,也就是“多走一圈”,就相当于从 i 所在的剩余类中的每个点都走了 k 步,也就可以更新完这个剩余类中的每个点。

时间复杂度:O(nm)