扩展欧几里得和中国剩余定理

贴一个dalao的证明：欧几里德算法与扩展欧几里德算法

1.欧几里得算法 O(logn)

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

2.贝祖定理

正整数a,b，那么一定存在非零整数x,y使得ax+by=gcd(a,b) 构造法证明——扩展欧几里得算法

贝祖定理扩展：如果 ax+by=n有解，充要条件是 n 是gcd(a,b)的倍数

a是最大公约数的倍数，b是最大公约数的倍数，他们乘以一个系数相加，也应该是gcd(a,b)的倍数

3.扩展欧几里得算法

$ax+by=d$，$d=(a,b)=(b,a\%b)=(b,a-\lfloor \frac{a}{b}\rfloor \ast b)=(b,a-\lfloor \frac{a}{b}\rfloor \ast b+\lfloor \frac{a}{b}\rfloor \ast b)=(b,a)$

$by+(a\%b)x=d$
$by+(a-\lfloor \frac{a}{b}\rfloor \ast b)x=d$
$ax+b(y-\lfloor \frac{a}{b}\rfloor \ast x)=d$，所以x是不需要变的，只需要y=y-a/b*x;

int exgcd(int a,int b,int &x,int &y)
{
	if(b==0)//如果b为0
	{
		x=1,y=0;
		return a;//(a,0)=a a*x+0*y=a x=1,y=0就是其中一组解
	}
	int d=exgcd(b,a%b,y,x);
	y=y-a/b*x;
	return d;
}

所求的x,y是不唯一的，只要求出了一个 $x_0,y_0$ 就可以求出所有的解

$a(x-\frac{b}{d})+b(y+\frac{a}{d})=d$

通解：$x=x_0-\frac{b}{d}k，y=y_0+\frac{a}{d}k(k\in Z)$

求解线性同余方程

输入

2
2 3 6
4 3 5

输出

impossible
7

思路

eg.    2x≡3(mod 6)无解
    4x≡3(mod 5)   x=2;

$ax\equiv b(modb)<=>存在整数y，使得ax=my+b<=>ax-my=b$

令y’=y，则有$ax+m{y}^{\prime }=b$，该方程有解的等价条件是 (a,m)|b

$ax+m{y}^{\prime }=gcd(a,m)，b=k\ast d，d=gcd(a,m)$

$ax+m{y}^{\prime }=d$可用扩展欧几里得算法来求，再让等式两边同时扩大$\frac{b}{d}$倍

#include <iostream>
#include <algorithm>
#include <unordered_map>
#include <vector>
const int maxn=1e6+5;
const int mod=1e9+7;
typedef long long ll;
using namespace std;
int exgcd(int a,int b,int &x,int &y)
{
	if(b==0)//如果b为0
	{
		x=1,y=0;
		return a;//(a,0)=a a*x+0*y=a x=1,y=0就是其中一组解
	}
	int d=exgcd(b,a%b,y,x);
	y=y-a/b*x;
	return d;
}
int main()
{
	int n;
    scanf("%d",&n);
	while(n--)//数据量有点大，建议用scanf
	{
		int a,b,m;
		scanf("%d%d%d",&a,&b,&m);
		int x,y;
		int d=exgcd(a,m,x,y);//注意这里应该对应a,m
		ll ans=1ll*b/d*x%m;//注意不能写成x/d*b，因为d不一定能整除x
		if(b%d)
		printf("impossible\n");
		else
		printf("%lld\n",ans);
	}
    return 0;
}

注意ll ans=1ll*b/d*x%m;要1ll*，不然会溢出。

4.中国剩余定理

所以$m_i$和$n_i$互质，$m_i^{-1}$表示$m_i$模$n_i$的逆

$x=a_1{m}_1{m}_1^{-1}+a_2{m}_2{m}_2^{-1}+...+a_k{m}_k{m}_k^{-1}$，用扩展欧几里得来求逆($ax\equiv 1modn$)

上式成立原因：
对于$n_1$来说，$m_1和n_1$互质，$(m_1{m}_1^{-1})mod{n}_1=1$，后面其他项$(m_i{m}_i^{-1})mod{n}_1$均为0，因为$m_i$中包含$n_i$，所以就可以得到第一项$x\equiv a_{1}mod{n}_1$

其他同理

例题

输入

2
8 7
11 9

输出

31

思路

$xmod{a}_1=m_1$  $x=k_1{a}_1+m_1$
$xmod{a}_1=m_1$  $x=k_2{a}_2+m_2$
要满足所有方程首先要满足前两个方程，然后我们可以两两合并地去求解

$k_1{a}_1+m_1=k_2{a}_2+m_2=>k_1{a}_1-k_2{a}_2=m_2-m_1$，已知$a_1,a_2,m_2-m_1$，要求$k_1,k_2$用扩展欧几里得，上式有解等价于 $(a_1,a_2)|m_2-m_1$

一旦我们求出了一组解$k_1,k_2$，通解就可表示为$k_1^{\prime }=k_1+k\frac{a_2}{d}，k_2^{\prime }=k_2+k\frac{a_1}{d}$ 前面加一个后面加一个，一减就抵消了

$x=k_1{a}_1+m_1=(k_1+k\frac{a_2}{d})a_1+m_1=a_1{k}_1+m_1+k\frac{a_1{a}_2}{d}$ x的所有解

$x_0=a_1{k}_1+m_1$，$a_1,a_2$的最小公倍数$a=\frac{a_1{a}_2}{d}$  => $x=x_0+ka$（这就是前两个式子合并后的结果）合并n-1次，就可以把n个方程组合成一个了

当只有一个方程时，要求最小非负整数x：
最后要求的就是$xmoda\equiv x_0$，$x和x_0$模a同余，就是求$x_{0}moda$的正余数

最后的a会是a1...an的最小公倍数，很有可能会溢出整数范围，所以一般方程式的个数不会太多

#include <iostream>
#include <algorithm>
#include <unordered_map>
#include <vector>
const int maxn=1e6+5;
const int mod=1e9+7;
typedef long long ll;
using namespace std;
ll exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0)//如果b为0
	{
		x=1,y=0;
		return a;//(a,0)=a a*x+0*y=a x=1,y=0就是其中一组解
	}
	ll d=exgcd(b,a%b,y,x);
	y=y-a/b*x;
	return d;
}
int main()
{
	int n;
    scanf("%d",&n);
	ll x=0,m1,a1;
	cin>>a1>>m1;
	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)
		{
			x=-1;
			break;
		}
		//exgcd求的是k1a1-k2a2=d的解，要把右部变成(m2-m1)的解，等式两边同乘(m2-m1)/d
		k1*=(m2-m1)/d;
		ll t=a2/d;//k1'=k1+a2/d*k
		k1=(k1%t+t)%t;//把k1变成这个方程的最小正整数解
		x=k1*a1+m1;//x0=k1a1+m1
		ll a=abs(a1/d*a2);//最小公倍数 a=a1a1/d
		//x=x0+ka   x≡m(%a) 构建新的迭代
		m1=k1*a1+m1;//m1=x0;
		a1=a;
	}
	if(x!=-1)
	x=(x%a1+a1)%a1;//取x%a
	
	cout<<x<<endl;
    return 0;
}

这道题数据比较极限，要让每个k都变得尽可能地小，所以要k1=(k1%t+t)%t;