acwing数学知识（二）欧拉函数 欧拉定理 快速幂 扩展欧几里得算法 中国剩余定理

1.欧拉函数

1 ~ N 中与 N 互质的数的个数被称为欧拉函数

欧拉函数的证明

利用容斥原理，求1 ~N-1中与N互斥的数的个数s
拆分出N的质因子p1、p2、p3…
s = N - N/p₁ - N/p₂ - N/p₃-…-N/p_k
+N/(p₁*p₂) + N/(p₁*p₂) + … + N/(p₁*p_k) + … + N/(p_k-1*p_k)
- N/(p₁*p₂*p₃) - N/(p₁*p₂*p₄) - … -N/(p_k-2*p_k-1*p_k)
+N/(p₁*p₂*p₃*p₄) + … + N/(p_p-3*p_p-2*p_p-1*p_p)
…（以此类推）

1. p1的倍数由N/p1个，p2的倍数有N/p2个…p_k的倍数有N/p_k个，因此要减去。
2. 但是有的数同时是p1和p2的倍数，所以会被重复减去，因此在第二行加回。
3. 但是有的数同时是p1、p2、p3的倍数，在第一行会被重复减去三次，而在第二行会重复加回三次，相当于没变，因此要在第三行减去。
4. 以此类推。。。。。。。

最后欧拉函数的公式如下：s = N * (1-1/p₁) * (1-1/p₂) * (1-1/p₃) … * (1-1/p_k)

接下来用公式来求欧拉函数 O(n*sqrt(n))

#include <iostream>
#include <algorithm>

using namespace std;

typedef long long ll;//数据范围较大记得用long long
const int N = 110;

int n;

int main()
{
    cin >> n;
    while(n--)
    {
        int cnt = 0;
        ll primes[N];
        int x;
        cin >> x;
        ll ans = x;
        
        for(int i = 2 ; i <= x / i ; i++)
        {
            if(x % i == 0)
            {
                ans = ans / i * (i - 1);
                while(x % i == 0) x /= i;
            }
        }
        
        if(x > 1) ans = ans / x * (x - 1);

         
        cout << ans << endl;
    }
    return 0;
}

线性筛法求欧拉函数之和O(n)
核心：在找到质数和筛掉和筛掉合数时，利用欧拉函数的公式直接计得出其欧拉函数

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 1000010;


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


void get_eulers(int n)
{
    euler[1] = 1;
    for (int i = 2; i <= n; i ++ )
    {
        if (!st[i])
        {
            primes[cnt ++ ] = i;
            euler[i] = i - 1;
			//如果一个数a是质数，那么它的欧拉函数是a-1
        }
        for (int j = 0; primes[j] <= n / i; j ++ )
        {
            int t = primes[j] * i;
            st[t] = true;
            if (i % primes[j] == 0)
   //当pj是i的质因子时，pj的质因子都在i的质因子中，因此pj*i的质因子与i的质因子相
   //同，而一个数的欧拉函数与其质子的次数无关，因此推导公式得：φ（pj*i)=pj*φ(i)
            {
                euler[t] = euler[i] * primes[j];
                break;
            }
  //当pj不是i的质因子时，pj*i的质因子只比i的质因子多一个pj，因此φ(pj*i)的展开式在
  //φ(i)的展开式中多乘上一项(1-1/pj)和pj即可，推导化简得：φ（pj*i)=(pj-1)φ(i)          
            euler[t] = euler[i] * (primes[j] - 1);
        }
    }
}


int main()
{
    int n;
    cin >> n;

    get_eulers(n);

    LL res = 0; //最后的答案要开long long 数据是1e6
    for (int i = 1; i <= n; i ++ ) res += euler[i];

    cout << res << endl;

    return 0;
}

2.欧拉定理

定理：若a与n互质，则a^Φ(n) ≡ 1(mod n)

证明：
假如1 ~ n中，所有与n互质的数是A数组：a₁、a₂、a₃…a_Φ(n)
设M数组：aa₁、aa₂、aa₃…aa_Φ(n)(同样与n互质)

下面先证明两个推理：
一、M数列中不存在两个数模n同余。
用假设法证明：
1.假如M数列存在m_a与m_b模n同余，即m_a≡m_b(mod n)
2.移项得m_i - m_j ≡ 0 (mod n)
3.即aa_i-aa_j ≡ 0 (mod n)
4.提取公因式得:a(a_i-a_j) ≡ 0 (mod n)
5.因为a与n互质，所以要是等式左边与n互质，那么必然:(a_i-a_j) ≡ 0 (mod n)
6.则a_i ≡ a_j(mod n),这显然与A数组中的数的性质相违背（A数组中各个数都小于n，并且不同，因此不可能模n同余）。
7.所以假设不成立！即：M数列中不存在两个数模n同余！

