《算法设计编程实验：大学程序设计课程与竞赛训练教材》——3.2 求解不定方程和同余方程的实验范例...

3.2　求解不定方程和同余方程的实验范例

公约数和同余问题是初等数论的两个核心内容，是求解许多数论问题的基础。本节将围绕这两个问题展开实验：
1)计算最大公约数和最大公约数的线性组合，在此基础上介绍求解不定方程的方法，并为求解同余方程作铺垫。
2)介绍求解同余方程和同余方程组的基本方法。
3.2.1　计算最大公约数和不定方程
整数a和b的最大公约数可通过欧几里得公式计算：gcd(a,b)=ba=0
gcd(b mod a,a)否则 　　证明：关键是证明gcd(a，b)与gcd(b，a mod b)可互相整除。设d＝gcd(a，b)，将a mod b化成a与b的线性组合 a mod b＝a-ab*b。由于d能整除a和b，因此d一定能整除a与b的线性组合，即d能整除(a mod b)。又因为ｄ能整除b和(a mod b)，显然d能整除gcd(b，a mod b)。同理可证gcd(b，a mod b)能整除gcd(a，b)。由于gcd(a，b)与gcd(b，a mod b)可互相整除，蕴涵着gcd(a，b)＝±gcd(b，a mod b)，因此欧几里得公式成立。
读者可直接按照欧几里得公式的递归定义编写计算最大公约数的程序。下面，我们将给出欧几里得公式的推广形式：
用一个线性组合来表示最大公约数，即d=gcd(a,b)=ax+by。计算出其中的整系数x和y (x与y可能为0或负数)。
分析：若b=0,则gcd(a,b)=a，x=1，y=0；若b≠0,则首先计算出d′=gcd(b,a mod b)和满足d′=bx′+(a mod b)y′的x′、y′。在这种情况下, 有d=gcd(a,b)=d′=gcd(b,a mod b)。
为了得到满足d=ax+by的x和y，我们将d′=gcd(b,a mod b)转化成线性组合d′=bx′+(a mod b)y′=bx′+a-ab*by′
=ay′+bx′-aby′　　按照d和d′两个线性组合对应的法则，选择x=y′＼y=x′-aby′就可以满足等式。显然，欧几里得推广公式是递归定义的，读者可以根据上述分析编写计算最大公约数的线性组合的递归程序。
应用欧几里得推广公式，可以求解如下形式的“二元一次”方程：
现给出不定方程ax+by+c=0中的a、b、c和变量x、y的值域x∈［x1,x2］、y∈［y1,y2］，问该范围内有多少组整数解？
由于x和y在值域内的解可能有多个，所以称为“二元一次”不定方程。
方法1：简单枚举
一一枚举值域范围内的每一个整数值，找出其中使不定方程成立的x和y。显然，这种方法需要计算（x2-x1）*（y2-y1）次不定方程，计算量完全取决于值域的大小。
方法2：应用欧几里得算法
初等数论给出了一种更为有效的方法，这种方法将整个计算过程分为两个步骤。
步骤1：运用欧几里得算法的推广形式计算ax+by+c=0的初始解x0、y0和gcd(a，b)。gcd(a,b)=b,x=0,y=cb(a=0)&&(c%b=0)
gcd(b%a,a)其他 将欧几里得算法对应于不定式，可得:(b mod a)x′+ay′=c将(b mod a) 转化为线性组合后得出：(b-baa)x′+ay′=c整理上式后得出：x′b+a(y′-bax′)=c上式与原不等式对应，得出：y=x′,x=y′-bax′　　如此一直递归下去,直至a=0、c mod b=0为止。此时gcd(a,b)=b,x=0,y=cbb 。
步骤2：根据不定方程的通解形式计算所有解。
如果排除值域范围的因素，则不定方程ax+by+c=0的通解形式为x=x0+t*Δx
y=y0+t*Δy 其中,
x0和y0为ax+by=c的初始解；Δx=lcm(a,b)a 或bgcd(a,b)；Δy=lcm(a,b)b 或agcd(a,b)称之为通解形式。将上述通解代入后，方程ax+by+c=0成立。显然，在x和y公共范围内系数t的个数即为值域范围内解的个数。将这些t代入通解后，即可得出该范围内的所有解。由此得出计算不定方程ax+by+c=0的一般方法：
1)运用欧几里得算法的推广形式先求出满足方程的一组特殊解(x0, y0)，这组解是否满足范围限制条件并不重要。于是设Δx=lcm(a,b)a, Δy=lcm(a,b)b。方程的通解就可以表示为x=x0+tΔx
y=y0－tΔy 　　　 t∈Z　　2)根据x和y取值范围的限制算出t的取值范围。
x的区间为［lx,rx］,其中lx=min{(x1-x0),(x2-x0)}Δx；rx=max{(x1-x0),(x2-x0)}Δx。
y的区间为［ly,ry］,其中ly=min{-(y1-y0),-(y2-y0)}Δy；ry=max{-(y1-y0),-(y2-y0)}Δy。
显然，系数t的取值范围为x区间和y区间的公共区间［max{lx,ly}，min{rx,ry} ］，该范围内的整数个数即为（x,y）的组数，将该取值范围内的t依次代入x=x0+tΔx
y=y0-tΔy 后，即可得出不定方程ax+by+c=0的所有解。
【3.2.1.1　The equation】
【问题描述】
给出方程ax+by+c=0，求出有多少对整数解(x, y)，满足X1≤x≤X2，Y1≤y≤Y2。
输入：
按顺序给出a、b、c、X1、X2、Y1、Y2，所有数的绝对值都小于108。
输出：

