ACM：数论专题(6)——模线性方程组

题目描述：

    给定n组除数Mi和余数Ri (1≤i≤n, Ri<Mi), 定义方程组如下：
        x % M1 = R1
        x % M2 = R2
        x % M3 = R3
            ...
        x % Mn = Rn
求解满足以上方程组的x的最小正整数解。

解答：
    题目中定义的方程组称为“模线性方程组”，我们可以采用一种迭代的方式进行求解。

·基本情况：
    首先考虑方程组中只包含2个方程的情况：
        x % M1 = R1..................................................(F1)
        x % M2 = R2..................................................(F2)
假设： x = k1 * M1 + R1, x = k2 * M2 + R2，此时：
        k1 * M1 + R1 = k2 * M2 + R2
    =>  k1 * M1 - k2 * M2 = R2 - R1..................................(T)
M1, M2, R1, R2 都是已知的常数，这样就得到了一个以k1, k2为未知数的二元一次方程。此时，只要找到(T)式的一组整数解，然后代入，就可以找到满足由方程(F1) 和 (F2)的一个x解。
    
    求解(iii)式的过程可以利用拓展欧几里得算法。 简要说明其步骤如下：
    记：A = M1, B = M2, C = R2 - R1, x = k1, y = k2那么方程(iii)就可以抽象为：
             Ax + By = C..........................................(E1)
根据拓展欧几里得定理，方程E1有整数解(x 和 y均为整数)的充要条件是：C可以被A和B的最大公约数整除，即：C % gcd(A, B)。利用这个条件就可以判定方程是否有整数解。
    如果通过以上判定方式判定方程有整数解，则可以利用下面的方法，找出方程的一组解：
    记：方程E2: 
            Bx + (A%B)y = gcd(B, A%B)..............................(E2)
由辗转相除法可知：gcd(A, B) = gcd(B, A%B), 因此：方程E1和E2的等号右边项是相同的，因此：
            Ax + By = Bx + (A%B)y = Bx + [A - (A/B)*B]y = Bx + Ay - B*(A/B)y
         => Ax + By = Ay + B * [x - (A/B)y]
此时，如果已经找到方程E2的一组解(x2, y2), 代入上式可得：
            Ax + By = Ay2 + B * [x2 - (A/B)y2]
因此，x = y2, y = x2 - (A/B)y2 就是方程E1的一组解。
    根据以上特性，给出方程E1后，我们可以将其转化为方程E2， 然后求解E2后，得到E1的解。至于方程E2的求解，则可以递归地进行，令：E2中：A = B, B = A%B, 将其转化为E3,再通过求解E3, 得到E2，如此循环往复，直到某一轮迭代的方程Ei中，参数A以被参数B整除，即：A%B = 0, 那么此时，gcd(A, B) = B, 方程就化为：Ax + By = B. 显然，这个方程有一组解为：x = 0, y = 1, 然后再依次回溯至方程：Ei-1, Ei-2 ... E1, 就可以求解出原方程的一组解。
    求解完成后，将解得的k1的值代入：k1*M1 + R1, 就能得到一个满足方程(F1)和(F2)的一个x值。
    以上方法只能保证可以找到满足方程组的某一个解，但却不能保证题目中要求的最小正整数解。为了得到最小正整数解，还需要以下的处理：假设已经找到了方程E1的一组解：x1, y1, 并记此时计算得到的满足方程(F1)和(F2)的x值为x0。将x1，y1代入方程E1得：
        Ax1 + By1 = gcd(A, B) 
    =  Ax1 + By1 + sAB - sAB 
    =  A(x1-sB) + B(y1+sA) = A(x1+sB) + B(y1-sA)
其中s为任意整数。
    以上式子表明：如果x1, y1是原方程的解，那么x1±sB， y1±sA 也是原方程的解。对应于：A = M1, B = M2, C = R2 - R1, x = k1, y = k2， 我们可以得到一组x值：
        x = k1 * M1 + R1 = (x1±B)*M1 + R1 = x1*M1 + R1 ± BM1 = x0±s*M2*M1 
因此通过调整上式中的s的值，就一定可以找到满足原方程组的最小正整数解。

