# Polynomial GCD over a ring (with composite characteristic)

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.

• Would be useful to have an actual example, in Mathematica syntax, of input and desired output. – Daniel Lichtblau Mar 10 '14 at 16:52

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[AlgebraPolynomialGCDModList[CoefficientList[f[x], x],
CoefficientList[g[x], x], 7*13]], x]

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


If the extended gcd is needed there is also AlgebraPolynomialExtendedGCDModList.

--- end edit 2 ---

• 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? – M.S. Dousti Mar 11 '14 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 '14 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? – M.S. Dousti Mar 11 '14 at 18:17