浅谈基础数论算法
浅谈基础数论算法
数论，作为离散数学一个分支，是计算机科学的重要组成部分。在网络安全等方面有着广泛的应用。今天笔者将从基础内容出发，简要介绍最大公约数，解二元一次不定方程，解同余方程组等算法知识。
最大公约数
最大公约数，即两个整数的最大公因子。在计算机科学中，我们通常采取辗转相除的方法来计算两个数的最大公约数，代码如下。
int gcd(int a,int b)
{
if(!b)return a;
else return gcd(b,a%b);
}
这一算法的正确性源于一个简单的数论定理\(gcd(a,b)=gcd(b,a-b)\)，由此基本定理进行简单归纳可知\(gcd(a,b)=gcd(b,a\%b)\)。这两个定理的证明都较为基础，留给读者作为练习。
读者读到这里可能会对该算法的时间复杂度抱有一些疑惑。事实上，我们进行简单分析即可发现，在\((a,b)->(b,a\%b)\)的过程中，我们有\(a\%b \leq \lfloor \frac{a}{2} \rfloor\)，即其中一个数字会衰减到原来的一般。由此，我们可以得知，最多经过\(log_2{(a+b)}\)轮，两个数字其中之一就会变为\(0\)，此时算法也就运行结束，故该算法的时间复杂度为\(O(log_2n)\)
二元一次不定方程
二元一次不定方程指的是形如\(ax+by=c\)的方程，其中\(a,b,c\)为参数，\(x,y\)为变量，取值范围均为整数。
在计算机科学中，经常会遇到求解这类方程的特解，通解，正整数解，解的数量等问题。而扩展欧几里得算法\((exgcd)\)给出了一个求解该类方程通解的方法。
具体的算法流程是这样的，我们考虑沿用求两个数字最大公约数时的辗转相除法，迭代构造求解。
首先，在递归出口\(b=0\)时，\(d=gcd(a,b)=a\)，原方程转化为\(ax=a\)，此时我们取\(x_0=1,y_0=0\)即可得到一组特解。我们设在递归上一层时的参数为\(A,B\)，我们有\(a=B,b=A\%B\)。由刚才构造的特解可知，我们有方程\(B*x_0+(A\%B)*y_0=d\)成立。我们尝试构造\(x,y\)使得新的方程\(A*x+B*y=d\)。经过推导，我们可以得到以下表达式：
代入验证不难发现这组新的\((x,y)\)一定可以使得\(A*x+B*y=d\)成立。因此，我们只需要回溯即可递推出原方程的一组特解，设其为\(X_0,Y_0\)，我们可以证明，对于原方程\(ax+by=1\)，通解形式如下：
证明较为复杂，读者可以自行查看相关资料。
如下程序展示了求解二元一次不定方程非负整数解的方法，供以参考。
int X,Y;
int exgcd(int a,int b)
{
if(!b){X=1;Y=0;return a;}
int d=exgcd(b,a%b),t=X;X=Y;Y=t-(a/b)*Y;
return d;
}
int main()
{
int a=read(),b=read(),c=read();//ax+by=c
int d=exgcd(a,b);//d is the gcd of a and b
if(c%d!=0)printf("No solution");
else
{
X*=c/d;Y*=c/d;
int p1=b/d,p2=a/d;//x=k*X+t*p1 y=k*Y-t*p2
X=(X%p1+p1)%p1;Y=(c-a*X)/b;
if(Y<0)printf("There is no positive integer solution");
else printf("%d %d\n",X,Y);
}
return 0;
}
扩展欧几里得算法还可以同时计算正整数解数量，\(x\)最大的正整数解，\(y\)最大的正整数解等问题，限于篇幅限制，本文不再赘述，具体可以参考这道例题https://www.luogu.com.cn/problem/P5656。
利用这一算法，我们可以解决许多其他问题。比如求整数\(a\)在模\(p\)意义下的除法逆元\(x\)，即\(ax \equiv 1 \pmod{p}\)，我们只需要把该式转化为\(ax+kp=1\)即可。此外，该算法在接下来介绍的解同余方程组中也有应用。
同余方程组
同余方程组指的是对于许多条形如\(x\equiv a \pmod{p}\)的方程，去计算它们解集的交集。
我们考虑对于两条方程
等价于：
联立消元得：
这显然是一个关于\(k_1,k_2\)的二元一次不定方程，套用扩展欧几里得算法即可得到一组\(k_1,k_2\)的特解。
然后根据中国剩余定理，我们可以得知，\(x \equiv a_1+k_1*p_1 \pmod{lcm(p_1,p_2)}\)，这样我们就将两个方程合并成了一个方程。只需要不断沿用此方法，即可将方程组两两合并成一个方程，得到方程组的通解。
代码如下：
ll ksc(ll x,ll y,ll p){return ((x*y-(ll)(((ldb)x/p*y))*p)%p+p)%p;}//防止溢出的乘法
ll X,Y;
ll exgcd(ll a,ll b)
{
if(!b){X=1,Y=0;return a;}
ll d=exgcd(b,a%b),t=X;X=Y;Y=t-(a/b)*Y;
return d;
}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b){return (a/gcd(a,b))*b;}
int main()
{
// x=a (mod p)
ll n=read(),p,a;
for(ll i=1;i<=n;i++)
{
ll tp=read(),ta=read();
if(i==1){p=tp;a=ta;continue;}
ll np=lcm(p,tp),d=exgcd(p,tp);
X=ksc(X,(ta-a)/d,tp/d);
a=(ksc(X,p,np)+a)%np;p=np;//合并结果
}
printf("%lld",a);
return 0;
}