试题来源：《ACM》杂志15卷9期
在线测试地址：SGU 106
试题解析
方程若有整数解，必须满足如下条件：
若a=b=0，则c必须为0；否则c必须能够整除gcd(a，b)且x和y存在满足条件的整数区间，且两个区间有交集。
由此得出算法：
首先使用欧几里得算法的推广形式计算初始解x0、y0和gcd(a，b)。
1)直接求解a=b=0的特例（即gcd(a，b)==0）：
若常数c=0，则根据乘法原理，x和y取值范围内整数个数的乘积（（X2-X1+1)*(Y2-Y1+1)）为整数解的对数。
若常数c≠0，则方程不成立，无解退出。
2)若c不能整除gcd(a，b)，则方程无整数解，失败退出。
3)将方程ax+by+c=0转换为ax+by=-c。x0=x0-cgcd(a,b)，y0=y0-cgcd(a,b)
Δx=lcm(a,b)a或bgcd(a,b)，Δy=-lcm(a,b)b或-agcd(a,b)　　计算x所在的区间［lx，rx］=X1－x0Δx,X2－x0ΔxΔx>0
X2－x0Δx,X1－x0ΔxΔx≤0 。若lx>rx，则x无解，失败退出。
计算y所在的区间［ly，ry］=Y1－y0Δy,Y2－y0ΔyΔy>0
Y2－y0Δy,Y1－y0ΔyΔy≤0 。若ly>ry，则y无解，失败退出。
计算两个区间的交集［lp，rp］=［max{lx，ly}，min{rx,ry}］。若两个区间无交集(lp>rp)，则失败退出；否则整数解的对数为公共区间中整数的个数rp-lp+1。
程序清单

# include <iostream>
# include <cmath>
using namespace std;
typedef long long int64;

int64 exgcd(int64 a,int64 b,int64 &x,int64 &y)// 使用欧几里得算法的推广形式计算初始解x、y
// 和gcd(a，b)
　　if(!b){
　　　　　x=1,y=0;
　　　　　return a;
　　}else{
　　　　　int64 ret=exgcd(b,a%b,y,x);
　　　　　y-=a/b*x;
　　　　　return ret;
　　}
}