二、M中的数对n的余数全部与n互质（即m_i%n 与 n 互质）
证明：
1.已知a与n互质，a_i与n互质，则aa_i也与n互质，也可得M数列中的数与n互质。
2.然后带入欧几里得算法中推导即可：gcd(aa_i , n) = gcd(n , aa_i%n ) = 1 ， 即m_i%n 与 n 的最大公约数是1
（a与n互质，p_i与n互质，因此ap_i与n互质）
3.得证：M中的数对n的模全部与n互质！

根据这两个推理，可以开始推导欧拉定理。

根据推理二，我们可以知道M中的数模n全部与n互质，而且M中不

存在两两对n取模相等的数，而且A、M数组长度相同，因此一个M中的数一定

可以在A中找到唯一一个数和它模n同余（即两个数组在模n的意义下是等价的），因此可得两个数列分别的乘积对n同余

根据A数组的含义，A数组模n后的数组就是A数组本身

即：m₁m₂m₃m₄…m_Φ(n) ≡ a₁a₂a₃…a_Φ(n) (mod n)

将每项m展开的：a * a₁ * a * a₂ * a * a₃…a * a_Φ(n) ≡ a₁a₂a₃…a_Φ(n) (mod n)

化简得：a^Φ(n)(a₁a₂a₃…a_Φ(n)) ≡ a₁a₂a₃… a_Φ(n) (mod n)

因A数列中的每一项都跟n互质，那么整个数列的乘积也和n互质，所以

可以两边同除a₁a₂a₃…a_Φ(n)得：a^Φ(n) ≡ 1 (mod n)

证毕！

费马小定理
由欧拉定理可以马上推得费马定理。
当n为质数时，Φ(n) = n - 1
因此当p是质数时，a^p-1≡ 1 (mod p) 即费马定理！

3.快速幂O(logb)

算法核心：将b用二进制表示，然后对每一位进行01判断，进而累乘即可。

i)题目：求a^b%p的值

#include <iostream>

using namespace std;

typedef long long ll;

int qmi(int a , int b , int p)
{
    int ans = 1;
    while(b)
    {
        if(b & 1) ans = (ll)ans * a % p; //因为数据很大，所以边循环边取模
        b >>= 1;                         //如果末位是1，就将ans累乘
        a = (ll)a * a % p;               //同时对a取平方
    }
    
    return ans;
}

int main()
{
    int n;
    
    cin >> n;
    
    while(n--)
    {
        int a, b , p;
        cin >> a >> b >> p;
        cout << qmi(a , b , p) << endl;
    }
    return 0;
}
 

ii)快速幂求逆元

题目：给定n组a_i , p_i，其中p_i是质数,求a_i模p_i的乘法逆元，若逆元不存在则输出impossible。

逆元定义：若整数b，p互质，并且对于任意的整数 a，如果满足b|a=0，则存在一个整数x，使得a/b≡a∗x(mod p)，则称x为b的模p乘法逆元，记为b⁻¹(mod p)。

接下来进行推导：
若：a/b≡a∗x(mod p)
两边同除a,得：1/b ≡ x (mod p)
移项得：x*b ≡ 1 (mod p) ，即x满足该式子即可！

由费马小定理得：若p是质数，且b不是p的倍数， 则b^p-1 ≡ 1 (mod p)
因此x = b ^p-2 就是解。
当b%p = 0时 ，x*b ≡ 1 (mod p)显然无解。

//代码中的a是上面文字中的b
#include <iostream>

using namespace std;

typedef long long ll;
//求a的b次方
ll qmi(ll a , ll b , ll p)
{
    ll ans = 1;
    while(b)
    {
        if(b & 1) ans = ans * a % p;
        b >>= 1;
        a = a * a % p;
    }
    
    return ans;
}

int main()
{
    int n;
    cin >> n;
    while(n--)
    {
        ll a , p;
        cin >> a >> p;
        ll t = qmi(a , p-2 , p);
        if(a % p) cout << t << endl;
        else puts("impossible");
    }
    return 0;
}

4.扩展欧几里得算法

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

裴蜀定理：若a,b是整数,且gcd(a,b)=d，那么对于任意的整数x,y，ax+by都一定是d的倍数，特别地，一定存在整数x,y，使ax+by=d成立。

显然因为d|a , d|b , 则d|ax+by。

有一个直接的应用就是 如果ax+by=1有解，那么gcd(a,b)=1；

根据裴蜀定理，我们知道a_i∗x_i+b_i∗y_i=gcd(a_i,b_i)一定存在解，要是想求出x，y就得用到扩展欧几里得算法。
需要找出递归关系
求解方程 ax+by=gcd(a,b)
当 b=0 时， ax+by=a 故而 x=1,y=0。
当b != 0 时，因为gcd(a,b)=gcd(b,a%b)
而gcd(a,b) = xa+yb , gcd(b,a%b) = yb+x(a%b) . 其中 a%b = a - (a/b)b
所以 gcd(b , a%b) = yb+x(a - (a/b)b) , 化简得：gcd(b , a%b) = xa+(y-(a/b)x)b
推得：xa+yb = xa+(y-(a/b)x)b
根据恒等定理得：y = (y-(a/b) * x) 这就是递归关系。