·拓展：
    上文中提供了一种求解仅包含2个方程的模线性方程组的解的方法。而对于包含多个方程的方程组，则可以采用迭代的策略进行求解，方法如下：
    假设原方程组中包含n个方程，分别记作：F1, F2, ... Fn
    联立方程F1和F2， 根据上文中的方法，F1和F2组成的方程组求得的解：
            x = x0 + s*M2*M1
此时，上式可以转换为：
            x % (M2*M1) = x0 .................................(G2)
    记上面的方程为G2， 显然，满足方程G2的x值一定满足方程F1和方程F2。此时联立G2和F3，求解方程组后，得到的解就一定同时满足方程：F1, F2, F3。此时，利用类似的方法，可以将G2和F3组成的方程组的解转换为方程G3，然后联立G3和F4求解方程组......以此类推。
    记方程F1为G1，我们可以发现：
        求解方程组G1和F2，得到的解满足F1和F2，同时得到方程G2；
        求解方程组G2和F3，得到的解满足F1, F2, F3, 同时得到方程G3;
        求解方程组G3和F4，得到的解满足F1, F2, F3, F4, 同时得到方程G4;
            .....
        求解方程组Gn-2和Fn-1，得到的解满足F1, F2, F3, ... , Fn-1, 同时得到方程Gn-1
        求解方程组Gn-1和Fn，得到的解满足F1, F2, F3, ... , Fn，此时不需要得到方程Gn, 因为得到的解满足原方程组中的所有方程，求解完毕。  
         
输入输出格式：
    输入：第1行为1个正整数N,表示方程组中的方程的个数；接下来的N行，每行2个数字，表示Mi和Ri
    输出：输出一个整数，表示方程组的最小正整数解，如果方程组无解，则输出-1。  

数据范围：
     2 ≤ Mi ≤ 20,000,000
     0 ≤ Ri < Mi 

程序代码：

/****************************************************//* File        : Hiho_Week_97                       *//* Author      : Zhang Yufei                        *//* Date        : 2016-05-09                         *//* Description : HihoCoder ACM program. (submit:g++)*//****************************************************/#include<stdio.h>#include<stdlib.h>/* * Define the structure to record the result. * Parameters: *@x & @y: The result (x, y). */typedef struct node {long long x;long long y;} result;/* * This function computes the greatest common divisor of 2 numbers a and b. * Parameters: *@a & @b: The 2 numbers to compute. * Returns: *The greatest common divisor of a and b. */long long gcd (long long a, long long b) {if(a % b == 0) {return b;} else {return gcd (b, a % b);}}/* * This function computes the lowest common multiple of 2 numbers a and b. * Parameters: *@a & @b: The 2 numbers to compute. * Returns: *The lowest common multiple of a and b. */long long lcm (long long a, long long b) {long long g = gcd (a, b);return (a / g) * b;}/* * This function calculates the result of the equation: *Ax + By = 1 using extend euclidean theory. * Parameters: * @A & @B: The parameter A, B in the equation. * Returns: * One result of this equation, if the equation don't have  *a legal result, returns NULL. */ result* extend_euclidean (long long A, long long B) {result *r; if(A % B == 0) {r = (result*) malloc (sizeof(result));r->x = 0;r->y = 1;} else {r = extend_euclidean (B, A % B);long long x = r->x;long long y = r->y;r->x = y;r->y = x - A / B * y;}return r;}/* * The main program. */int main (void) {int N;scanf ("%d", &N);long long m1, r1, m2, r2;long long A, B, C;scanf ("%lld %lld", &m1, &r1);result *r;long long answer = 0;for(int i = 1; i < N; i++) {scanf ("%lld %lld", &m2, &r2);if(answer == -1) {continue;}A = m1;B = m2;C = r2 - r1;long long g = gcd (A, B);if(C % g != 0) {answer = -1;continue;}long long lc = lcm (A, B);A /= g;B /= g;C /= g;r = extend_euclidean (A, B);r->x *= C;r->x %= B;answer = m1 * r->x + r1;while(answer < 0) {answer += lc;}while(answer > lc) {answer -= lc;}m1 = lc;r1 = answer;} printf("%lld\n", answer);return 0;}