求解线性同余方程--扩展欧几里得

本文主要是介绍求解线性同余方程--扩展欧几里得，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

资料来源：https://blog.csdn.net/

//求解ax=b(mod m) 返回0为无解，否则返回gcd(a,m)个mod m意义下的解,用X[]存
int mod(int a, int b, int m)
{return equation(a, m, b);
}

先看一道题目：

poj 2115 C Looooops 【扩展欧几里得】

Description:

A Compiler Mystery: We are given a C-language style for loop of type

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

I.e., a loop which starts by setting variable to value A and while variable is not equal to B, repeats statement followed by increasing the variable by C. We want to know how many times does the statement get executed for particular values of A, B and C, assuming that all arithmetics is calculated in a k-bit unsigned integer type (with values 0 <= x < 2k) modulo 2k.

Input：

The input consists of several instances. Each instance is described by a single line with four integers A, B, C, k separated by a single space. The integer k (1 <= k <= 32) is the number of bits of the control variable of the loop and A, B, C (0 <= A, B, C < 2k) are the parameters of the loop.

The input is finished by a line containing four zeros.

Output：

The output consists of several lines corresponding to the instances on the input. The i-th line contains either the number of executions of the statement in the i-th instance (a single integer number) or the word FOREVER if the loop does not terminate.

Sample Input

Sample Output

0
2
32766
FOREVER

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

资料来源1：https://blog.csdn.net/

题意1：求解最小的y使得 A + C*y = B (mod 1<<k)。

思路1：化简等式 (1<<k)*x + y*C = B-A。当且仅当(B-A) % gcd(1<<k, C) == 0才有解。

令temp = (B-A)/gcd(1<<k, C);

化简——(1<<k)*x/temp + y*C/temp = gcd(1<<k, C);

求出x和y，y *= temp。

由扩展欧几里得定理三——若gcd(a, b) = d，则方程ax ≡ c (mod b)在[0, b/d - 1]上有唯一解。

求解最小的y时，只需对(1<<k) / gcd(1<<k, C)取余即可。

AC代码1：

#include <cstdio>
#include <cstring>
#include <algorithm>
#define LL long long
using namespace std;LL A, B, C, k;
void exgcd(LL a, LL b, LL &d, LL &x, LL &y)
{if(b == 0) {d = a, x = 1, y = 0;}else{exgcd(b, a%b, d, y, x);y -= x * (a / b);}
}
void solve()
{B -= A;LL d, x, y;LL a = (1LL<<k);exgcd(a, C, d, x, y);if(B % d)printf("FOREVER\n");else{x *= (B / d); y *= (B / d);y = (y % (a / d) + a / d) % (a / d);printf("%lld\n", y);}
}int main()
{while(scanf("%lld%lld%lld%lld", &A, &B, &C, &k) != EOF){if(A == 0 && B == 0 && C == 0 && k == 0)break;solve();}return 0;
}

资料来源2：https://www.cnblogs.com/

本题和poj1061青蛙问题同属一类，都运用到扩展欧几里德算法，可以参考poj1061，解题思路步骤基本都一样。
一，题意：
　　对于for(i=A ; i!=B ;i+=C)循环语句，问在k位存储系统中循环几次才会结束。
　　　　比如：当k=4时，存储的数 i 在0-15之间循环。(本题默认为无符号)
　　若在有限次内结束，则输出循环次数。
　　否则输出死循环。
二，思路：
本题利用扩展欧几里德算法求线性同余方程，设循环次数为 x ,则解方程 (A + C*x) % 2^k = B ;求出最小正整数 x。
　　1，化简方程化为求线性同余方程标准式 ax ≡ b (mod n);
　　2，扩展欧几里德算法求解线性同余方程 C*x ≡ B-A (mod 2^k);
　　3，求出最小非负整数解。
三，步骤：
　　1，化简：(A + C*x) mod 2^K = B --> C*x mod 2^k = B-A --> C*x ≡ B-A (mod 2^k);
　　2，求线性同余方程 C*x ≡ B-A (mod 2^k) , 就相当于求二元一次方程 C*x + 2^k * y = B-A
　　　　i,代入扩展欧几里德算法，求解方程 C*x + 2^k * y = gcd(C , 2^k) ;
　　　 ii,利用方程 C*x + 2^k * y = gcd(C , 2^k)的解 x0 以及公式 x1 = x0 * c/d 求出原方程 a*x + b*y = c 的解 x1 ;前提是：d|c (c 能被 d 整除);
　　3，利用周期性变化求最小的非负整数解公式: x1 = (x1 % (b/d) + (b/d) ) % (b/d);

