【算法笔记】数论合集

[TOC]

线性筛

欧拉函数

定义

定义：1-n 中与 n 互质的数的个数称为欧拉函数，记为 \varphi(n)。写出 n 的唯一分解式 n = \prod p_i^{c_i}，则欧拉函数的计算式为
\varphi(n) = n\prod(1 – p_i)
证明可用容斥原理。

性质

• $\forall n \gt 1$，$[1,n]$ 中与 $n$ 互质的数之和为 $n\varphi(n)/2$。简证：与 $n$ 不互质的数 $x,n-x$ 成对出现，所以平均值为 $n/2$
• 若 a\perp b，则 \varphi(ab) = \varphi(a)\varphi(b)。显然，不证。
• 若 p\mid n 且 p^2\mid n，则 \varphi(n) = \varphi(n/p)\cdot p。简证：若 p\mid n 且 p^2\mid n，则 n,n/p 包含相同质因子，写出欧拉函数计算式可得结论
• 若 p\mid n 但 p^2 \not\mid n，则 \varphi(n) = \varphi(n/p)(p-1)。简证：说明 p\perp n/p，易知 \varphi(n) =\varphi(n/p)\varphi(p)=\varphi(n/p)(p-1)
• $\displaystyle\sum_{d\mid n}\varphi(d) = n$。简证：设 $f(n) = \displaystyle\sum_{d\mid n}\varphi(d)$。由 $\varphi$ 为积性函数，得到 $f(nm) = \sum_{d\mid nm}\varphi(d) = (\sum_{d\mid n}\varphi(d))(\sum_{d\mid m}\varphi(d)) = f(n)f(m)$。所以 $f$ 为积性函数。对于单个质因子 $f(p^m) = \varphi(1)+\varphi(p)+\varphi(p^2)+\cdots+\varphi(p^m) = p^m$。所以 $f(n) = n$

积性函数：若 a\perp b 且 f(ab) = f(a)f(b)，则 f 为积性函数。

同余类与剩余系

对于 \forall a\in [0, m – 1]，有集合 \lbrace a + km\rbrace, k\in\mathbb{Z} 中的所有元素模 m 同余，称这个集合为一个模 m 的同余类，记为 \overline{a}。模 m 的同余类有 m-1 个，构成 m 的完全剩余系。[1, m] 中与 m 互质的数有 \varphi(m) 个，构成 m 的简化剩余系。

简化剩余系关于模 m 乘法封闭。因为若 a\perp m, b\perp m，则一定有 ab\perp m。

费马小定理与欧拉定理

费马小定理

若 p 为质数，则对于 \forall a\in\mathbb{Z}，有
a^p\equiv a\pmod p
或者记为
a^{p-1} \equiv 1\pmod p
上式常用于求乘法逆元，要求 a\perp p

欧拉定理

若 a,b\in\mathbb{Z}^+,a\perp b，则
a^{\varphi(n)}\equiv 1\pmod n

证明

设 n 的简化剩余系为 \lbrace\overline{a_1},\overline{a_2},\cdots,\overline{a_{\varphi(n)}} \rbrace。对 \forall a_i, a_j，若 aa_i\equiv aa_j\pmod n，则 a(a_i-a_j)\equiv 0\pmod n。因为 a\perp n，所以 a_i\equiv a_j。故当 a_i\not=a_j 时，aa_i,aa_j 也代表不同同余类。

而简化剩余系关于模 n 乘法运算封闭 ，所以简化剩余系 \lbrace\overline{a_1},\overline{a_2},\cdots,\overline{a_{\varphi(n)}} \rbrace 与 \lbrace\overline{aa_1},\overline{aa_2},\cdots,\overline{aa_{\varphi(n)}} \rbrace 是等价的。综上：
a^{\varphi(n)}a_1a_2\cdots a_{\varphi(n)}\equiv (aa_1)(aa_2)\cdots(aa_{\varphi(n)})\equiv a_1a_2\cdots a_{\varphi(n)}\equiv 1\pmod n
欧拉定理得证。