int64 a,b,c,X1,X2,Y1,Y2;
int main(){
　　cin>>a>>b>>c>>X1>>X2>>Y1>>Y2;// 方程系数为a、b、c；x的数值区间为［X1,X2］；y的
// 数值区间为［Y1,Y2］
　　c*=-1;// 常数取反
　　int64 x,y,r,dx,dy;
　　r=exgcd(a,b,x,y);// 计算gcd(a，b)和方程初始解
　　if(r==0){// 在a和b全为0的情况下，若常数c为0，则整数解的
// 对数为(X2-X1+1)*(Y2-Y1+1)；否则方程不成立
　　　　　if(c==0)　cout<<(X2-X1+1)*(Y2-Y1+1)<<endl;
　　　　　else cout<<0<<endl;
　　　　　return 0;
　　}
　　if(c%r!=0){// 若c不能整除gcd(a，b)，则方程无整数解
　　　　　cout<<0<<endl;
　　　　　return 0;
　　}
　　x*=c/r;y*=c/r;
　　dx=b/r,dy=a/r*-1;// 计算x,y的改变量

　　int64 lx,ly,rx,ry;
　　X1-=x,X2-=x,Y1-=y,Y2-=y;
　　if(dx>0)// 计算x所在的区间［lx，rx］
　　　　　　　lx=ceil(double(X1)/dx),rx=floor(double(X2)/dx);
　　　　else lx=ceil(double(X2)/dx),rx=floor(double(X1)/dx);
　　if(dy>0)// 计算y所在的区间［ly，ry］
　　　　　　　ly=ceil(double(Y1)/dy),ry=floor(double(Y2)/dy);
　　　　　else　ly=ceil(double(Y2)/dy),ry=floor(double(Y1)/dy);

　　int64 lp=max(lx,ly),rp=min(rx,ry);// 计算两个区间的交集［lp，rp］
　　if(rx<lx||ry<ly||rp<lp)　cout<<0<<endl;// 若x或y不存在满足条件的整数区间，或交集为空，
// 则方程不成立；否则返回交集的整数个数
　　　　else　cout<<rp-lp+1<<endl;
　　return 0;
}

