Number Theory -- Solve problems on reminders
其实这个东西我也一直很晕，只知道算法，具体的东西从来没仔细推过......（太懒了）
这两天才仔细把这个问题研究了一下，为了防止以后忘掉，记下来好了。
Part I 最大公约数和辗转相除法
所有的问题都是从求两个数的最大公约数开始的。
两个整数的公约数定义如下：d是a , b的公约数，当且仅当d | a && d | b。
a , b的最大公约数就是a和b公约数里最大的一个，记作gcd(a , b)或(a , b)。
最大公约数的一种非常流行的求法就是欧几里德辗转相除法。
算法基于这样一个事实
gcd(a , b) = gcd(b , a % b)
具体一些来说，我们要求(a0 , b0)的最大公约数，那么设
a0 = q0 * b0 + b1 ... (1)，其中q0是a0除以b0的商，那么0 <= b1 < b0。
设d = (a0 , b0)，那么观察等式(1)，d | a0 , d | q0 * b0 ，所以d | b1，所以gcd(a0 , b0) = gcd(b0 , b1)
同样的，我们可以得到
b0 = q1 * b1 + b2 ... (2)
b1 = q2 * b2 + b3 ... (3)
...
bk = qk+1 * bk+1 ... (k+1)
不断的迭代下去，一定存在一个等式，最后一项为0，不妨设是第k + 1个等式，为什么呢？
0 = bk+1 < bk < bk-1 < ... < b1 < b0，由于小于b0的正整数是有限的，所以必存在bk+1 = 0。
那么(a0 , b0) = (b0 , b1) = ... = (bk , bk+1) = bk+1
这样就得到了gcd(a0 , b0)。
写个简单的程序
int gcd(int a , int b)
{
if (!b)
return a ;
else
return gcd(b , a % b)
}
写成迭代的形式，效率会更高一些。
Part II ax + by = (a , b)
一个新的问题：若(a , b) = d，必存在整数x , y满足ax + by = d，并求出x , y。
可以仿照辗转相除法迭代求x , y出来。
现在我们在迭代的过程中，知道了如下事实，我们设返回的最大公约数为d：
B = q * C + D ... (1) , x * B + y * C = d ... (2)
我们要求的就是x' , y'满足
A = q' * B + C ... (3) , x' * A + y' * B = d ... (4)
由(2)和(4) 有：
x' * A + y' * B = x * B + y * C = d ... (5)
将(3)代入，得到：
x' * (q' * B + C) + y' * B = x * B + y * C
x' * q' * B + x' * C + y' * B = x * B + y * C
现在这样构造，令：
x' = y
则有
y' * B = x * B - x' * q' * B
y' = x - x' * q' = x - y * q'
这样我们得出新的x' , y'
x' = y
y' = x - y * [A / B]
这样我们在迭代中就求出了x和y。
int extended_gcd(int a , int b , int& x , int& y)
{
if (!b)
{ x = 1 ; y = 0 ; return a ; }
else
{ int res = extended_gcd(b , a % b , y , x) ; y -= x * (a / b) ; return res ; }
}
这里有一个问题，我们在分析的时候，余数为0的最后一层，并不是我们的程序中的递归的最后一层。按照我们的分析，最后一层x = 0 , y = 1，而我们的程序判断的时候，是b = 0为最后的出口，所以x , y要交换一下，变成x = 1 , y = 0。
这个答案不是唯一的，当x+=b , y-=a，同样满足。
Part III 线性同余方程ax ≡ b(mod n)
我们的新问题又来了：解一个单独的线性同余方程ax ≡ b(mod n)
显然，它等价于求出这样的x，使得存在整数y，使
ax - by = b ... (1)
记d = (a , n) , a = a' * d , n = n' * d，那么必有d | b，否则方程无解。
这样方程就变成了
a'x - n'y = b/d ... (2)，(a' , n') = 1
我们在Part II中得到的x0 , y0是方程a'x0 - n'y0 = 1 ... (3)的解，
观察(2)和(3)，于是得到
x = x0 * (b / d) , y = y0 * (b / d)
注意的是，还要对n取一下模。
这个x同样不是唯一的，所有解属于同一个模n'剩余系（注意是n'）。
void modular_linear_equation(int a , int b , int n)
{
int e , x , y , d = extended_gcd(a , n , x , y) , ans[] ;
if (b % d)
report_no_solution ;
else
{
e = x * (b / d) % n ;
for (int i = 0 ; i < n ; i++)
ans[i] = (e + i * n / d) % n ;
}
}
Part IV 线性同余方程组
我们这里只考虑两个方程的情况。（再多我就不会了）
x ≡ a1 (mod m1) ... (1)
x ≡ a2 (mod m2) ... (2)
有解的充要条件是(m1 , m2) | a1 - a2。如果其成立，则方程组仅有一个小于m的非负整数解。
(1) 等价于
x = a1 + m1y
(2) 等价于
x = a2 + m2z
于是：
a1 + m1y = a2 + m2z
设c = (m1 , m2)，则：
(a1 - a2) / c = m2 / c * z - m1 / c * y
只要求得p , q使得：
(m2 / c) * p + (m1 / c) * q = 1
只需令：
z = p * (a1 - a2) / c
就得到一个解：
x = a2 + m2z
Part V 中国剩余定理
设m1 , m2 , ... , mk两两互素，则下面的方程组：
x ≡ a1 (mod m1)
x ≡ a2 (mod m2)
...
x ≡ ak (mod mk)
在0 <= x < M = m1m2...mk内有唯一解。
记Mi = M / mi，则(Mi , mi) = 1，故有pi , qi，使得Mipi + miqi = 1，记ei = Mipi，那么必有
ei ≡ 0(mod mj) (j != i)
ei ≡ 1(mod mi)
这样就可以看出，e1a1 + ... + ekak是一个解，不过不一定是最小的。
加减M的整数倍可以得到其他解。
int ChineseReminderTheory(int a[] , int m[] , int k)
{
int M , Mi , i , answer , pi , qi ;
for (M = 1 , i = 0 ; i < k ; i++) M *= m[i] ;
for (answer = i = 0 ; i < k ; i++)
{
Mi = M / m[i] ;
extended_gcd(Mi , m[i] , pi , qi) ;
answer = (answer + Mi * pi * a[i]) % M ;
}
if (answer < 0) answer += M ;
return answer ;
}
Postscript
浙公网安备 33010602011771号