codeforces724C. Ray Tracing（扩展欧几里得详解）

传送门：C. Ray Tracing。

题意：一个矩形内有k个点，t=0从原点45度射出一束光线，遇到矩形的边就折射（入射角等于出射角），遇到顶点就停止，每秒前进根二单位长度，问矩形内的点第一次遇到光线是第几秒。

解法：首先把矩形展开成一个无边界平面，这时我们可以发现当x==y==n*m/gcd(n,m)时光线停止，设原坐标为a，b，则当展开后横坐标为a,2*n±a……2*x*n±a。同理纵坐标为2*y*m±b（这一点自己画图扩展一下便知）。因为光线45度角射出去，可知当横坐标等于纵坐标

即2*x*n±a==2*y*m±b（x、y<n*m/gcd(n,m)）的时候点会受到光照。因为是第一次遇到这里的x、y应为最小。即所求问题转化成2*x*n±a==2*y*m±b，即2*x*n - 2*y*m == ±(a ± b)，使2*x*n和2*y*m最小，即求使等式成立的最小x，y。

就可以转化成典型的扩展欧几里得算法可以解决的问题ax+by=c，已知a和b，求使等式成立的x，y（整数）。

扩展欧几里得算法解决的最基本问题是gcd（a，b）=ax+by。证明转载至下：

对于不完全为 0 的非负整数 a，b，gcd（a，b）表示 a，b 的最大公约数，必然

存在整数对 x，y ，使得 gcd（a，b）=ax+by。

求解 x，y的方法的理解

设 a>b。

1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0；

2，a>b>0 时

设 ax₁+ by₁= gcd(a,b);

bx₂+ (a mod b)y₂= gcd(b,a mod b);

根据朴素的欧几里德原理有 gcd(a,b) = gcd(b,a mod b);

则:ax₁+ by₁= bx₂+ (a mod b)y₂;

即:ax₁+ by₁= bx₂+ (a - [a / b] * b)y₂=ay₂+ bx₂- [a / b] * by₂;

也就是ax₁+ by1 == ay₂+ b(x₂- [a / b] *y₂);

根据恒等定理得：x₁=y₂; y₁=x₂- [a / b] *y₂;

这样我们就得到了求解 x₁,y₁ 的方法：x₁，y₁ 的值基于 x₂，y₂.

上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以结束。

c++实现：

[cpp] view plain copy

1. int exgcd(int a,int b,int &x,int &y){  
2.     if (b==0){  
3.         x=1,y=0;  
4.         return a;  
5.     }  
6.     int q=gcd(b,a%b,y,x);  
7.     y-=a/b*x;  
8.     return q;  
9. }  

或者

[cpp] view plain copy

1. int exGcd(int a,int b,int &x,int &y)  
2. {  
3.     if(b==0)  
4.     {  
5.         x=1;y=0;  
6.         return a;  
7.     }  
8.     int r=exGcd(b,a%b,x,y);  
9.     int t=x;x=y;y=t-a/b*y;  
10.     return r;  
11. }  

这样求出来的x0和y0是gcd（a，b）=ax+by的一个特解，通解可以用下面的方法表示：

x = x0 + (b/gcd)*t
y = y0 – (a/gcd)*t   证明如下：

首先把 a/gcd(a,b)化为a'，b/gcd(a,b)化为b'，这样ax+by=1，gcd(a,b)=1（注意这里a，b实际上写作a'和b'，为简单起见省掉撇号）
由于
ax+by = ax0+by0
所以
a(x-x0) = -b(y-y0) = P
所以
a|P，b|P（“|”表示整除关系）
因为a，b互素，所以
ab|P
即
P=kab
所以
x-x0=kb，y-y0=-ka
整理得到
x=x0+kb，y=y0-ka

再由gcd（a，b）=ax+by的通解可以推导出更为一般的ax+by=c的通解，此处若c mod gcd（a，b）==0则有整数解，否则无解，证明如下：

由于我们知道，存在一组x与y使得a*x+b*y=gcd(a,b)。
将等式两边同时乘以整数k，即a*x*k+b*y*k=gcd(a,b)*k。如果c mod gcd(a,b)=f，则0<=f<gcd(a,b)。
那么可以令c=gcd(a,b)*k+f。这样一来，就有a*x*k+b*y*k+f=c。
若f  0，由于f<gcd(a,b)<=a<=b（假设a<=b），所以不存在f=a*m（m为整数），也就不存在a*(x*k+m)+b*y*k=c。也就是说，不存在a*x+b*y=c的整数解x与y。
所以f=0，即只有当c mod gcd(a,b)=0时，a*x+b*y=c有正整数解。得证。