3.2.2　计算同余方程和同余方程组
对于一组整数Z，Z里的每一个数都除以同一个数m，得到的余数可以为0，1，2，…m-1，共m种。我们就以余数的大小作为标准，将Z分为m类，每一类都有相同的余数。在每一类下的任意两个数a和b都关于m同余，记为a≡b(mod m)。
下面讨论同余方程和同余方程组的计算方法。
1.计算同余方程
诸如下述形式的方程ax≡b(mod m)（其中a、b和m都是整数）称为同余方程，方程中ax和b关于m同余。
显然，变量x的解可能有零个、一个或多个。可以简单地尝试，依次用x=O，1，…，m-1来代入该方程，找出其中在模m时满足该方程的整数x。但这种盲目搜索的效率完全取决于m的大小。初等数论给出了一种更为有效的方法:
1）通过欧几里得推广算法求出d=gcd(a,m)和两个满足d=ax′+my′的值x′和y′，表明x′是方程ax′≡d(mod m)的一个解。
2）d不能被b整除，则方程ax≡b(mod m)没有解；否则有d个解，其中第1个解的值为x0=x′bdmod m，其余d-1个解可以通过对模m加上md的倍数来得到，即xi=x0+imdmod m(1≤i≤d-1)。
证明：
证明中引用了以下三个定理。
【定理3.2.1】若ac≡bc(mod m)且gcd(c,m)=d, 则a≡bmodmd。
【定理3.2.2】设d≠0且ad≡bd(mod md), 则a≡b(mod m)。
【定理3.2.3】若gcd(a,m)=1，则一次同余式ax+b≡0(mod m)有解。
上述三个定理比较简单，我们省略了证明过程。读者可自行证明或查阅相关资料。下面，分两个步骤推导证明同余方程ax≡b(mod m)的计算方法。
步骤1：证明adx≡bd(modmd)与ax≡b(mod m)同解。
由定理3.2.1和定理3.2.2知, 同余式ax≡b(mod m)和同余式adx≡bdmodmd同解。由于gcd(a,m)=d＞1,因此ad与md互质。由定理3.2.3知,adx≡bdmodmd有解，即admodmd有同余类［x］。若取0≤x＜md,则［x］=x+kmdk=0,±1,±2,…且［x］还是ax≡b(mod m)的解。
鉴于［x］、x+md、x+2md、…、x+(d-1)md都在［x］之中,而0≤x+imdmd、…、x+(d-1)*md是ax≡b(mod m)的d个不同的解。
步骤2：证明ax≡b(mod m)只有［x］、 x+md、x+2md、 …、x+(d-1) md这d个不同的解。
从以上讨论知,可设x+tmd为ax≡b(mod m)的一个解,由于t≡i(mod d),i∈{0,1,2,…,d-1}，据定理3.2.2,在t≡i(mod d)两边同乘md便有tmd≡imd(mod m),这就是说,x+tmd为［x］、x+md、x+2md、…、x+(d-1) md中之一。
2.计算同余方程组
同余方程组a≡a1(mod n1)
a≡a2(mod n2)
...
...
...
a≡ak(mod nk）中每个模ni>0且每对模ni和nj互为质数(1≤i,j≤k,i≠j) ,a1～ak和n1～nk已知,要求解出a值，a值的唯一性可由中国剩余定理证明。问题是如何求a。
显然，若通过盲目搜索k个同余式来计算a的话，则时间效率太低。下面给出另一种算法：
应用数论知识将同余方程组转化为一个多项式a=(a1c1+…aici+…+akck)mod（n1n2…nk）由该多项式直接计算出a值。
那么如何证明这个多项式与同余方程组对应，如何计算式中变量ci（1≤i≤k）的值呢？
证明：因为正整数n1，n2，…，nk两两互素, 所以当i≠j时有gcd(ni,nj)=1, 且由于mi=nni, 故gcd(ni,mi)=1成立,i=1，2，…，k。于是由整数a、b的最大公约数d=gcd(a,b)可以表示为a、b的倍数和的性质（即存在整数x和y，使得d=ax+by）知,存在整数n′i和m′i,使得mim′i+nin′i=1,即存在整数m′i满足mim′i≡1(mod ni)　i=1,2,…,k①另一方面, 当i≠j时，由gcd(ni,nj)=1和mi=nni，得ni|mj, 因此有ajmjm′j≡0(mod ni)　 i,j=1,2,…,k②而由式①与②得a1m1m′1+a2m2m′2+…+akm′kmk≡aimim′i(mod ni)，aimim′i≡ai(mod ni), i=1、2、…、k由此证明了x=a1m1m′1+a2m2m′2+…+akmkm′k(mod n)是同余方程组的解，并引出了计算同余方程组的4个步骤。
步骤1：计算mi。设:n=n1 n2 ...*nk
m1=nn1=n2 n3…nk；m2=nn2=n1 n3n4…*nk ； …
mi=nni=n1 …ni-1 ni+1…nk ； …
mk=nnk=n1 …nk-2 nk-1 　　步骤2：计算mi中关于ni的乘法逆元m-1i（即满足mim-1i≡1(mod ni)的m-1i）。有两种方法计算m-1i：
1)利用同余方程。
由于mi与ni互为质数(即最大公约数为1)，因此可以通过m1m-11≡1(mod n1)；…；mim-1i≡1(mod ni)；…；mkm-1k≡1(mod nk)得出m-11,…,m-1i,…,m-1k。方程mim-1i≡1(mod ni)对模ni有唯一解m-1i(1≤i≤k)。
2)利用欧几里得算法。
求出欧几里得推广公式gcd(ni,mi)=nix+miy＝1中的x和y,此时的y即为m-1i(1≤i≤k)。
步骤3：进行k次模运算ci=mi*(m-1imod ni) (1≤i≤k),得出c1,…，ck。
步骤4：最后将c1，…，ck代入方程a=(a1c1+…aici+…+ak*ck)mod n,便可求出a值。
【3.2.2.1　C Looooops】
【问题描述】
给出一个C语言风格类型的循环：

for (variable = A; variable != B; variable += C)

