Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'd like to implement the "Franklin-Reiter Related Message Attack" (see section 4.3 of Boneh's paper). As part of the implementation, I require to compute the GCD of two polynomials over $\mathbb{Z}_N[x]$, where $N$ is a composite integer (whose factorization is unknown). However, Mathematica's built-in PolynomialGCD allows the GCD to be taken over $\mathbb{Z}_p[x]$, where $p$ is a prime.

How can I compute polynomial GCDs over a ring (with composite characteristic)?

I think this question has similarities with one of my earlier questions (regarding the computation of the inverse of a polynomial in a polynomial ring). Unfortunately, I don't have enough math background to solve this question using the answers in the previous one.

Edit: Example, as per @DanielLichtblau's request:

Let $N = 91$ ($=7\times 13$), and consider the following polynomials:

f[x_] := PolynomialMod[(x + 1) (30 x + 40), 91]
g[x_] := PolynomialMod[(x + 1) (50 x + 60), 91]

Since $(x+1)$ divides both polynomials, I expect it to divide their GCD as well. Simply executing:

PolynomialGCD[f[x], g[x]]

will return 1, as it does not compute the GCD over the ring $\mathbb{Z}_{91}[x]$. Furthermore, running:

PolynomialGCD[f[x], g[x], Modulus -> 91]

returns error PolynomialGCD::modp: Value of option Modulus -> 91 should be a prime number or zero, as the ring's characteristic ($ = 91$) is a composite number, not a prime.

So, the question is, how can I actually compute the GCD of f[x] and g[x] over $\mathbb{Z}_{91}[x]$ using Mathematica?


Edit2: As pointed out above, the factorization of $N$ is unknown. Therefore, the algorithm cannot use the factors of $N$. (However, we know that $N$ is an RSA modulus, meaning that it is of the form $p\times q$, where $p$ and $q$ are primes.)

A general strategy is to use the extended Euclidean algorithm. Example:

r0 = g[x]
(* = 60 + 19 x + 50 x^2 *)
r1 = f[x]
(* = 40 + 70 x + 30 x^2 *)

We now find the inverse of the leading coefficient of r1 over $\mathbb{Z}_N$:

PowerMod[30, -1, 91]

which is 88. Using extended Euclidean algorithm:

r2 = PolynomialMod[r0 - 50*88*r1, 91]
(* = 54 + 54 x *)

Next, the inverse of the leading coefficient of r2 over $\mathbb{Z}_N$ is obtained:

PowerMod[54, -1, 91]

which is 59. Using extended Euclidean algorithm:

r3 = PolynomialMod[r1 - 30*59x*r2, 91]
(* = 40 + 40 x *)

The inverse of the leading coefficient of r3 over $\mathbb{Z}_N$ is 66:

PowerMod[40, -1, 91]

Using extended Euclidean algorithm:

r4 = PolynomialMod[r2 - 54*66*r3, 91]
(* = 0 *)

Therefore, the GCD of f[x] and g[x] equals r3 = 40 + 40 x = 40 (1+x).

Notice that if h[x] is the GCD of f[x] and g[x] over $\mathbb{Z}_N$, and if $a \not\mid N$, then a*h[x] is also a GCD of f[x] and g[x] over $\mathbb{Z}_N$. (That is, the GCD is not uniquely defined.) Consequently, 1+x is a valid GCD as well.


Note: If in any step above, the leading coefficient was not coprime to $N$, its inverse would not exist. In such cases, it seems reasonable to say that the GCD does not exist over the ring.

share|improve this question
1  
Would be useful to have an actual example, in Mathematica syntax, of input and desired output. –  Daniel Lichtblau Mar 10 at 16:52
add comment

1 Answer 1

For moduli that are square-free one can use Chinese remaindering on the coefficient lists to get a result valid for the moduli product.

cfs[p1_, p2_, x_, p_] := 
 Reverse[CoefficientList[PolynomialGCD[p1, p2, Modulus -> p], x, 
   1 + Min[Exponent[p1, x], Exponent[p2, x]]]]

FromDigits[
 ChineseRemainder[
  Transpose[{cfs[f[x], g[x], x, 7], cfs[f[x], g[x], x, 13]}], {7, 
   13}], x]

(* Out[19]= 1 + x *)

For a given pair of polynomials there may be finitely many primes for which this fails. I believe this corresponds to there not being a well defined GCD in such cases. Have not thought that through too carefully though.

--- edit ---

One can use GroebnerBasis in many cases. I show an example below.

Some background, albeit tailored to the case of Hensel lifting, can be found here.

f[x_] := 
 PolynomialMod[(x + 135)^2*(30 x + 40)*(x^2 + 44 x + 8), 91]
g[x_] := PolynomialMod[(x + 135)*(50 x + 60)*(x^2 + 2 x + 41), 91]

If we get a basis with the modulus and only one other generator then the method succeeded.

ideal = 
 GroebnerBasis[{f[x], g[x], 7 13}, x, CoefficientDomain -> Integers]

(* Out[85]= {91, 44 + x} *)

That is promising. The gcd is then the second element.

gcd = Last[ideal]

(* Out[81]= 44 + x *)

One can check that this is the same result that the previous method would give. Also we explicitly test that the this divide both inputs.

Map[PolynomialMod[#, {gcd, 91}] &, {f[x], g[x]}]

(* Out[87]= {0, 0} *)

--- end edit ---

--- edit 2 ---

One can get faster code using an internal function that operates on univariate polynomials represented as lists of coefficients.

FromDigits[
 Reverse[Algebra`PolynomialGCDModList[CoefficientList[f[x], x], 
   CoefficientList[g[x], x], 7*13]], x]

(* Out[96]= 44 + x *)

If the extended gcd is needed there is also Algebra`PolynomialExtendedGCDModList.

--- end edit 2 ---

share|improve this answer
    
Thanks a lot for the answer. I think it uses the factorization of $N$, which as pointed out in the question, is unknown. Please see my edited question for a suggestion on how the GCD can be computed (it uses the extended Euclidean algorithm). Could you please help me implement it using Mathematica? –  Sadeq Dousti Mar 11 at 3:44
    
See edit for method that does not factor the modulus. It is pretty much identical to Franklin-Reiter's use of a Euclidean remainder sequence. I implemented it using GroebnerBasis simply to save on coding. –  Daniel Lichtblau Mar 11 at 15:22
    
Thanks again. Your answers are always inspiring! GroebnerBasis works like a charm, but I guess it's slower than extended Euclidean algorithm. In the real-world problem I'm facing, f[x_]:=x^e and g[x_]:=(x+c)^e for some constant c. Here, $e \approx 16,000$ and $N$ is a 1024-bit integer. When I executed GroebnerBasis on such large values, it took Mathematica a long time, and after a while I aborted the execution. I guess extended Euclidean algorithm would execute much faster. Do you agree? If so, could you please provide me with the faster code? –  Sadeq Dousti Mar 11 at 18:17
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.