GCD of Polynomials over a finite field for NFS
 

I'm stuck at the square root step of the number field sieve. Specifically the part that verifies that primes are inert. I'm supposed to calculate a GCD of two polynomials modulo a prime to verify that such prime is usable for Chinese remaindering. I can't seem to get the code right.

These numbers were taken directly from Briggs paper at page 58 and the code was based on [url]http://stackoverflow.com/a/26178457/3933560[/url] . I've actually implemented my own, but I chose to post this one because it's cleaner. Either way, I get the same wrong results. What am I missing?

EDIT: The code is python.

[CODE]from itertools import izip

def normalize(poly):
while poly and poly[-1] == 0:
poly.pop()
if poly == []:
poly.append(0)


def poly_remainder_GFp(num, den, p):
#Create normalized copies of the args
num = num[:]
normalize(num)
den = den[:]
normalize(den)

if len(num) >= len(den):
#Shift den towards right so it's the same degree as num
shiftlen = len(num) - len(den)
den = [0] * shiftlen + den
else:
return num

quot = []
divisor = den[-1]
for i in xrange(shiftlen + 1):
#Get the next coefficient of the quotient.
mult = (num[-1] * modinv(divisor,p)) % p
#quot = ([mult] + quot) % p

#Subtract mult * den from num, but don't bother if mult == 0
#Note that when i==0, mult!=0; so quot is automatically normalized.
if mult != 0:
d = [mult * u for u in den]
num = [(u - v)%p for u, v in zip(num, d)]

num.pop()
den.pop(0)

normalize(num)
return num

def egcd(a, b):
if a == 0:
return (b, 0, 1)
else:
g, y, x = egcd(b % a, a)
return (g, x - (b // a) * y, y)

def modinv(a, m):
g, x, y = egcd(a, m)
if g != 1:
raise Exception('modular inverse does not exist')
else:
return x % m

if __name__ == '__main__':
N = [7301, 1477, 7726]
D = [8, 29, 15, 1]
while (D != [0]):
(N,D) = (D,poly_remainder_GFp(N, D,9923))
print "%s, %s" % (N,D)

print ""

N = [5984,4697,7449]
D = [8,29,15,1]
while (D != [0]):
(N,D) = (D,poly_remainder_GFp(N, D,9929))
print "%s, %s" % (N,D)[/CODE]

PS: I can't seem to edit the post anymore. Anyway, the difference of the code from the stackexchange link is that I reduce the coefficients modulo p, and I used an inverse modulo operation instead of division.

The code that runs first is the one at the bottom. I used a while loop to calculate the GCD. It terminates when the remainder becomes zero. Running this gives me [9858, 7744] and [1548] (which corresponds to 9858 + 7744x and 1548) which are both wrong.

And one more thing, I've implemented the actual Chinese remainder part, it's just the prime checking that's missing now.

Running in Pari, I get
[code]
? f=8 + 29*x + 15*x^2 + x^3
%1 = x^3 + 15*x^2 + 29*x + 8
? f*(Mod(1,9923))
%2 = Mod(1, 9923)*x^3 + Mod(15, 9923)*x^2 + Mod(29, 9923)*x + Mod(8, 9923)
? g=7726*x^2+1477*x+7301
%3 = 7726*x^2 + 1477*x + 7301
? g*Mod(1,9923)
%4 = Mod(7726, 9923)*x^2 + Mod(1477, 9923)*x + Mod(7301, 9923)
? gcd(f*Mod(1,9923),g*Mod(1,9923))
%5 = Mod(7744, 9923)*x + Mod(9858, 9923)
[/code]

which looks the same as the result you get and claim is wrong. Obviously any multiple of that linear polynomial would work ...

[code]
? f2 = f-(4214*x)*g
%10 = Mod(7581, 9923)*x^2 + Mod(4838, 9923)*x + Mod(8, 9923)
? Mod(7726/7581,9923)
%18 = Mod(6081, 9923)
? g2 = g-f2*6081
%23 = Mod(3294, 9923)*x + Mod(8268, 9923)
? f3 = f2-(6464*x)*g2
%46 = Mod(5764, 9923)*x + Mod(8, 9923)
? 5995*f3-g2
%55 = Mod(0, 9923)
[/code]

I'm confused now. The two polynomials (7726*x^2+1477*x+7301 and x^3 + 15*x^2 + 29*x + 8) both have the same root 847 mod 9923 since:
(7726*847*847 + 1477*847 + 7301) mod 9923 = 0
(847*847*847 + 15*847*847 + 29*847 + 8) mod 9923 = 0

So, x-847 should be the GCD, again from Brigg's paper (page 58, [url]http://scholar.lib.vt.edu/theses/available/etd-32298-93111/unrestricted/etd.pdf[/url] ).

[QUOTE=paul0;392530]I'm confused now. The two polynomials (7726*x^2+1477*x+7301 and x^3 + 15*x^2 + 29*x + 8) both have the same root 847 mod 9923 since:
(7726*847*847 + 1477*847 + 7301) mod 9923 = 0
(847*847*847 + 15*847*847 + 29*847 + 8) mod 9923 = 0

So, x-847 should be the GCD, again from Brigg's paper (page 58, [url]http://scholar.lib.vt.edu/theses/available/etd-32298-93111/unrestricted/etd.pdf[/url] ).[/QUOTE]

EDIT: It turns out, 9858 + 7744x = x - 847 mod 9923 since 9858/7744 = -847 mod 9923. I'm still trying to figure out the other polynomial (7449x^2 + 4697x + 5984) whose GCD should result to 1 instead of 1548.

EDIT2: It looks like my GCD code return a polynomial degree 0 (a coefficient) if the GCD is 1. The proper GCD is returned otherwise.

The GCD is only meaningfully defined up to multiplication by a unit (IE in this case a non-zero field element); so if it's a degree-0 polynomial it might as well be 1. For polynomial GCD, the "greatest" in the names means 'largest degree'.

The degree of the GCD is well-defined, the set of roots of the GCD is well-defined, but you can pick any leading coefficient you want.

[QUOTE=fivemack;392566]The GCD is only meaningfully defined up to multiplication by a unit (IE in this case a non-zero field element); so if it's a degree-0 polynomial it might as well be 1. For polynomial GCD, the "greatest" in the names means 'largest degree'.

The degree of the GCD is well-defined, the set of roots of the GCD is well-defined, but you can pick any leading coefficient you want.[/QUOTE]

Thanks! :)