而ax+by=c的整数解，只需将x * a+y * b = Gcd(a, b)的每个解乘上 c/Gcd(a, b) 即可，但是所得解并不是该方程的所有解，找其所有解的方法如下：
在找到x * a+y * b = Gcd(a, b)的一组解x0,y0后，可以
得到x * a+y * b = c的一组解x1 = x0*(c/Gcd(a,b)),y1 = y0*(c/Gcd(a,b))，则x * a+y * b = c的其他整数解满足：
x = x1 + b/Gcd(a, b) * t
y = y1 - a/Gcd(a, b) * t(其中t为任意整数)
x 、y就是x * a+y * b = c的所有整数解。

至于这个通解为什么是这样的证明，我也没搞太懂，个人认为可以类比上面求gcd（a，b）=ax+by的通解证明，还望各路数论大神指点。

扩展欧几里得定理常用于求解中国剩余定理相关问题和求乘法逆元的问题：

中国剩余定理：即求同余方程组的解，可以做如下转化：

求解方程组c≡x(mod m)，c≡b(mod n)的c，即

求解同余方程a≡c(mod m)，b≡c(mod n)

将其化为a=c+mx，b=c+ny，作差得

mx-ny=a-b，从而化为线性不定方程利用扩展欧几里得定理解决。

模m乘法逆元：
定义：对于整数a，m，如果存在整数b，满足ab ≡ 1(mod m)，则说，b是a的模m乘法逆元。
定理：a存在模m的乘法逆元的充要条件是gcd(a,m) = 1

证明：充分性：
因为
gcd(a,m) = 1
根据欧拉定理，有
a^φ(m) ≡ 1(mod m)
因此
a * a^(φ(m)-1) mod m = 1
所以存在a的模m乘法逆元，即a^(φ(m)-1)

必要性：
假设存在a模m的乘法逆元为b，则
ab ≡ 1 (mod m)
所以
ab = km +1
所以
1 = ab - km
由欧几里得定理，有
gcd(a,m) = 1
由定理知：
对于ax + my = 1，可以看出x是a模m的乘法逆元，y是m模a的乘法逆元。
反过来，要计算a模m的乘法逆元，就相当于求ax + my = 1的x的最小正整数解，从而化为线性不定方程解决。

下面是如何求解最小正整数解（转）：一般，我们能够找到无数组解满足条件，但是一般是让你求解出最小的那组解，怎么做？我们求解出来了一个特殊的解 x’ 那么，我们用 x' % m其实就得到了最小的解了（更广泛情况应为x0%（m/gcd））。为什么？
可以这样思考：
    x 的通解不是 x' + m*t （t为任意整数）吗？（更广泛情况应为x0+（m/gcd）*t）(个人理解：x的通解可以由任何特解表示出来，因此也可以用最小正整数特解表示出来，假设最小正整数解为x0，则x=x0+m*t，那么x’=x0+m*t)
    那么，也就是说， a 关于 m 的所有逆元是关于 m 同余的，那么根据最小整数原理，一定存在一个最小的正整数，它是 a 关于m 的逆元，而最小的逆元肯定是在（0 , m）之间的，而且只有一个，这就好解释了（更广泛情况应为a 关于 m 的所有逆元是关于 m/gcd 同余的，因为x%（m/gcd）==x0，x为a的任意逆元，x0为最小正整数解，x0在（0，m/gcd）之间，括号外面最后一句话解释了x0的存在性且唯一性）
。

    可能有人注意到了，这里，我写通解的时候并不是 x0 + (m/gcd)*t ，但是想想一下就明白了，gcd = 1，所以写了跟没写是一样的，但是，由于问题的特殊性，有时候我们得到的特解 x0 是一个负数，还有的时候我们的 m 也是一个负数这怎么办？
    当 m 是负数的时候，我们取 m 的绝对值就行了，当 x0 是负数的时候，他模上 m 的结果仍然是负数（在计算机计算的结果上是这样的，虽然定义的时候不是这样的），这时候，我们仍然让 x0 对abs(m) 取模，然后结果再加上abs(m) 就行了，于是，我们不难写出下面的代码求解一个数 a 对于另一个数 m 的乘法逆元：

[cpp] view plain copy