当 p 为质数时，\varphi(p) = p – 1，且只有 p 的倍数与 p 不互质。显然当 a\perp p 时 a^{p – 1} \equiv 1\pmod n，两边同乘 a 得到费马小定理。而当 a\not\perp p 时费马小定理显然成立。

Q.E.D.

推论

若 a\perp n，则 \forall b\in\mathbb{Z}^+，有 a^b\equiv a^{b\bmod \varphi(n)}\pmod n

证明：

设 b = q\varphi(n) + r，即 r = b\bmod{\varphi(n)}。所以 a^b\equiv a^{q\varphi(n) + r}\equiv (a^{\varphi(n)})^qa^r\equiv a^r\pmod n

拓展欧几里得算法

回顾欧几里得算法的过程。

int gcd(int a, int b)
{
    if (!b)
        return a;
    return gcd(b, a % b);
}

这是一个递归的过程。思考：我们还可以利用这个算法做些什么。

拓展欧几里得算法可以求解形如 ax + by = \gcd(a, b) 的一类不定方程。具体的

乘法逆元

定义

若

ax\equiv 1\pmod p

且

a\perp p (\gcd(a,p)=1)

则称 x 是 a 在模 p 意义下的逆元，可表示为 a^{-1}。

费马小定理

假如 a 是一个整数，p 是一个质数，则 a^p-a 是 p 的倍数，表示为

a^p\equiv a\pmod p

或

a^{p-1}\equiv 1\pmod p

求法一：费马小定理

由 a^{p-1}\equiv a^{p-2}\cdot a\equiv1\pmod p

即当 p 是质数时，可以直接快速幂求 a^{p-2}\bmod p 即为所求逆元

求法二：拓欧

方程 ax\equiv 1\pmod b 可以转化为 ax+by=1（其中 y 为整数），若保证有解的话说明 \gcd(a,b)=1，我们就可以使用拓展欧几里得算法进行求解。

若题目要求输出最小正整数解，可以将求出来的 x 先加上 b 后对 b 取模。

#include <cstdio>

int a,b;

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

int main()
{
    scanf("%d %d",&a,&b);
    int x,y;
    exgcd(x,y,a,b);
    printf("%d\n",(x+b)%b);
    return 0;
}

求法三：线性求逆元

令 p=ki+r，其中 k=\lfloor \frac p i\rfloor，r=p\bmod i，且 1<r<i<p。则在模 p 意义下有

ki+r\equiv 0\pmod p

两边同时乘以 r^{-1}i^{-1}，则

kr^{-1}+i^{-1}\equiv0\pmod p

移项，

i^{-1}\equiv-kr^{-1}\pmod p

代入 k=\lfloor \frac p i\rfloor，r=p\bmod i，有

i^{-1}\equiv -\lfloor \frac p i\rfloor(p\bmod i)^{-1}\pmod p

由于要保证 i^{-1}>0，在最终式子的右边加上 p（p\equiv 0\pmod p），最终的式子就是：

i^{-1}\equiv p-\lfloor \frac p i\rfloor(p\bmod i)^{-1}\pmod p

令 inv[i] 表示 i^{-1} 则递推式如下：

for(int i=2;i<=n;i++)
    inv[i]=(p-p/i)*inv[p%i]%p;

代码（P3811）：

#include <cstdio>

typedef long long ll;
const int maxn=3e6+5;
ll n,p;
ll inv[maxn]={0,1};

int main()
{
    scanf("%lld %lld",&n,&p);
    for(ll i=2;i<=n;i++)
        inv[i]=(p-p/i)*inv[p%i]%p;
    for(int i=1;i<=n;i++)
        printf("%lld\n",inv[i]);
    return 0;
}

中国剩余定理

是什么

中国剩余定理（孙子定理，Chinese remainder theorem，CRT）是中国古代求解一次线性同余方程组的定理。

