算法基础（11）

欧拉函数

1 ~ N 中与 N 互质的数的个数被称为欧拉函数，记为$\phi(N)$。若在算数基本定理中，$N=p_1^{a_1}p_2^{a_2}\dots p_k^{a_k}$（分解质因数） ，则：

$$\phi(N)=N \cdot (1-\frac{1}{p_1}) \cdot (1-\frac{1}{p_2}) \dots (1-\frac{1}{p_k})$$

互质是公约数只有1的两个整数。

证明利用了容斥原理：

1. 从1到N中去掉$p_1,p_2,\dots,p_k$的所有倍数；

2. 有些数会被减掉两次，需要再加回来，即加上所有$p_i * p_j$的倍数；

3. 若某些数是$p_1, p_2, p_3$的公倍数，那么它在第一步被减去三次，在第二步被加上三次，相当于是没处理掉，因此要再减去所有$p_ip_jp_k$的倍数；
4. 按此规律继续下去，加上所有四个质因数的倍数；减去所有5个质因数的倍数；……

$$\begin{align*} N &-\frac{N}{p_1}-\frac{N}{p_2}\dots-\frac{N}{p_k} \\ &+\frac{N}{p_1p_2}-\frac{N}{p_1p_3}\dots-\frac{N}{p_{k-1}p_k} \\ &+\frac{N}{p_1p_2p_3}-\frac{N}{p_1p_2p_4}\dots-\frac{N}{p_{k-2}p_{k-1}p_k} \\ &-\frac{N}{p_1p_2p_3p_4}-\dots\\ &\dots \end{align*}$$

而上面两式是相等的，即证。

用上面的公式计算欧拉数的时间复杂度为$O(\sqrt{n})$，瓶颈在分解质因数，而分解质因数的时间复杂度为$O(\sqrt{n})$。

给定n个正整数$a_i$，请你求出每个数的欧拉函数。

数据范围：$1 \le n \le 100, 1 \le a_i \le 2∗10^9$

#include<iostream>
using namespace std;

int main()
{
    int n;
    cin >> n;
    
    while(n--)
    {
        int a;
        cin >> a;
        
        int res = a;   
        for(int i = 2; i <= a / i; i++)
            if(a % i == 0)
            {
                res = res / i * (i -1);   //套用欧拉函数的公式，注意这里要先除再乘，防止int溢出
                while(a % i == 0)
                    a /= i;
            }
        if(a > 1)  res = res / a * (a - 1);
        
        cout << res << endl;
    }
    return 0;
}

• 用筛法求欧拉函数

给定一个正整数n，求1~n中每个数的欧拉函数之和。

数据范围：$1 \le n \le 10^6$。

若是用上面的公式就1~n中每个数的欧拉函数，那么时间复杂度就是$O(n\sqrt{n})$，若是借用之前讲的筛质数的线性筛的思路，可以将时间复杂度优化到$O(n)$。

#include<iostream>
using namespace std;

const int N = 1000010;

typedef long long LL;

int primes[N], cnt;
int phi[N];
bool st[N];

LL get_euler(int n)
{
    phi[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!st[i])
        {
            primes[cnt ++] = i;
            phi[i] = i - 1;      //i是质数，i的互质数（从1到i之中的）的个数是i - 1
        }
        for(int j = 0; primes[j] <= n / i; j++)
        {
            st[primes[j] * i] = true;      //primes[j]一定是primes[j] * i的最小质因数
            if(i % primes[j] == 0)  //如果primes[j]是i的最小质因数
            {
                phi[primes[j] * i] = phi[i] * primes[j];
                break;
            }
            phi[primes[j] * i] = phi[i] * (primes[j] - 1);
        }
    }
    
    LL res = 0;
    for(int i = 1; i <= n; i++)  res += phi[i];
    
    return res;
}
int main()
{
    int n;
    cin >> n;
    
    cout << get_euler(n) << endl;
    return 0;
}