1. LL cal(LL a,LL m)    
2. {    
3.     LL x,y;    
4.     LL gcd=exgcd(a,m,x,y);    
5.     if(1%gcd!=0) return -1;    
6.     x*=1/gcd;     
7.     if(m<0) m=-m;    
8.     LL ans=x%m;    
9.     if(ans<=0) ans+=m;    
10.     return ans;    
11. }    

这个地方可能会有人和我刚看到这段代码时一样的疑问，既然明知道要判断是否gcd==1了，那干嘛要写1%gcd！=0啊，直接写gcd！=1不就好了，其实这是因为这个子函数是一个模板，一般我们要求的等式为a*x+b*y=c（此处b就相当于上面代码中m），c==1时我们可以像上面这么写（当然可以更简洁），如果c为任意数，代码就变成了如下的样子：

[cpp] view plain copy

1. LL equation(LL a, LL b, LL c, LL &x, LL &y)  //a*x+b*y==c  
2. {    
3.     LL g = extend_Euclid(a, b, x, y);   
4. //返回值为gcd（a，b），不明白请向上翻看exgcd代码，再不明白就去搜辗转相除法（欧几里得）的c++实现吧。。  
5.      if(c % g) return -1; //判断是否有整数解，证明向上翻   
6.     LL ran = b / g; //求"m%gcd"(见上面红色字体)，这里的b相当于上面m，因此ran就是"m%gcd"   
7.     x *= c/g;  //将x扩大c/gcd（a，b）倍，至于为什么要扩大，依然向上翻至绿色部分   
8.    if(ran < 0) ran = -ran;  //将m%gcd变为正数  
9.     x = (x%ran + ran) % ran;  //让x对ran取余求最小正整数解，即上面红色字体x对（m%gcd）取余，这一步操作还保证了x为正数  
10.     return 0;  }    

这才是模板最基本的样子。

下面附上本题代码：

[cpp] view plain copy

1. #include <iostream>  
2. #include <cstdio>  
3. #include <cstring>  
4. #include <math.h>  
5. #include <algorithm>  
6. #define ll long long  
7. using namespace std;  
8. int gcd(int a,int b)//辗转相除法  
9. {  
10.     return b?gcd(b,a%b):a;  
11. }  
12. ll n,m;  
13. ll M;  
14. ll extend_Euclid(ll a,ll b,ll &x,ll &y)//扩展欧几里得算法  
15. {  
16.     if(b==0)  
17.     {  
18.         x=1;y=0;  
19.         return a;  
20.     }  
21.     ll ans=extend_Euclid(b,a%b,y,x);  
22.     y-=(a/b)*x;  
23.     return ans;  
24. }  
25. ll equation(ll a,ll b,ll &x,ll &y,ll c)//这个函数见上面详解  
26. {  
27.     ll gcd=extend_Euclid(a,b,x,y);  
28.     if(c%gcd) return -1;  
29.     ll k=abs(b/gcd);  
30.     x*=c/gcd;  
31.     x=(x%k+k)%k;  
32.     return 0;  
33. }  
34. ll judge(ll x,ll y)  
35. {  
36.     ll p,q;  
37.     if(equation(2*n,-2*m,p,q,-x+y)==-1)//判断点是否会被照射到  
38.     return M+1;  
39.     ll t=2*p*n+x;  
40.     if(t<0||t>M)//判断求出来的时间是否在范围内，即点应该在光线停止传播之前被照射到  
41.     return M+1;  
42.     return t;  
43. }  
44. ll solve(ll x,ll y)  
45. {  
46.     ll g=gcd(n,m);  
47.     M=1ll*m/g*n;//此处M为光线结束传播的时间  
48.     ll ans=M+1;  
49.     ans=min(ans,judge(x,y));//枚举2*x*n - 2*y*m == ±(a ± b)中等号右边情况  
50.     ans=min(ans,judge(x,-y));  
51.     ans=min(ans,judge(-x,y));  
52.     ans=min(ans,judge(-x,-y));  
53.     if(ans==M+1)  
54.     return -1;  
55.     return ans;  
56. }  
57. int main()  
58. {  
59.     int k;  
60.     scanf("%d%d%d",&n,&m,&k);  
61.     ll x,y;  
62.     for(int i=0;i<k;i++)  
63.     {  
64.         scanf("%lld%lld",&x,&y);  
65.         printf("%lld\n",solve(x,y));  
66.     }  
67. return 0;  
68. }  

以上内容部分转载于：http://blog.csdn.net/synapse7/article/details/9901195；百度百科；http://blog.csdn.net/zhjchengfeng5/article/details/7786595