倒水问题C++实现

量水问题方案

[量水问题]：
有三个分别装有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;