在线性筛法的代码上作添加即可，用phi[i]存数i的欧拉函数（从1到i的互质数的个数）：

1. 如果st[i]=false，即i是质数，那么i的欧拉函数就是i-1；
2. 如果$i \, mod \, p_j = 0$，则$p_j$是$i$的最小质因数，也是$p_ji$的最小质因数；而$p_ji$的分解质因数的结果只比$i$分解质因数的结果多了一项：$p_j$，又因为$p_j$是$i$的质因数，因此在$\phi(i)$的公式中已经计算过了$(1-\frac{1}{p_j})$这一项，那么有$\phi(p_ji)= \phi(i) p_j$；
3. 如果$i \, mod \, p_j \ne 0$，$p_j$是$p_ji$的最小质因数，但不是$i$的最小质因数，若设$\phi(i)= i \cdot (1-\frac{1}{p_1}) \cdot (1-\frac{1}{p_2}) \dots (1-\frac{1}{p_k})$，则$\phi(p_j i)= p_j i \cdot (1-\frac{1}{p_1}) \cdot (1-\frac{1}{p_2}) \dots (1-\frac{1}{p_k}) (1-\frac{1}{p_j})$，那么有$\phi(p_ji)=p_j \phi(i) (1-\frac{1}{p_j})=\phi(i) * (p_j-1)$。

欧拉函数的一个用处——欧拉定理：若$a$与$n$互质，则有$a^{\phi(n)}\equiv 1(mod \, n)$。

同余：两个整数$a, b$，若它们除以正整数$m$所得的余数相等，则称$a, b$对于模$m$同余，记作$a \equiv b( \mod m)$。读作$a$同余于$b$模$m$，或读作$a$与$b$关于模$m$同余。

快速幂

给定n组$a_i,b_i,p_i$，对于每组数据，求出 $a_i^{b_i} \mod p_i$ 的值。

数据范围：$1 \le n \le 100000, 1 \le a_i,b_i,p_i \le 2∗10^9$

快速幂可以在$O(\log n)$的时间复杂度，求出$a^k \, mod \,p$的值，其中$1 \le a, p, k \le 10^9$。

思路是：预处理出这些值，$a^{2^0} \, mod\, p$，$a^{2^1} \, mod\, p$，$a^{2^2} \, mod\, p$，……，$a^{2^{\log k}} \, mod\, p$，然后让：

$$a^k=a^{2^{x_1}} \cdot a^{2^{x_2}} \dots a^{2^{x_t}} =a^{2^{x_1}+2^{x_2}+\dots +2^{x_t}}$$

而具体的$x_1, x_2, \dots, x_t$则可以由$k$的二进制表示所有为1的位获得，如$(k)_{10}=(110110)_2$， 则有$k=2^1+2^2+2^4+2^5$。

#include<iostream>
using namespace std;

typedef long long LL;   //数论中的很多问题都会爆int， 会用到long long 

//返回 a^k mod p 的结果
int qmi(int a, int k, int p)
{
    int res = 1;
    while(k)
    {
        if(k & 1)  res = (LL)res * a % p;
        a = (LL)a * a % p;
        k >>= 1;
    }
    return res;
}

int main()
{
    int n;
    scanf("%d", &n);
    
    while(n--)
    {
        int a, k, p;
        scanf("%d%d%d", &a, &k, &p);
        
        printf("%d\n", qmi(a, k, p));
    }
    
    return 0;
}

• 快速幂求逆元

乘法逆元的定义：若整数$b，m$互质，并且对于任意的整数$ a$，如果满足$b|a$，则存在一个整数$x$，使得$a/b≡a∗x(mod \,m)$，则称$x$为$b$的模$m$乘法逆元，记为$b^{−1}(mod \,m)$。$b$存在乘法逆元的充要条件是$b$与模数$m$互质。当模数$m$为质数时，$b^{m−2}$即为$b$的乘法逆元。

