直线上的点(扩展欧几里得算法)

来源:互联网 发布:office mac 2011 登录 编辑:程序博客网 时间:2024/06/05 22:46

描述下问题:求直线ax+bx+c=0上有多少个整点(x,y)满足x-[x1,x2],y-[y1,y2]。
根据数论的理论,要求ax+bx=-c的一个解(x0,y0),先要求出ax+bx=gcd(a,b)的一个解(x1,y1),然后根据gcd(a,b),记作g,判断它与c的倍数,由(x1,y1)来得到解(x0,y0),扩展欧几里得算法就是求ax+bx=gcd(a,b)的一个解的。
边界条件和欧几里得算法一样是b=0时,a=gcd,且此时x=1,y=0.
这里参考一篇博客给出想x,y状态更新的方程的解释:
参考的博客:http://blog.csdn.net/zhjchengfeng5/article/details/7786595
欧几里德算法停止的状态是: a= gcd , b = 0 ,那么,这是否能给我们求解 x y 提供一种思路呢?因为,这时候,只要 a = gcd 的系数是 1 ,那么只要 b 的系数是 0 或者其他值(无所谓是多少,反正任何数乘以 0 都等于 0 但是a 的系数一定要是 1),这时,我们就会有: a*1 + b*0 = gcd
当然这是最终状态,但是我们是否可以从最终状态反推到最初的状态呢?
假设当前我们要处理的是求出 a 和 b的最大公约数,并求出 x 和 y 使得 a*x + b*y= gcd ,而我们已经求出了下一个状态:b 和 a%b 的最大公约数,并且求出了一组x1 和y1 使得: b*x1 + (a%b)*y1 = gcd , 那么这两个相邻的状态之间是否存在一种关系呢?
我们知道: a%b = a - (a/b)*b(这里的 “/” 指的是整除,例如 5/2=2 , 1/3=0),那么,我们可以进一步得到:
gcd = b*x1 + (a-(a/b)*b)*y1
= b*x1 + a*y1 – (a/b)*b*y1
= a*y1 + b*(x1 – a/b*y1)
对比之前我们的状态:求一组 x 和 y 使得:a*x + b*y = gcd ,是否发现了什么?
这里:
x = y1
y = x1 – a/b*y1
所以每次进行状态改变就把y1赋给x,x1传入到y,然后更新y-=a/b*x
代码如下:

void gcd(int a, int b, int&d,int &x, int &y){    if (!b){        d = a; x = 1; y = 0;    }    else {        gcd(b, a%b, d,y, x);        y -= x*(a / b);    }}

根据得到的x,y的值和d(gcd的值)可以得到ax+bx=-c的解(x,y),然后根据公式(x+kb’,y-ka’)得到该直线的通解,k为整数,即可解决直线上的点的问题。其中a’=a/gcd,b’=b/gcd
a*x + b*y = c 有解的充要条件: c % gcd(a , b) == 0也就是c是gcd(a,b)的整数倍_


扩展欧几里德有什么用处呢?
求解形如 a*x +b*y = c 的通解,但是一般没有谁会无聊到让你写出一串通解出来,都是让你在通解中选出一些特殊的解,比如一个数对于另一个数的乘法逆元
什么叫乘法逆元?
a*x≡1(mod m),注意这里的记号≡,表示a与1关于模m同余,amod m=1mod m,即a-1是m的整数倍
这里,我们称 x 是 a 关于 m 的乘法逆元
这怎么求?可以等价于这样的表达式: a*x + m*y = 1
由之前的结论,1必须是gcd(a,m)的整数倍,所以呀,a,m必须是互质的。
接着乘法逆元讲,一般,我们能够找到无数组解满足条件,但是一般是让你求解出最小的那组解,怎么做?我们求解出来了一个特殊的解 x0 那么,我们用 x0 % m其实就得到了最小的解了。为什么?
因为x=x0+kb’=x+k*m/gcd=x+k*m
那么,也就是说, a 关于 m 的逆元是一个关于 m 同余的,那么根据最小整数原理,一定存在一个最小的正整数,它是 a 关于m 的逆元,而最小的肯定是在(0 , m)之间的,而且只有一个。由上式可以知道x0%m就是我们要求的最小的同余解x

有时候我们得到的特解 x0 是一个负数,还有的时候我们的 m 也是一个负数这怎么办?因为我们要求的是正整数,所以

当 m 是负数的时候,我们取 m 的绝对值就行了,当 x0 是负数的时候,他模上 m 的结果仍然是负数(在计算机计算的结果上是这样的,虽然定义的时候不是这样的),这时候,我们仍然让 x0 对abs(m) 取模然后结果再加上abs(m) 就行了,于是,我们不难写出下面的代码求解一个数 a 对于另一个数 m 的乘法逆元:

int calc_inversemin(int a, int b){    int x, y, d;    gcd(a, -b, d, x, y);    if (abs(1%d))return -1;    x *= 1 / d;    int m = abs(b);    int ans = x%m;    if (ans <= 0) ans += m;    return ans;}

注意判断是否存在逆元时,要在1%d处加绝对值,d代表最大公约数,但是最大公约数可能是负数,比输入参数是2,3的时候,转换成将2x≡1(mod 3)转换成方程为2x-3y=1,那么我们所求的是2和-3的最大公约数,所以第二个参数是-b,且2,-3的最大公约是-1,也说明是有逆元的。
但是为了记忆的简便,可以按下面模板进行记忆,即使对模线性方程组没考虑公约数的正负问题(不知道为什么。。尝试了几个例子,发现针对模线性方程而言,我要使用的是x,而x的系数没带负号,所以对结果按gcd(a,n)和gcd(a,-n)两者是一样的),也可以得到正确的答案:

typedef long long LL;LL cal(LL a,LL b,LL c)  {      LL x,y,gcd_num;      LL gcd(a,b,gcd_num,x,y);      if(c%gcd_num!=0) return -1;      x*=c/gcd_num;      b/=gcd_num;      if(b<0) b=-b;      LL ans=x%b;      if(ans<=0) ans+=b;      return ans;  }