#include <iostream>

using namespace std;

void exgcd(int a , int b , int &x , int &y)
{
    if(!b)//终止条件
    {
        x = 1 , y = 0;
        return;
    }
    exgcd(b , a % b , y , x);//这里将y写前面，x写后面可以简化代码，否则要进行交换
    y -= a / b * x;
    /*这是不交换x，y的写法
	exgcd(b, a%b, x, y);
    int t = y;
    y = x - (a/b) * y;
    x = t;
   */
    return ;
}


int main()
{
    int n;
    cin >> n;
    
    while(n--)
    {
        int a , b , x , y;
        cin >> a >> b;
        exgcd(a , b , x , y);
        printf("%d %d\n" , x , y);
    }
    return 0;
}

通解 = 特解 + 齐次解
而齐次解即为方程 ax+by=0的解
通解为 x=x′+k∗b/d , y=y′−k∗a/d ，k∈z

利用扩展欧几里德算法求解线性同余方程

题目：给定n组数据a_i,b_i,m_i，对于每组数求出一个x_i，使其满足a_i∗x_i≡b_i(mod m_i)，如果无解则输出impossible。
算法核心：将原方程进行等价代换：
1.ax≡b (mod m) 等价于 ax+my≡b (mod m)
2.即 (ax+my)%m=b
3.只要利用扩展欧几里德算法求出gcd(a,m)，如果b是gcd(a,m)的倍数，那么显然可以通过乘积的方式，来使ax+my等于b，否则返回impossible。

#include <iostream>

using namespace std;

typedef long long ll;

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;
    cin >> n;
    
    while(n--)
    {
        int a , b , m , x , y;
        cin >> a >> b >> m;
        int d = exgcd(a , m , x , y);//求gcd(a,m)
        if(b % d) puts("impossible");
        else cout << (ll) x * b / d % m << endl;//如果b是最大公约数d的倍数，则可以通过乘上一个b/d来使求得的结果由d变成b
    }
    return 0;
}

中国剩余定理

描述：给定 2n 个整数a₁,a₂,…,a_n和m₁,m₂,…,m_n,求一个最小的非负整数 x，满足∀_i∈[1,n],x≡m_i(mod a_i)。

思路：
合并：将第一条和第二条合并，再合并进第三条……一直合并n-1次，就只剩一条了。
具体实现：1.由题意得x % a₁ = m₁ , x % a₂ = m₂,
2.即x = k₁a₁ + m₁ , x = k₂a₂+m₂,
3.取出两式等号右边得：k₁a₁ + m₁ = k₂a₂+m₂
4.移项得:k₁a₁-k₂a₂ = m₂-m₁.( * )
5.接下来利用扩展欧几里得算法求得k1的解和最大公约数d，判断d是否能整除m₂-m₁，若不行说明无解，否则继续
6.再将k1自乘上一个(m2-m1)/d才是( * )的一个解，可以知道k1的通解是k₁+ka₂/d，k2的通解是k₂+ka₁/d
7.为了防止数据溢出，我们需要求得k1的最小正整数解，在c++中，取模会得到负值，因此我们需要以下操作
8.令t=a₂/d，k₁ = (k₁ % t + t) % t
9.将k₁= k₁+ka₂/d，代回x = k₁a₁ + m₁中得：x = k₁a₁+m₁+ka₁a₂/d
10.发现此时x的表达式和原来的很像，一部分常量，一部分变量，所以我们重新赋值 , m₁ = k₁a₁+m₁，a₁= abs(a₁ / d * a₂)
11.这样一来原来的两条式子就合并为一条x=k₁a₁+m₁（此时a₁和m₁已经变了），继续合并即可。
12.最后出m₁就是解，在10中可以发现m₁=k₁a₁+m₁就是x的值，不过在输出也要对其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;
    }
    int d = exgcd(b , a % b , y , x);
    y -= a / b * x;
    return d;
}

int main()
{
    int n;
    cin >> n;
    
    ll a1 , m1;
    cin >> a1 >> m1;
    
    bool flag = true;
    for(int i = 0 ; i < n - 1 ; i++)
    {
        ll a2 , m2;
        cin >> a2 >> m2;
        ll k1 , k2;
        ll d = exgcd(a1 , a2 , k1 , k2);
        
        if((m2 - m1) % d)
        {
            flag = false;
            break;
        }
        
        k1 *= (m2 - m1) / d;//此时求出来的k1只是众多k1的解中的一个而已
        ll t = a2 / d;//防止溢出，要求得k1的最小正整数解，易得k1+k*(a2/d)是k1的通解，因此使其对t 取模
        k1 = (k1 % t + t) % t;//在c++中，取模的值可得负数，在这样的操作之后会变成同等意义的正数
        
        m1 = a1 * k1 + m1;
        a1 = abs(a1 / d * a2);
    }
    
    if(flag) cout << (m1 % a1 + a1) % a1 << endl;
    else puts("-1");
    return 0;
}