$b \cdot x \equiv 1(mod \, p)$，由费马小定理：$b^{p-1} \equiv 1(mod \, p)$，即$b \cdot b^{p-2} \equiv 1(mod \, p)$，因此我们要求的逆元$x$就是$b^{p-2}$，当然前提是$b$与模数$p$互质，其$p$是质数，这就转换了求快速幂的问题，qmi(a, p - 2, p)

给定$n$组$a_i,p_i$，其中$p_i$是质数,求$a_i$模$p_i$的乘法逆元，若逆元不存在则输出impossible。注意：请返回在0∼p−1之间的逆元。

数据范围：$1 \le n \le 10^5,1 \le a_i,p_i \le 2∗10^9$。

#include<iostream>
using namespace std;

typedef long long LL;

//返回 a^k mod p 的结果
int qmi(int a, int k, int p)
{
    int res = 1;
    while(k)
    {
        if(k & 1)  res = (LL)res * a % p;
        a = (LL)a * a % p;
        k >>= 1;
    }
    return res;
}

int main()
{
    int n;
    scanf("%d", &n);
    
    while(n--)
    {
        int a, p;
        scanf("%d%d", &a, &p);
        
        int res = qmi(a, p - 2, p);
        if(a % p)  printf("%d\n", res);
        else puts("impossible");
    }
    
    return 0;
}

扩展欧几里得算法

裴蜀定理：对于任意正整数$a, b$，一定存在非零整数$x, y$，使得$ax+by=(a, b)$。

$(a, b)$表示$a,b$的最大公约数

要证明存在可以使用构造法，扩展欧几里得算法就提供了一种构造的思路。

给定$n$对正整数$a_i,b_i$，对于每对数，求出一组$x_i,y_i$，使其满足$a_i∗x_i+b_i∗y_i=gcd(a_i,b_i)$。

数据范围：$1 \le n \le 10^5 , 1 \le a_i,b_i \le 2∗10^9$

#include <iostream>
using namespace std;

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

int main()
{
    int n;
    scanf("%d", &n);

    while (n -- )
    {
        int a, b, x, y;
        scanf("%d%d", &a, &b);

        exgcd(a, b, x, y);
        printf("%d %d\n", x, y);
    }

    return 0;
}

设$a x+by=d $，由欧几里得算法知，$(a, b)=(b, a \,mod \, b)$，则有：$by+(a \, mod \, b)=d$；又因为$(a \, mod\, b)=a-\lfloor\frac{a}{b} \rfloor \cdot b$，代入得：$ax+b(y-\lfloor \frac{a}{b} \rfloor \cdot x)=d$，因此在递归d=exgcd(b, a % b, y, x)后要令y -= a \b * x。

扩展欧几里得算法的一个应用：求解线性同余方程

给定$n$组数据$a_i,b_i,m_i$，对于每组数求出一个$x_i$，使其满足$a_i∗x_i≡b_i(mod\, m_i)$，如果无解则输出impossible。

数据范围：$1 \le n \le 10^5 , 1 \le a_i,b_i, m_i \le 2∗10^9$

若存在一个$x$，使得$ax \equiv b(mod\, m)$，即是存在一个$y$，使得$ax=my+b$，即$ax-my=b$，令$y’=-y$，则等价于方程$ax+my’=b$有解，这就是扩展欧几里得算法的形式了，上式有解的充分必要条件是$b$是$a$和$m$的最大公约数的倍数，即$(a, m)|b$。

#include <iostream>
using namespace std;

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

int main()
{
    int n;
    scanf("%d", &n);

    while (n -- )
    {
        int a, b, m;
        scanf("%d%d%d", &a, &b, &m);
        
        int x, y;
        int d = exgcd(a, m, x, y);
        if(b % d)  puts("impossible");   //如果b不是gcd(a, b)的倍数，那一定无解
        else printf("%d\n", (long long) x * (b / d) % m);
    }

    return 0;
}

中国剩余定理

可参考：中国剩余定理

中国剩余定理 (Chinese Remainder Theorem, CRT) 可求解如下形式的一元线性同余方程组（其中$m_1, m_2, \dots, m_k$ 两两互质）：