怎么用

CRT 可以求解如下的一次线性同余方程组

\begin{cases}
x &\equiv a_1\pmod{m_1} \
x &\equiv a_2\pmod{m_2} \
&\vdots \
x&\equiv a_n\pmod{m_n}
\end{cases}

其中 m_i，m_j（i,j\in[1,n] 且 i\neq j）要求两两互质。

流程如下：

1. 计算 \displaystyle M=\prod_{i=1}^nm_i，\displaystyle M_i=\frac M {m_i}
2. 对于每个 M_i，计算其在模 m_i 意义下的逆元 t_i，即 M_it_i\equiv1\pmod{m_i}
3. 方程组的一个特解 x_0=\displaystyle\sum_{i=1}^na_iM_it_i
4. 最小正整数解即为 (x_0+M)\bmod M，方程组的解集为 \left\lbrace x|x=\displaystyle\sum_{i=1}^na_iM_it_i+kM,k\in\mathbb{Z}\right\rbrace

为什么

下面给出 CRT 的证明：

首先，我们易知 \forall j\in[1,n]，当 i\neq j 时有

a_jM_jt_j\equiv 0\pmod{m_i}

当 i=j 时有

a_iM_it_i\equiv a_iM_iM_i^{-1}\equiv a_i\pmod{m_i}

所以对于我们得到的解 x=\displaystyle\sum_{i=1}^na_iM_it_i，对于任意的 m_i 都有 x\equiv a_i\pmod{m_i}，定理得证。

怎么实现

模板题 https://www.luogu.com.cn/problem/P1495

代码实现如下：

#include <cstdio>

typedef long long ll;
ll n,a[15],m[15],Mul,M[16],inv[16],x,y;

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

int main()
{
    Mul=1;
    scanf("%lld",&n);
    for(int i=1;i<=n;i++)
    {
        scanf("%lld %lld",m+i,a+i);
        Mul*=m[i];
    }
    ll ans=0;
    for(int i=1;i<=n;i++)
    {
        M[i]=Mul/m[i];
        x=0,y=0;
        exgcd(M[i],m[i],x,y);
        inv[i]=x<0?x+m[i]:x;
        ans+=a[i]*M[i]*inv[i];
    }
    printf("%lld\n",ans%Mul);
    return 0;
}

拓展中国剩余定理

问题引入

当我们需要解形如

\begin{cases}
x\equiv r_1\pmod{m_1}\
x\equiv r_2\pmod{m_2}\
\cdots\
x\equiv r_n\pmod{m_n}
\end{cases}

线性同余方程组时，我们也许会想到 CRT，但是回忆一下 CRT 的过程，其中有一步是这样的：

对于每个 M_i，计算其在模 m_i 意义下的逆元 t_i

回忆一下 M_i：\displaystyle M_i=\prod_{j = 1,i\ne j,}^n m_j。要满足其在模 m_i 意义下的逆元存在，首先就需要满足所有的模数 m_i 两两互质。然而这里的例子中模数都不一定两两互质，怎么办呢？

转化问题

这时就不能再局限于 CRT 的构造解的思维了，先从只有两个方程的情况考虑起：

\begin{cases}
x\equiv r_1\pmod{m_1} \
x\equiv r_2\pmod{m_2} \
\end{cases}

这个方程可以改写为：

x=k_1m_1+r_1=k_2m_2+r_2

其中 k_1,k_2 为未知量，然后方程可以进一步改写为关于 k_1,k_2 的形式：

k_1m_1 – k_2m_2=r_2-r_1

由 Bézout 定理得这个方程有解的条件是 \gcd(m_1,m_2)\mid r_2-r_1，若有解， 则令 \displaystyle d = \gcd(m_1,m_2),p_1 = \frac{m_1}{d},p_2=\frac{m_2}{d}，方程化为

k_1p_1-k_2p_2 = \frac{r_2-r_1}{d}

