数论相关：同余方程与同余方程组的解法

同余方程

形如 (ax equiv b pmod n) 的式子称为线性同余方程。对于这样的式子有解的充要条件是 (gcd(a,n) mid b) .

于是扩展gcd求解
将原方程化为一次不定方程 (ax+ny = b) .
利用扩展欧几里得算法求解不定方程 $ ax + ny = b$ 的整数解的求解全过程，步骤如下：

1、先计算 (gcd(a,n))，若 (b) 不能被 (gcd(a,n)) 整除，则方程无整数解；否则，在方程右边除以 (b/gcd(a,n))，记 得到新的不定方程 (ax_0 + ny_0 = gcd(a,n)).

2、利用扩展欧几里德算法求出方程 $ax_0 + ny_0 = gcd(a, b) $的一组整数解 (x_0) , (y_0)；

3、根据数论中的相关定理，记 (k=b/gcd(a,n)),可得方程 (ax + ny = b) 的所有整数解为：

[ x = k*x_0 + n/gcd(a,n)* t \ y =k* y_0 –a/gcd(a,n)* t \ (t=0,1,2，……)]

调整得到关于 (x) 的正整数解
注意因为解有多个，而我们要求最优解(正整数中最小的)，所以 ((x+=n/gcd(a,n)\%(n/gcd(a,n)));
加法是为了保证正数,取模是为了最小

青蛙的约会

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cstdlib>
#include <cmath>
using namespace std;
long long  init(){
	long long  rv=0,fh=1;
	char c=getchar();
	while(c<'0'||c>'9'){
		if(c=='-') fh=-1;
		c=getchar();
	}
	while(c>='0'&&c<='9'){
		rv=(rv<<1)+(rv<<3)+c-'0';
		c=getchar();
	}
	return fh*rv;
}
long long x,y,m,n,l;
long long exgcd(long long a,long long b,long long &x,long long &y){
	if(b==0){
		x=1;y=0;return a;
	}
	long long  t=exgcd(b,a%b,y,x);
	y-=(a/b)*x;
	return t;
}
int main(){
	freopen("in.txt","r",stdin);
	x=init();y=init();m=init();n=init();l=init();	
	if(m - n< 0) swap(m, n), swap(x, y);
	long long  a=0,b=0;
	long long  t=exgcd(m-n,l,a,b);
	if(n==m||(x-y)%t!=0){
		printf("Impossible");
		return 0;
	}
	(a*=(y-x)/t)%=(l/t);
	(a+=l)%=(l/t);//以保证最优解
	cout<<a;
	fclose(stdin);
	return 0;
}

exgcd可以用来求逆元

(ax equiv 1 pmod n) 已知 (a) , (n) 求 (x)
因为 (n) 是个素数，所以 (gcd(a,n)==1) ;
原方程可化为 (ax equiv gcd(a,n) pmod n)
用exgcd求解即可。

同余方程组

(x\%p_1 = b_1)
(x\%p_2 = b_2)
(x\%p_3 = b_3)
(x\%p_4 = b_4)
求 (x) 的最小正整数解
小范围数据直接枚举

对于模数互质的情况，使用中国剩余定理（CRT）
令 (m=p_1p_2p_3 ldots p_n)
构造出
//定义 (ni(k,p)) 是 (k) 在模 (p) 意义下逆元

[x = (m/p_1*ni(m/p_1, p_1)*a1 + ldots) \% m ]

CRT求解同余方程组

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
long long a1,a2,a3,a4,b1,b2,b3,b4;
long long m;
long long exgcd(long long a,long long b,long long &x,long long &y){
	if(!b){
		x=1;y=0;
		return a;
	}
	long long t=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return t;
}
long long ni(long long a,long long b){
	long long x=0,y=0;
	long long t=exgcd(a,b,x,y);
	if(t!=1) return -1;
	else return((x+b)%b);
}
int main(){
	freopen("in.txt","r",stdin);
	cin>>a1>>b1>>a2>>b2>>a3>>b3>>a4>>b4;
	m=a1*a2*a3*a4;
	cout<<(m/a1*ni(m/a1,a1)*b1+m/a2*ni(m/a2,a2)*b2+m/a3*ni(m/a3,a3)*b3+m/a4*ni(m/a4,a4)*b4)%m;
	fclose(stdin);
	return 0;
}

对于一般情况采用exgcd两两合并,
(x+a_1k_1=b_1)
(x+k_2a_2=b_2)
(a_1k_1-k_2a_2=b_1-b_2)
$t=exgcd(a_1,-a_2,k_1,k_2) $ 实际上 (-a_2) 可以写作 (a_2)
合并：
(k_1=(k_1*(b_1-b_2)/t)\%a2); //此处是模 (a_2) ,因为可以看成是模 (a_2)意义下的同余方程
(b_1-=a_1*k_1) // (b_1) 就是原式中的 (x)
(a_1=a_1/t*a_2) //把 $ a_1(变成)lcm(a1,a2)( )b_1%=a_1$ //把 (b_1) 调整至新式子的B
此时就把两个式子合并为了一个，待所有的都合并完后，结果就是 (b_1) 调整好的最小正整数 ((b_1+=a_1)\%=a_1)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <cstdlib>
#define LL long long 
using namespace std;
LL init(){
	LL rv=0,fh=1;
	char c=getchar();
	while(c<'0'||c>'9'){
		if(c=='-') fh=-1;
		c=getchar();
	}
	while(c>='0'&&c<='9'){
		rv=(rv<<1)+(rv<<3)+c-'0';
		c=getchar();
	}
	return rv*fh;
}
LL a,b,a1,b1;
LL exgcd(LL a,LL b,LL &x,LL &y){
	if(!b){
		x=1;y=0;return a;
	}
	LL t=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return t;
}
int main(){
	freopen("in.txt","r",stdin);
	a=init();b=init();
	for(int i=1;i<=3;i++){
		a1=init();b1=init();
		LL x=0,y=0;
		LL t=exgcd(a,a1,x,y);
		x=(x*(b-b1)/t)%a1;
		b-=a*x;
		a=a/t*a1;
		b%=a;
	}
	(b+=a)%=a;
	cout<<b;
	fclose(stdin);
	return 0;
}

求逆元

(ax equiv 1 pmod n)
逆元存在的充要条件是(gcd(a,n)==1);
一般采用exgcd求逆元，若p是质数，也可使用费马小定理，快速幂
模质数 (P) 意义下

从(1! sim n!)

先用快速幂处理出 (n!^{-1}),并预处理出来,(1! sim n!),那么

[n^{-1} = n!^{-1} * (n - 1)! ]

[(n-1)!^{-1} = n!^{-1} * n ]

由此递推即可

从 $1 sim n $

此方法不需要求阶乘,代码简单
首先 (1^{-1} = 1pmod p)
设 $p = k * i + r $, (1 < i < p, r < i)
所以 (k * i + r equiv 0 pmod p)
所有两边同乘以 (i^{-1}*r^{-1})得

[k * r^{-1} + i^{-1} equiv 0 pmod p\ i^{-1} equiv -k * r^{-1} pmod p\ i^{-1} equiv -lfloor p/r floor * ( pmod i)^{-1} ]

N[i] = (p -p / i) * N[p % i] %p;