GCD of Polynomials over a finite field for NFS

[paul0]

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 http://stackoverflow.com/a/26178457/3933560 . 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.

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)

[paul0]

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.

[fivemack]

Running in Pari, I get

? 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)

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

? 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)

[paul0]

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, http://scholar.lib.vt.edu/theses/ava...ricted/etd.pdf ).

[paul0]

Originally Posted by paul0

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, http://scholar.lib.vt.edu/theses/ava...ricted/etd.pdf ).

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.

[fivemack]

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.

[paul0]

Originally Posted by fivemack

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.

Thanks! :)