使用 exgcd 求出方程 p_1\lambda_1 + p_2\lambda_2 = 1 的解 \lambda_1,\lambda_2，则

\begin{cases}
k_1=\lambda_1\cdot\dfrac{r_2-r_1}{d} \
k_2=\lambda_2\cdot\dfrac{r_2-r_1}{d}
\end{cases}

把 k_1 还原回去，得到满足原方程的一个特解

x^* = m_1\lambda_1\cdot\dfrac{r_2-r_1}{d} + r_1

则这个方程所有的通解为 x = x^* + k\cdot\operatorname{lcm}(m_1,m_2),k\in\mathbb{Z}。

所以两个方程就合并为一个：

x\equiv x^*\pmod{\operatorname{lcm}(m_1,m_2)}

这样子两两合并下去就可以将所有的方程合为一个，得解。

实现

注意乘法过程可能溢出，需要龟速乘， 然后就是负数取模的问题

#include <cstdio>
#include <cctype>
#define il inline

typedef long long ll;

inline ll read()
{
    char c = getchar();
    ll s = 0;
    bool x = 0;
    while (!isdigit(c))
        x = x | (c == '-'), c = getchar();
    while (isdigit(c))
        s = 10 * s + c - '0', c = getchar();
    return x ? -s : s;
}

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

ll mul(ll a, ll b, ll mod)
{
    ll ans = 0, x = a;
    b = (b + mod) % mod, a = (a + mod) % mod;
    for (; b; b >>= 1ll)
    {
        if (b & 1ll)
            ans = (ans + x) % mod;
        x = (x << 1) % mod;
    }
    return ans;
}

const int maxn = 1e5 + 5;
ll n, m[maxn], r[maxn];

int main()
{
    n = read();
    for (int i = 1; i <= n; ++i)
        m[i]= read(), r[i] = read();
    for (int i = 2; i <= n; ++i)
    {
        ll tmpx, tmpy;
        ll gcd = exgcd(m[i], m[i - 1], tmpx, tmpy);
        ll p1 = m[i - 1] / gcd, p2 = m[i] / gcd;
        exgcd(p1, p2, tmpx, tmpy);
        m[i] = m[i] / gcd * m[i - 1];
        r[i] = (r[i - 1] + mul(m[i - 1], mul((r[i] - r[i - 1]) / gcd, tmpx, m[i]), m[i])) % m[i];
    }
    printf("%lld\n", r[n] % m[n]);
    return 0;
}

BSGS

问题引入

给定 a,b,p\in\mathbb{Z}，其中 a\perp p，求一个非负整数 x 使得

a^x\equiv b\pmod p

该问题也被称为离散对数。

暴力算法

首先 a^0\equiv 1\pmod p，由费马小定理，a^{p – 1}\equiv1\pmod p。所以只需要在 [1, p – 1] 之间枚举答案。

优化

设 x = i\sqrt p – j，则有

a^{i\sqrt p – j}\equiv b\pmod p

两边同时除以 b

\frac{a^{i\sqrt p}}{a^jb}\equiv 1\pmod p

考虑预处理分母然后枚举 i，所以如果在模 p 意义下 a^{i\sqrt p} = a^jb，那么方程有解 i\sqrt p – j。

同时由于解的形式是 i\sqrt p – j，所以 j\in[0, \sqrt p)。预处理出 a^jb\bmod p 在哈希表（可用 unordered_map 代替）中即可。

#include <cstdio>
#include <cmath>
#include <unordered_map>
#define R register
#define FOR(i, a, b) for(R int i = a; i <= b; ++i)
#define int long long

std::unordered_map<int, int> hash;

int ksm(int base, int p, int mod)
{
    int ans = 1;
    for (; p; p >>= 1)
    {
        if (p & 1)
            ans = ans * base % mod;
        base = base * base % mod;
    }
    return ans;
}

