Search code examples
pythonpython-3.xsympypolynomials

Changing coefficients modulo p in a SymPy polynomial


I took a cryptography course this semester in graduate school, and once of the topics we covered was NTRU. I am trying to code this in pure Python, purely as a hobby. When I attempt to find a polynomial's inverse modulo p (in this example p = 3), SymPy always returns negative coefficients, when I want strictly positive coefficients. Here is the code I have. I'll explain what I mean.

import sympy as sym
from sympy import GF

def make_poly(N,coeffs):
    """Create a polynomial in x."""
    x = sym.Symbol('x')
    coeffs = list(reversed(coeffs))
    y = 0
    for i in range(N):
        y += (x**i)*coeffs[i]
    y = sym.poly(y)
    return y

N = 7
p = 3
q = 41
f = [1,0,-1,1,1,0,-1]

f_poly = make_poly(N,f)

x = sym.Symbol('x')

Fp = sym.polys.polytools.invert(f_poly,x**N-1,domain=GF(p))
Fq = sym.polys.polytools.invert(f_poly,x**N-1,domain=GF(q))

print('\nf =',f_poly)
print('\nFp =',Fp)
print('\nFq =',Fq)

In this code, f_poly is a polynomial with degree at most 6 (its degree is at most N-1), whose coefficients come from the list f (the first entry in f is the coefficient on the highest power of x, continuing in descending order).

Now, I want to find the inverse polynomial of f_poly in the convolution polynomial ring Rp = (Z/pZ)[x]/(x^N - 1)(Z/pZ)[x] (similarly for q). The output of the print statements at the bottom are:

f = Poly(x**6 - x**4 + x**3 + x**2 - 1, x, domain='ZZ')
Fp = Poly(x**6 - x**5 + x**3 + x**2 + x + 1, x, modulus=3)
Fq = Poly(8*x**6 - 15*x**5 - 10*x**4 - 20*x**3 - x**2 + 2*x - 4, x, modulus=41)

These polynomials are correct in modulus, but I would like to have positive coefficients everywhere, as later on in the algorithm there is some centerlifting involved, so I need to have positive coefficients. The results should be

Fp = x^6 + 2x^5 + x^3 + x^2 + x + 1
Fq = 8x^6 + 26x^5 + 31x^4 + 21x^3 + 40x^2 + 2x + 37

The answers I'm getting are correct in modulus, but I think that SymPy's invert is changing some of the coefficients to negative variants, instead of staying inside the mod.

Is there any way I can update the coefficients of this polynomial to have only positive coefficients in modulus, or is this just an artifact of SymPy's function? I want to keep the SymPy Poly format so I can use some of its embedded functions later on down the line. Any insight would be much appreciated!


Solution

  • This seems to be down to how the finite field object implemented in GF "wraps" integers around the given modulus. The default behavior is symmetric, which means that any integer x for which x % modulo <= modulo//2 maps to x % modulo, and otherwise maps to (x % modulo) - modulo. So GF(10)(5) == 5, whereas GF(10)(6) == -4. You can make GF always map to positive numbers instead by passing the symmetric=False argument:

    import sympy as sym
    from sympy import GF
    
    def make_poly(N, coeffs):
        """Create a polynomial in x."""
        x = sym.Symbol('x')
        coeffs = list(reversed(coeffs))
        y = 0
        for i in range(N):
            y += (x**i)*coeffs[i]
        y = sym.poly(y)
        return y
    
    N = 7
    p = 3
    q = 41
    f = [1,0,-1,1,1,0,-1]
    
    f_poly = make_poly(N,f)
    
    x = sym.Symbol('x')
    
    Fp = sym.polys.polytools.invert(f_poly,x**N-1,domain=GF(p, symmetric=False))
    Fq = sym.polys.polytools.invert(f_poly,x**N-1,domain=GF(q, symmetric=False))
    
    print('\nf =',f_poly)
    print('\nFp =',Fp)
    print('\nFq =',Fq)
    

    Now you'll get the polynomials you wanted. The output from the print(...) statements at the end of the example should look like:

    f = Poly(x**6 - x**4 + x**3 + x**2 - 1, x, domain='ZZ')
    Fp = Poly(x**6 + 2*x**5 + x**3 + x**2 + x + 1, x, modulus=3)
    Fq = Poly(8*x**6 + 26*x**5 + 31*x**4 + 21*x**3 + 40*x**2 + 2*x + 37, x, modulus=41)
    

    Mostly as a note for my own reference, here's how you would get Fp using Mathematica:

    Fp = PolynomialMod[Algebra`PolynomialPowerMod`PolynomialPowerMod[x^6 - x^4 + x^3 + x^2 - 1, -1, x, x^7 - 1], 3]
    

    output:

    1 + x + x^2 + x^3 + 2 x^5 + x^6