POJ 2115 线性模方程 (归根到底还是扩展欧几里得) 转载请标明出处http://www.cnblogs.com/MasNoon/p/5733264.html
题意:给你4个数A,B,C,K,求出(i=A;i!=B;i+=c)这个循环在int k 位里能执行几次。int K位其实就是在2^k这个范围内执行几次,超过2^k,又会重新从0计数。
根据题意列方程:A+X*C=B->x*C=(B-A)%2^K就是一个线性模方程。
那怎么解线性模方程呢?先给出几个定理
定理1:方程ax=b(mod n)对于未知量x有解,当且仅当gcd(a, n)|b(就是gcd(a,b)能整除b,即b%gcd(a,b)=0)
定理2:方程ax=b(mod n)或者对模n有d个不同的解其中d=gcd(a, n),或者无解.
定理3:假设方程ax=b(mod n)有解(即有d|b,其中d=gcd(a, n)),x0是该方程的任意一个解,则该方程对模n恰有d个不同的解,分别为:xi=x0+i(n/d)(i = 1, 2, …, d-1)
再来看方程x*C=(B-A)%2^K,由定理1可知,如果(B-A)%gcd(C,2^K)不等于0,则方程无解,即输出"FOREVER".
如果(B-A)%gcd(C,2^K)等于0,那方程可以写成X*C+n*2^K=B-A,学过扩展欧几里得的孩子都应该觉得很眼熟吧。下面讲讲扩展欧几里得,学过的同学就可以跳着看了。
也是先给出定理
定理:对于不完全为0的非负整数a,b,gcd(a,b)表示a,b的最大公约数d,必然存在整数对x,y,使得gcd(a,b)=d=ax+by
对于gcd(a,b) = d,对(a, b)用欧几里德辗转相除会最终得到(d, 0)。此时对于把a =d, b = 0 代入a*x + b*y = d,显然x = 1,y可以为任意值。 我们可以用a = d, b = 0的情况逆推出来任何gcd(a, b) = d 满足a*x + b*y = d的解。如果x0,y0是b*x + (a%b)*y = d 的解,那么对于a*x + b*y = d的解呢?
b*x + (a%b)*y = d → b*x + (a - [a/b]*b)*y = d → a*y + b*(x - [a/b]*y) = d 所以a*x + b*y = d的解 x1 = y0,y1= x0 - [a/b]*y0;
现在知道(B-A)%gcd(C,2^K)等于0,令L=(B-A)/gcd(C,2^K)
x1*C+y1*2^k=gcd(C,2^k) ①
方程①*L等于方程X*C+n*2^K=B-A
即x=x1*L(其实从这里可以推出线性模方程的定理一,如果(B-A)%gcd(C,2^K)不等于0,L肯定是个小数,而我们知道X是整数,x1也是整数,L*x1肯定不会等于X.所以无解)
又从线性模方程的定理三可知,xi=x0+i(n/d)(i = 1, 2, …, d-1)(n/d就是这里的2^k/gcd(C,2^k))
我们要求的肯定要是最小解,(废话:第一次相等的时候循环就停止了),xi=x0+i(n/d)->x0=xi%(n/d)这样最小解就出来了。
所以答案是x1*L%(2^k/gcd(C,2^k)).
下面贴代码
#include<stdio.h>
#include<cstring>
#include<iostream>
using namespace std;
__int64 ex_Euc( __int64 a, __int64 b, __int64 &x, __int64 &y)//扩展欧几里得,递归
{
if(b==0)
{
x=1,y=0;
return a;
}
__int64 d=ex_Euc(b,a%b,x,y);
__int64 t=x;
x=y;
y=t-a/b*y;
return d;
}
__int64 quickpow( __int64 a, __int64 b)
{
__int64 res=1;
while(b!=0)
{
if(b&1)
{
res*=a;
}
a=a*a;
b>>=1;
}
return res;
}
int main(void)
{
__int64 a,b,c,k;
while ( scanf("%I64d%I64d%I64d%I64d",&a,&b,&c,&k)!=EOF)
{
if(a==0&&b==0&&c==0&&k==0)
{
break;
}
__int64 m=quickpow(2,k);
__int64 x;
__int64 y;
__int64 g=ex_Euc(c,m,x,y);
__int64 n=b-a;
if(n%g!=0)
{
puts("FOREVER");本身就有空行不用再\n
continue;
}
x=(x*(n/g))%m;//(x*n%(m/g))
x=(x%(m/g)+m/g)%(m/g);//避免有负数,mod 运算,(a+b)%b=a%b
printf("%I64d\n",x);
}
return 0;
}