int BSGS(int a, int b, int p)
{
    hash.clear();
    if (a % p == 0)
        return -1;
    int t = ceil(sqrt(p)), val, tmp;
    hash[val = b] = 0;
    FOR(j, 1, t)
        hash[val = val * a % p] = j;
    a = ksm(a, t, p), val = 1;
    FOR(i, 1, t)
    {
        val = val * a % p;
        if (tmp = hash[val])
            return val = i * t - tmp, (val % p + p) % p;
    }
    return -1;
}

signed main()
{
    int p, a, b;
    scanf("%lld %lld %lld", &p, &a, &b);
    int ans = BSGS(a, b, p);
    if (ans >= 0)
        return printf("%lld\n", ans), 0;
    else return printf("no solution\n"), 0;
} 

拓展

还是求解

a^x\equiv b\pmod p

当 a\not\perp p 时，原来的法子行不通了因为 a^jb 关于模 p 的逆元不存在了，所以考虑让他们变得互质。原方程可化为

a^x + kp = b

两边同时除以 g = \gcd(a, p)，得到

a^{x-1}\frac{a}{d} + k\frac p d = \frac b d

显然当 g\not\mid b 时方程无解。所以原方程就变为了

\frac a d a^{x-1} \equiv \frac b d\pmod{\frac p d}

记所有的 \displaystyle\frac{b}{\prod d} = B, \frac{p}{\prod d} = P,\prod\frac a d = A，继续套娃递归至最后 a\perp \frac p D 即可，记一共套了 k 次娃，则方程变为了

Aa^{x – k}\equiv B\pmod P

此时 A\perp P，所以我们就可以愉快地求出 A^{-1}\pmod P 然后把他乘到右边来，方程就化为

a^{x – k}\equiv A^{-1}B\pmod P

用 BSGS 求出该方程的解即可。

细节：如果 a = 1\or p = 1，最小自然数解必然为 0，然后如果在递归的过程中发现 A=B，则说明 x = k，可直接返回答案。

代码：

#include <cstdio>
#include <cmath>
#include <unordered_map>
#define R register
#define FOR(i, a, b) for(R int i = a; i <= b; ++i)
#define int long long

std::unordered_map<int, int> hash;

int ksm(int base, int p, int mod)
{
    int ans = 1;
    for (; p; p >>= 1)
    {
        if (p & 1)
            ans = ans * base % mod;
        base = base * base % mod;
    }
    return ans;
}

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

int inv(int a, int mod)
{
    int x, y;
    exgcd(a, mod, x, y);
    return (x % mod + mod) % mod;
}

int BSGS(int a, int b, int p)
{
    hash.clear();
    if (a % p == 0)
        return -1;
    int t = ceil(sqrt(p)), val, tmp;
    hash[val = b] = 0;
    FOR(j, 1, t)
        hash[val = val * a % p] = j;
    a = ksm(a, t, p), val = 1;
    FOR(i, 1, t)
    {
        val = val * a % p;
        if (tmp = hash[val])
            return val = i * t - tmp, (val % p + p) % p;
    }
    return -1;
}

int exBSGS(int a, int b, int p)
{
    if (a == 1 || p == 1)
        return 0;
    int x, y, k = 0, na = 1, g = exgcd(a, p, x, y);
    while (g > 1)
    {
        if (b % g)
            return -1;//判无解
        ++k, b /= g, p /= g, na = (na * (a / g)) % p;//递归
        if (na == b)//特判可以直接返回答案的情况
            return k;
        g = exgcd(a, p, x, y);
    }
    int ans = BSGS(a, b * inv(na, p) % p, p);
    return ans == -1 ? -1 : ans + k;
}

signed main()
{
    int a, p, b;
    for (scanf("%lld %lld %lld", &a, &p, &b); a || b || p; scanf("%lld %lld %lld", &a, &p, &b))
    {
        int ans = exBSGS(a % p, b % p, p);
        if (ans > 0)
            printf("%lld\n", ans);
        else printf("No Solution\n");
    }
    return 0;
}