若方程的C*x + 2^k * y = B-A 的一组整数解为(x1 , y1)，则它的任意整数解为(x1 + k * (b/d) , y1 - k * (a/d) ) ( k取任意整数 ), T = b/d就为 x1 增长的周期

　　　　i，若x1为负值，取最大的非正值：x1 = x1 % T ; 若x1为正值，以下两步无影响;
　　　　ii，取正：x1 = x1 + T ;
　　　　iii，防止 i 中的 x1=0 即 ii 中的 x1=T ：x1 = x1 % T ;

代码如下：

#include<iostream>
using namespace std;void exgcd(long long a,long long b,long long& d,long long& x,long long& y){//int& a 是定义一个存放整形变量a的地址if(!b){ d=a ; x=1 ; y=0; }              // d用来存储gcd(a,b)的值 else { exgcd(b , a%b , d , y , x); y -= x* (a/b); }
}int main(){long long  A,B,C,d,x,y,T;int k ; while(cin>>A>>B>>C>>k){if(A==0&&B==0&&C==0&&k==0)break;long long n = 1LL<<k; //n = 1 * 2^k ;注意此处,若为__int64,则应该是n = (__int64)1 << k;exgcd(C,n,d,x,y);           if( (B-A) % d != 0 ){cout<<"FOREVER\n";}else {x = x * (B-A) / d ;  T = n / d;x = ( x%T + T ) % T ;cout<<x<<endl;}} return 0; 
}

AC代码3：

就是求一个同余方程组的最小解

代码如下：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;typedef long long LL;LL ex_gcd(LL a,LL b,LL &x,LL &y){if(!b){x=1,y=0;return a;}else{LL gcd = ex_gcd(b,a%b,x,y);LL tmp = x;x = y;y = tmp-a/b*y;return gcd;}
}int n;
LL r1,r2,a1,a2,x,y;LL solve(){LL r1,r2,a1,a2,x,y;bool tag = 1;cin>>a1>>r1;for(int i=1;i<n;i++){cin>>a2>>r2;LL a = a1,b=a2,c = r2-r1;LL d = ex_gcd(a,b,x,y);if(c%d) tag = 0;LL t = b/d;x = (x*(c/d)%t+t)%t;r1 = a1*x+r1;a1 = a1*(a2/d);}if(tag) return r1;else return -1;
}int main()
{while(~scanf("%d",&n)){printf("%lld\n",solve());}return 0;
}

AC代码4：

/*
题意：
(A+Cx)%2^k=B ==> A+Cx = 2^k*y+B
==> Cx-2^k*y=B-A;
形如 ax+by=c
a=C
b=2^k
c=B-A;*/#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
using namespace std;__int64 f[33];
void make_hxl()
{int i;f[0]=1;f[1]=2;for(i=2;i<=32;i++)f[i]=f[i-1]*2;
}__int64 Ex_gcd(__int64 a,__int64 b,__int64 &x,__int64 &y)
{if(b==0){x=1;y=0;return a;}__int64 g=Ex_gcd(b,a%b,x,y);__int64 hxl=x-(a/b)*y;x=y;y=hxl;return g;
}void make_ini(__int64 A,__int64 B,__int64 C,__int64 k)
{__int64 c,a,b,x,y,d;a=C;b=f[k];c=B-A;d=Ex_gcd(a,b,x,y);if(c%d!=0){printf("FOREVER\n");return;}x=c/d*x;b=b/d;x=(x%b+b)%b;printf("%I64d\n",x);
}
int main()
{__int64 A,B,C,k;make_hxl();while(scanf("%I64d%I64d%I64d%I64d",&A,&B,&C,&k)>0){if(A==0 && B==0 && C==0 &&k==0)break;make_ini(A,B,C,k);}return 0;
}

再来看一道题目：

青蛙的约会

Description

两只青蛙在网上相识了，它们聊得很开心，于是觉得很有必要见一面。它们很高兴地发现它们住在同一条纬度线上，于是它们约定各自朝西跳，直到碰面为止。可是它们出发之前忘记了一件很重要的事情，既没有问清楚对方的特征，也没有约定见面的具体位置。不过青蛙们都是很乐观的，它们觉得只要一直朝着某个方向跳下去，总能碰到对方的。但是除非这两只青蛙在同一时间跳到同一点上，不然是永远都不可能碰面的。为了帮助这两只乐观的青蛙，你被要求写一个程序来判断这两只青蛙是否能够碰面，会在什么时候碰面。

我们把这两只青蛙分别叫做青蛙A和青蛙B，并且规定纬度线上东经0度处为原点，由东往西为正方向，单位长度1米，这样我们就得到了一条首尾相接的数轴。设青蛙A的出发点坐标是x，青蛙B的出发点坐标是y。青蛙A一次能跳m米，青蛙B一次能跳n米，两只青蛙跳一次所花费的时间相同。纬度线总长L米。现在要你求出它们跳了几次以后才会碰面。

Input

输入只包括一行5个整数x，y，m，n，L，其中x≠y < 2000000000，0 < m、n < 2000000000，0 < L < 2100000000。

Output

输出碰面所需要的跳跃次数，如果永远不可能碰面则输出一行"Impossible"

Sample Input

1 2 3 4 5

Sample Output

资料来源：https://mp.csdn.net/

AC代码1：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define mem(a,x) memset(a,x,sizeof(a))
using namespace std;
typedef long long ll;
ll e_gcd(ll a,ll b,ll &x,ll &y)
{if (b == 0){x = 1,y = 0;return a;}ll ans = e_gcd(b,a%b,x,y);ll tmp = x;x = y;y = tmp - a/b*y;return ans;
}
ll cal(ll a,ll b,ll c)//求最小的x使ax+by=c
{ll x,y;ll gcd = e_gcd(a,b,x,y);if (c%gcd != 0) return -1;x *= c/gcd;b/=gcd;if (b < 0) b = -b;ll ans = x%b;if (ans <= 0) ans += b;return ans;
}
int main()
{ll xa,xb,va,vb,L;while (~scanf("%lld %lld %lld %lld %lld",&xa,&xb,&va,&vb,&L)){ll ans = cal(vb-va,L,xa-xb);if (ans == -1) puts("Impossible");else printf("%lld\n",ans);}return 0;
}

资料来源：https://blog.csdn.net/

AC代码2：

青蛙的约会题解：

假设走了t次相遇，则有等式(x+mt)-(y+nt)=pL成立,等价于求解同余方程(n-m)t≡(x-y) (mod L)的最小整数解

（a）对于一般同余方程ax=d mod b,方程有解，则有（a,d）| b ,所以问题第一步判断解的情况

（b）有（n-m）t+pL=x-y,t、p均为未知变量，将问题转化为求解ax+by=d的最小整数x，扩展欧几里得算法：

briefly:扩展欧几里得算法是辗转相除法求gcd的拓展，表现在ax+by=gcd(a,b)，函数extended_gcd()不仅能返回gcd(a,b),还能求出gcd的线性系数x,y,具体的操作步骤如下：

①首先化简，得到新的ax+by=d,注意此时（a,b）=1

②先求ax+by=1的解x0、y0（解具有唯一性），利用扩展欧几里得算法得到唯一解x0,则ax+by=d的解x=d*x0

③通解X=x+b*k(k为整数)

（c）通过(b)实际上可以得到同余方程的通解，但是题目要求最小整数解，利用min=(X%b+b)%b,X取正取负均满足最小，问题得解

以CX=B-A(mod M)为例：

#include 
#include 
#include 
#include 
#include 
#include 
#include
using namespace std;
typedef long long ll;
const int maxn=1000000+10;
ll x,y;
ll extended_gcd(ll a,ll b,ll &x,ll &y)//扩展欧几里得算法求ax+by=gcd(a,b)的线性系数x,y，返回gcd(a,b)
{if(b==0) {x=1;y=0;return a;}ll d=extended_gcd(b,a%b,y,x);y-=a/b*x;return d;
}
int main()
{// freopen("input.txt","r",stdin);ll a,b,m,A,B,C,k;while(cin>>A>>B>>C>>k&&(A||B||C||k)){a=C;b=(ll)1<

资料来源：http://www.cnblogs.com/

AC代码3：

/*
题意：两只青蛙同方向跳给出初始位置，和每一次跳的距离和总的长度
问第几天跳到一起，约会了.两只青蛙在同一个地方，满足<------------------x-------------y---------m*T-p*L = n*T-(x-y);   向西为正
转化成  T*(m-n)-p*L =y-x;
==>     ax+by=c;a=m-n;b=-L;c=y-x;
*/#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
using namespace std;__int64 Ex_GCD(__int64 a,__int64 b,__int64 &x,__int64 &y)
{if(b==0){x=1;y=0;return a;}__int64 g=Ex_GCD(b,a%b,x,y);__int64 hxl=x-(a/b)*y;x=y;y=hxl;return g;
}int main()
{__int64 x,y,m,n,l,a,b,c,k,x1,y2;while(scanf("%I64d%I64d%I64d%I64d%I64d",&x,&y,&m,&n,&l)>0){a=m-n;b=-l;c=y-x;k=Ex_GCD(a,b,x1,y2);if(c%k!=0){printf("Impossible\n");continue;}b=b/k;//   b/gcd(a,b);x1=x1*(c/k); // x=x*c/gcd(a,b);x1=x1%(b);while(x1<0) x1=x1+b;printf("%I64d\n",x1);}return 0;
}