$$\begin{cases} x &\equiv a_1(\mod m_1) \\ x &\equiv a_2(\mod m_2) \\ &\dots\\ x &\equiv a_k(\mod m_k) \\ \end{cases}$$

求解步骤：

1. 求所有模数的积， $M=m_1m_2\dotsm_k$
2. 对于第$i$个方程：$M_i=\frac{M}{m_i}$，求$M_i$在模$m_i$的逆元$M_i^{-1}$（可以扩展欧几里得算法解，令$b=1$）
3. 方程组的唯一解：$a= \sum^k_{i=1} a_iM_iM_i^{-1}(\mod n)$

表达整数的奇怪方式：

给定 $2n$ 个整数$a_1,a_2,…,a_n$和$m_1,m_2,…,m_n$,求一个最小的非负整数 $x$，满足$∀i∈[1,n],x≡m_i(\mod a_i)$。

数据范围：$1 \le a_i \le 2^{31}−1 , 0 \le m_i \le a_i, 1 \le n \le 25$。

注意本题中的$a_i, m_i$并没有任何限制，而中国剩余定理中要求$m_1, m_2, \dots,m_k$两两互斥。先来分析前两个式子：

$$\begin{align*} x &= k_1 * a_1 + m_1 \\ x &= k_2 * a_2 + m_2 \\ \end{align*}$$

因此有：

$$k_1*a_1-k_2*a_2=m_2-m_1$$

由扩展欧几里得算法可知，上式有解，等价于$a_1, a_2$的最大公约数能整除$m_2-m_1$，是$(a_1, a_2)|m_2-m_1$。不定方程的所有解为（假设已求出一组$k_1, k_2$，其中$k$是任意整数）：

$$\begin{cases} k_1+k \frac{a_1}{d}\\ k_2+k\frac{a_2}{d} \\ \end{cases}$$

则$x$的所有解为：

$$\begin{align*} x&=k_1a_1+m_1\\ &=(k_1+k \frac{a_1}{d})a_1+m_1\\ &=a_1k_1+m_1+k \frac{a_1a_2}{d} (前两项记为m_0，后一项记为ka)\\ &=m_0+ka \end{align*}$$

发现前两个方程$x$的解的形式一致，因此通过这个方法可以把两个不定方程合并为一个，若有$n$个不定方程，即通过$n-1$合并可以转化为一个方程$x=m_0+ka$，即$x \, mod\, a \equiv m_0$，即求$m_0 \, mod \, a$的正余数。

#include<iostream>
#include<algorithm>

using namespace std;
typedef long long LL;

LL exgcd(LL a, LL b, LL &x, LL &y)          //扩展欧几里得算法
{
    if(!b)
    {
        x = 1, y = 0;
        return a;
    }
    LL d =exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

int main()
{
    int n;
    cin >> n;
    
    bool has_answer = true;
    LL a1, m1;
    
    cin >> a1 >> m1;
    
    for(int i = 0; i < n - 1; i ++)
    {
        int a2, m2;
        cin >> a2 >> m2;
        
        LL k1, k2;
        LL d = exgcd(a1, a2, k1, k2);   //求a1, a2最大公约数，这时已经计算出来k1, k2的解
        if((m2 - m1) % d)    //如果a1, a2最大公约数不能整数m2 - m1，则无解
        {
            has_answer = false;
            break;
        }
        
        k1 *= (m2 - m1) / d;
        //将k1变成方程的最小整数解，防止溢出
        LL t = a2 / d;
        k1 = (k1 % t + t) % t;     //求k1 模 t 的正的余数
        
        m1 = a1 * k1 + m1;         //求合并方程的m
        a1 = abs(a1 / d * a2);     //求合并方程的a， 即是a1, a2的最小公倍数
    }
    
    if(has_answer)
        cout << (m1 % a1 + a1) % a1 << endl;   //如有解，则输出
    else 
        puts("-1");
        
    return 0;
}