mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Programming

Reply
 
Thread Tools
Old 2015-01-15, 09:56   #1
paul0
 
Sep 2011

5710 Posts
Default 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 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.

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)

Last fiddled with by paul0 on 2015-01-15 at 10:07
paul0 is offline   Reply With Quote
Old 2015-01-15, 11:27   #2
paul0
 
Sep 2011

718 Posts
Default

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.

Last fiddled with by paul0 on 2015-01-15 at 11:44
paul0 is offline   Reply With Quote
Old 2015-01-15, 13:03   #3
fivemack
(loop (#_fork))
 
fivemack's Avatar
 
Feb 2006
Cambridge, England

22×32×179 Posts
Default

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)
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)
fivemack is offline   Reply With Quote
Old 2015-01-15, 23:24   #4
paul0
 
Sep 2011

1110012 Posts
Default

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

Last fiddled with by paul0 on 2015-01-15 at 23:29
paul0 is offline   Reply With Quote
Old 2015-01-16, 01:07   #5
paul0
 
Sep 2011

1110012 Posts
Default

Quote:
Originally Posted by paul0 View Post
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.

Last fiddled with by paul0 on 2015-01-16 at 02:00
paul0 is offline   Reply With Quote
Old 2015-01-16, 09:52   #6
fivemack
(loop (#_fork))
 
fivemack's Avatar
 
Feb 2006
Cambridge, England

22·32·179 Posts
Default

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.

Last fiddled with by fivemack on 2015-01-16 at 09:53
fivemack is offline   Reply With Quote
Old 2015-01-16, 15:12   #7
paul0
 
Sep 2011

3·19 Posts
Default

Quote:
Originally Posted by fivemack View Post
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! :)
paul0 is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Finite differences for the win George M Miscellaneous Math 2 2018-01-01 21:37
Finite field extensions carpetpool Abstract Algebra & Algebraic Number Theory 3 2017-11-15 14:37
Polynomials defining the same field as cyclotomic polynomial order 5 carpetpool Abstract Algebra & Algebraic Number Theory 0 2017-04-19 20:33
Polynomials with constant terms replaced defining the same field carpetpool Miscellaneous Math 2 2017-02-25 00:50
Is there a finite-field equivalent to the DWT? fivemack Math 4 2008-03-27 17:58

All times are UTC. The time now is 12:41.


Mon Dec 6 12:41:41 UTC 2021 up 136 days, 7:10, 0 users, load averages: 1.41, 1.91, 1.72

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.