一、问题阐述：

         对形如：ax^2+bx+c = 0 (mod p)； (A)

         有 x, a, b, c, p皆为整数，求出所有可能的解x

        

         问题可抽象为：

         输入： int a, b, c ,p

         输出： int x， 使其满足方程(A)

 

二、算法思想：

         为简化问题，可先认为p为素数，且0<=a,b,c,x<p

         当a=0时，(A)问题可转化为如 ax+b = 0 (mod p)(B)

         当a!=0时，可使用配方法将(A)问题转化为求U^2+D = 0 (mod p)(C)的解，具体为:

1)       当p=2，由于x<p，故x只能取0或1；若c=0，则显然x=0是其中一个解，否则x=0不是问题的解；当a+b+c mod p为0则x=1为问题的解，否则x=1不是问题的解。

2)       当p为奇素数时，若a^(p-1)/2≡1 (mod p)不成立则原问题无解。否则一律将(A)转变为(C)，其中U=2Ax+B 而 D=(4AC-B²)。解问题(C)可采用Shank-Tonelli算法，当解出一个U时，在[0,p-1]范围内p-U为问题(C)的另外一个解，因此原问题也对应着两个解，当然存在这两个解相同的情况。剩下的步骤将U代入(B)求解x

3)       求解(B)使用Extended Euclid算法

 

Shank-Tonelli算法[1]如下：

1)       求k, 使得p-1 = 2ⁿk

2)       求出能满足q^(p-1)/2≡ -1 (mod p)最小的q

3)       t = a^(k+1)/2 % p and let r = a^k%p

4)       t为问题的一个解，如果能满足r≡1 (mod p) ，否则转5)

5)       求出最小的i使得满足r²ⁱ≡1(mod p)

6)       分别计算e = 2^(n-i-1) u = q^(ke)%P， t = (t*u)%P ， r = (r*u*u)%P

7)       转到4

 

当p为合数时，不能直接使用上面的解法。由孙子定理[3]可知(A)的解为将p分解成如p₁^k1p₂^k2…p_n^kn时对应各子同余式在扩展解空间情况下的公共解，其中p₁，p₂，…p_n为素数。

如：x^2 mod 6 = 0 解为2和4。而X^2 mod 2 = 0 解为0，X^2 mod 3 = 0 解为1和2，扩展解空间时，x+kp均为原问题的解，因此在[0,6]范围内，X^2 mod 2 = 0解为0和2，4，X^2 mod 3 = 0解为1、2、4、5，因此其公共解为2和4

故此时(A)问题转化为求ax^2+bx+c = 0 (mod p^k)(D)的解，求此类问题可使用Hensel's_lemma[2]，在此就不列出其算法。

因此当p为合数时算法步骤为

1)       将p分解为p₁^k1p₂^k2…p_n^kn

2)       针对每个p_i，使用Hensel’s lemma求出(D)的解

3)       在扩展解空间后求出所有pi的公共解

 

三、算法代码(p为素数部分)：

 

         Input: a, b, c, p

         Output: {X}

         procedure solve_all_problem (a, b, c, p)

                   if a== 0

                            if b==0&&c!=0

                                     return 0; // 在[0,p-1]内 0 是唯一的解

                            else if(b!=0){

                                     int result = solve_linear(b,-c,p); //利用Extended Euclid算法解线性方程

                                     return result;

                            else

                                     return NO_SOLUTION;

                            end if

         else if p == 2

                   int solution_a = -1,solution_b = -1; //初始化为无解

                   if (c%2) == 0

                     solution_a = 0;

                   end if

                   if (a + b + c)%p == 0

                     solution_b = 1;

                   return solution_a and solution_b //两个可能的解

                   end if       

         else

                   int u;

                   int d = ((long long)b*b-4LL*a*c)%p;

                   u = shank_tonelli(d,p);

                   if u==-1

                            return NO_SOLUTION

                   else

                            int solution_a = solve_linear(2*a,u-b,p);

                            int solution_b = solve_linear(2*a,p-u-b,p);

                            if(solution_a==solution_b){

                                     return solution_a //只有一个解

                            else

                                     return solution_a and solution_b

                            end if                 

                   end if

         end if

        

四、测试用例：

input:2 3 4 6

output: 2 2 4

-------------------

input:4 3 3 13

output: 1 11

-------------------

input:17 8 1 71

output: 0

-------------------

input:2 3 5 7

output: 2 5 4

-------------------

input:0 4 0 7

output: 1 0

-------------------

input:2 34 9 19

output: 2 6 15

-------------------

input:12 5 27 2

output: 1 1

-------------------

input:1 179344794 146367396 179424691

output: 2 1876 78021

 

 

五、参考文献

[1] Tutorial for Quadratic Equations:

 http://www.codechef.com/wiki/tutorial-quadratic-equations

[2] Hensel’s Lemma: http://en.wikipedia.org/wiki/Hensel's_lemma

[3] 孙子定理: http://baike.baidu.com/view/157384.htm