　　statement;即一个循环，在开始的时候将A 赋给变量，当变量不等于B 时，重复语句，然后对变量增加C。对于特定值A、B和C，我们要知道语句执行多少次，本题设定所有的算术运算都是在以2k为模的k位无符号整数类型（即值x满足0 ≤ x<2k）上进行。
输入：
输入包含若干测试用例，每个测试用例一行，给出用一个空格分开的4个整数A、B、C、k，整数k(1≤ k ≤32)是循环控制变量的二进制位数，而A、B、C(0≤ A, B, C < 2k)是循环参数。
输入以包含4个0的一行结束。
输出：
输出相应于输入实例由若干行组成。第i行给出第i 个测试用例中语句的执行次数（一个整数），如果循环不终止，则输出单词“FOREVER”。

试题来源：CTU Open 2004
在线测试地址：POJ 2115，ZOJ 2305
试题解析
按照题意，循环变量的初值为A，终值为B，每次循环的增量为C。由于所有算术运算都是在以2k为模的k位无符号整数类型（0≤x <2k）上进行，因此令D=(B-A)mod 2k 等价于 x*C≡D(mod 2k)显然，循环次数为0当且仅当D=(B-A)mod 2k=0。
要使得同余方程xC≡D(mod 2k) 有解，当且仅当gcd(c，2k)│D。通过欧几里得推广公式求方程xC+y2k=gcd(c，2k)中x的最小非负整数解，即x是同余方程Cx≡gcd(c，2k)(mod 2k)的一个解。若gcd(c，2k) 不能被D整除，则方程xC≡D(mod 2k) 无解，程序陷入死循环；否则(xD)mod 2k为方程xC≡D(mod 2k)的解，亦为语句的执行次数。
程序清单

#include<cmath>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#define ll long long
#include<iostream>
using namespace std;
ll exgcd(ll a,ll b,ll &x,ll &y){// 欧几里得推广：计算d=gcd(a,b)=ax+by的整系数
// x和y (X与Y可能为0或负数)
　　if(b==0){
　　　　　x=1;y=0;return a;
　　}
　　ll t=exgcd(b,a%b,y,x);
　　y-=a/b*x;
　　return t;
}
ll gcd(ll a,ll b){// 欧几里得公式：返回a和b的最大公约数
　　if(b==0)return a;
　　return gcd(b,a%b);
}
int main () {
　　int A,B,C,K;
　　ll i,j,ans;
　　while(1){
　　　　scanf("%d%d%d%d",&A,&B,&C,&K);// 反复输入循环参数和循环变量的二进制位数
　　　　if(!A&&!B&&!C&&!K)break;// 若输入4个0，则退出
　　　　ll a,b,c,k;
　　　　a=A,b=B,c=C,k=K;
　　　　ll d=b-a;// 计算(b-a)%2k的非负整数d。若d=0，则循环
// 次数为0；若d%gcd(c，2k) ≠0，则陷入死循环
　　　　k=(1ll)<<k;
　　　　d%=k;
　　　　if(d<0)d+=k;
　　　　if(d==0) puts("0");continue;
　　　　ll tem=gcd(c,k);
　　　　if(d%tem) puts("FOREVER");continue;
　　　　c/=tem,k/=tem,d/=tem;
　　　　exgcd(c,k,ans,j);// 计算gcd(c，k)=c*ans+k*j的一个解ans 
　　　　ans*=(d);// (ans*d)%k的非负整数即为语句执行的次数
　　　　ans%=k;
　　　　if(ans<0)ans+=k;
　　　　cout<<ans<<endl;
　　}
　　return 0;
}

