[量水问题]: 有三个分别装有a升水,b升水,c升水的量筒,其中a,b互质,c>b>a>0,现在c筒装满水,问能否在c筒中量出d升水(c>d>0)。若可以,给出方案。
解答:
所谓模数方程,就是模线性方程,即形如 ax ≡ b (mod c) 形式的方程,其中a,b,c是常数,x是自变量,这个方程表示ax mod c = b mod c,即ax和b模c同余。 这个量水问题,用模数方程解比较方便,具体算法分析如下。
量水过程实际上就是倒来倒去,每次倒的时候总有如下几个特点: 1。总有一个筒中的水没有变动; 2。不是一个筒被倒满酒是另一个筒被倒光; 3。c筒仅起到中转作用,而本身的容积除了必须足够装下a筒和b筒全部的水以外,别无其他的限制; 这样,假设整个倒水过程中对a筒倒满了x次,对b筒倒满了y次,则: ax + by = d, (1) 上式的x,y为整数,而且既可以是正整数(表示该筒(a或b)被c筒的水倒满的次数),也可以是负整数(表示该筒(a或b)被倒满后又倒回到c筒的次数)。 一般可以用穷举法来解方程(1),但是这种方法局限性很大。我们可以将方程转化为: ax = -by + d, 进一步变为 ax ≡ d (mod b), (20 这样问题就变成了求解模数方程(2)的解的问题。解x的个数就是可行方案的总数。其中xi表示第i种方案中a筒倒满的次数,xi代入方程(1)后求出来的yi表示b筒倒满的次数。
例如:有三个量筒,a=3,b=7,c=10,求c筒中量出5升水的所有方案。 解方程 : 3x ≡ 5 (mod 7) 得到 x1 = 4, y1=-1 和x2=-3, y2=2, 这两个解对应的倒水方案如下:
方案一: x1=4, y1=-1
次数 a b c 说明 1 0 0 10 初始状态 2 3 0 7 从c倒水到a中,把a倒满 3 0 3 7 从a倒水到b中,把a倒空 4 3 3 4 从c倒水到a中,把a倒满 5 0 6 4 从a倒水到b中,把a倒空 6 3 6 1 从c倒水到a中,把a倒满 7 2 7 1 从a倒水到b中,把b倒满 8 2 0 8 从b倒水到c中,把b倒空 9 0 2 8 从a倒水到b中,把a倒空 10 3 2 5 从c倒水到a中,把a倒满 11 0 5 5 从a倒水到b中,把a倒空
整个过程中,一共有4次“从c倒水到a中,把a倒满”,有1次“从b倒水到c中,把b倒空”;
方案二: x2=-3, y2=2
次数 a b c 说明 1 0 0 10 初始状态 2 0 7 3 从c倒水到b中,把b倒满 3 3 4 3 从b倒水到a中,把a倒满 4 0 4 6 从a倒水到c中,把a倒空 5 3 1 6 从b倒水到a中,把a倒满 6 0 1 9 从a倒水到c中,把a倒空 7 1 0 9 从b倒水到a中,把b倒空 8 1 7 2 从c倒水到b中,把b倒满 9 3 5 2 从b倒水到a中,把a倒满 10 0 5 5 从a倒水到c中,把a倒空
整个过程中,一共有3次“从a倒水到c中,把a倒空”,有2次“从c倒水到b中,把b倒满”;
至于其中解模数方程的方法,一些关于数论的书上有介绍,说起来也比较麻烦,有很多的定理公式,我直接给出一个程序吧:
========================================================== c 语言的程序: ==========================================================
/********************************************
求解模线性方程 ax=b (mod n) ,n>0 copyright starfish 2000/10/24
*********************************************/
void modular_linear_equation_solver(int a, int b, int n) { int e,i,d; int x,y; d = ext_euclid(a, n, x, y); if (b%d>0) { printf("No answer!\n"); } else { e=(x*(b/d))%n; for (i=0; i<d; i++) //notice! Here x maybe less than zero! printf("The %dth answer is : %ld\n",i+1,(e+i*(n/d))%n); } }
其中用到的ext_euclid 函数如下:
/*********************************************
扩展欧几里德算法求gcd(a,b)=ax+by
copyright starfish 2000/10/24
*********************************************/
int ext_euclid(int a,int b,int &x,int &y) { int t,d; if (b==0) {x=1;y=0;return a;} d=ext_euclid(b,a %b,x,y); t=x; x=y; y=t-a/b*y; return d; }
========================================================== pascal语言的程序: ==========================================================
{=== 扩展的欧几里德算法,求出gcd(a,b)和满足gcd(a,b)=ax+by的整数x和y ===} function extended_euclid(a,b:longint;var x,y:longint):longint; var t:longint; begin if b=0 then begin result:=a; x:=1; y:=0; end else begin result:=extended_euclid(b,a mod b,x,y); t:=x; x:=y; y:=t-(a div b)*y; end; end;
{=== 求解模线性方程 ax ≡ b (mod n) 其中n>0 ===} procedure modular_linear_equation_solver(a,b,n:longint); var d,x,y,e,i:longint; begin d:=extended_euclid(a,n,x,y); if b mod d>0 then writeln('No answer!') {输出无解信息} else begin e:=x*(b div d)mod n; for i:=0 to d-1 do writeln( (e+i*(n div d)) mod n ); {输出第i个解 } end; end; 
|