【3.2.2.2　Biorhythms】
【问题描述】
人生来就有三个生理周期，分别为体力、感情和智力周期，它们的周期长度为23天、28天和33天。每一个周期中有一天是高峰。在高峰这天，人会在相应方面表现出色。例如，在智力周期的高峰，人会思维敏捷，精力容易高度集中。因为三个周期的周长不同，所以通常三个周期的高峰不会落在同一天。对于每个人，我们想知道何时三个高峰落在同一天。对于每个周期，我们会给出从当前年份的第一天开始到出现高峰的天数（不一定是第一次高峰出现的时间）。你的任务是给定一个从当年第一天开始数的天数，输出从给定时间开始（不包括给定时间）下一次三个高峰落在同一天的时间（距给定时间的天数）。例如：给定时间为10，下次出现三个高峰在同一天的时间是12，则输出2（注意这里不是3）。
输入：
输入四个整数：p、e、i和d。 p、e、i分别表示体力、情感和智力高峰出现的时间（时间从当年的第一天开始计算）。d是给定的时间，可能小于p、e或 i。 所有给定时间是非负的并且小于365，所求的时间小于21252。
当p=e=i=d=-1时，输入数据结束。
输出：
从给定时间起，下一次出现三个高峰在同一天的时间（距离给定时间的天数）。
采用以下格式：

Case 1: the next triple peak occurs in 1234 days.

注意：即使结果是1天，也使用复数形式“days”。

试题来源：ACM East Central North America 1999
在线测试地址：POJ 1006，ZOJ 1160，UVA 756
试题解析
体力、感情和智力3个周期的长度为23天、28天和33天，这3个周期长度两两互质。假设第x天三个高峰同时出现，则可得到同余方程组x≡p(mod 23)
x≡e(mod 28)
x≡i(mod 33)根据中国剩余定理，x在232833=21252的范围内有唯一解。令ai和ni分别为三个高峰出现的时间和周期长度,即a1=p、a2=e、a3=i、n1=23、n2=28、n3=33。得到同余方程组x≡ai(mod ni)，(1≤i≤3)　　运用上述4个步骤求出s=∑3i=1miai(m－1imod ni)，其中m1=2833,m2=2333,m3=23*28，m－1i为mi中关于ni的乘法逆元，即mim－1i≡1(%ni)，乘法逆元可以通过欧几里得扩展公式求得。
由于三个高峰同时出现的时间是相隔给定时间d的天数，因此这个时间应为(s-d)mod n的最小正整数（n=232833）。
程序清单

#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<string>
using namespace std;
typedef long long ll;
ll power(ll a,ll p,ll mo){// 通过反复平方法计算ap%(mo)
　　ll ans=1;
　　　for(;p;p>>=1){
　　　　　if(p&1){
　　　　　　　　ans*=a;
　　　　　　　　if(mo>0)ans%=mo;
　　　　　}
　　　　　a*=a;
　　　　　if(mo>0)a%=mo;
　　　}
　　　return ans;
}
ll exgcd(ll a,ll b,ll &x,ll &y){ 　　　// 欧几里得推广公式：计算方程gcd(a,b)=ax+by中变量x的值
　　　if(b==0){
　　　　　x=1;y=0;return a;
　　　}
　　　ll t=exgcd(b,a%b,y,x);
　　　y-=a/b*x;
　　　return t;
}

ll niyuan(ll a,ll p){// 计算a-1%p
　　　ll x,y;
　　　exgcd(a,p,x,y);// 计算同余方程ax≡gcd(a,p)(% p)中的x
　　　return (x%p+p)%p;// x对p的模取正
}
int main(){
　　　int　a,b,c,d,i,j,k,u,v,te=0;
　　　while(1){
　　　　scanf("%d%d%d%d",&a,&b,&c,&d);// 反复输入体力、情感和智力高峰出现的时间和给定时间
　　　　if(a==b&&b==c&&c==d&&a==-1)break;// 结束标志
// 计算an=（∑3i=1(mi*ai*(m－1imod ni))-d)%（23*28*33）的非负整数，该数即为下一次三个高峰同天的时间
　　　　ll an=0;
　　　　an=28*33*a*niyuan(28*33,23)+23*33*b*niyuan(23*33,28)+23*28*c*niyuan(28*23,33);
　　　　an-=d;
　　　　an%=(28*33*23);
　　　　if(an<=0)an+=28*33*23;
　　　　printf("Case %d: the next triple peak occurs in %d days.＼n",++te,(int)an);
　　　}
}