Jsevillamol's profile - activity

Consider we have a polynomial $g$ and an integer $p$ defined as follows:

ZZX.<x> = ZZ[x]
p = 953
g = x+951

We want to find the unique polynomial which is congruent with g in $ZZ/ (p) [x]$ and has minimal max-norm, that is h = x-2.

How do I do that programatically in sagemath?

EDIT:

My first attempt is this code:

min_repr = lambda n, m : ZZ(mod(n,m)) if abs(ZZ(mod(n,m))) < m / 2 else ZZ(mod(n,m)) - m
min_norm = lambda g, m : ZZX([min_repr(coeff, m) for coeff in g])

Trying it out it seems to fail on some cases though.

EDIT: I found out I can use sage.rings.polynomial.polynomial_ring.is_PolynomialRing to selectively apply this solution

Any ideas on how to implement this without breaking the code for other types of euclidean domains which are not polynomial rings (eg integers and gaussians)?

Consider the following snippet, intended to compute the extended euclidean algorithm for polynomials in $F_q[x]$.

def extended_euclides(a,b,quo=lambda a,b:a//b):
    r0 = a; r1 = b
    s0 = 1; s1 = 0
    t0 = 0; t1 = 1

    while r1 != 0:
        q = quo(r0, r1)
        r0, r1 = r1, r0 - q * r1
        s0, s1 = s1, s0 - q * s1
        t0, t1 = t1, t0 - q * t1

    return r0, s0, t0

GFX.<x> = GF(5)[x]
g = x^3 + 2*x^2 - x - 2
h = x - 2
r,s,t = extended_euclides(f,g)
print("The GCD of {} and {} is {}".format(g,h,r))
print("Its Bézout coefficients are {} and {}".format(s,t))
assert r == gcd(g,h), "The gcd should be {}!".format(gcd(g,h))
assert r == s*g + t*h, "The cofactors are wrong!"

Executing this snippet produces the following output:

The GCD of x^3 + 2*x^2 + 4*x + 3 and x + 3 is 3*x + 1
Its Bézout coefficients are 2*x + 4 and 3*x^2 + 4

Error in lines 45-45
Traceback (most recent call last):
  File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute
    flags=compile_flags) in namespace, locals
  File "", line 1, in <module>
AssertionError: The gcd should be 1!

That is, my code disagrees with sagemath's gcd function, and I think sagemath's result is correct; g(2)!=0, so g and h should not have common factors.

Why is mine wrong?

I am trying to implement a factorization algorithm for polynomials with integer coefficients, based on Kronecker's principle.

# Kronecker Factorization Algorithm
def kronecker_splitting(f):
    """
    INPUT:
        a polynomial f in ZZ[x]
    OUTPUT:
        a non trivial factor of f, or f if it is irreducible
    """
    ZZX.<x> = ZZ[]
    if f.parent() != ZZ[x]: raise ValueError("f must be in ZZ[x]")
    if f.degree() == 1: return f
    a = []
    a.append(ZZ.random_element())

    # Check if a0 is a root of f
    if f.subs(a[0])==0: return x-a[0]
    div = f.subs(a[0]).divisors() # compute positive divisors
    M = [(d,) for d in div] # Turn m into a collection of tuples

    for e in range(1,ceil(f.degree()/2)+1):
        # Choose a not previously chosen integer
        while(True):
            aux = ZZ.random_element()
            if aux not in a:
                a.append(aux)
                break
        # Check that a[e] is not a root
        if f.subs(a[e])==0: return x-a[e]

        # compute divisors of a[e]
        pos_div = f.subs(a[e]).divisors()
        neg_div = [-d for d in pos_div]
        div = pos_div + neg_div

        Mp = M = [c_prev + (d,) for c_prev in M for d in div]
        while len(Mp)>0:
            c = Mp.pop()
            # solve linear system
            A = matrix.vandermonde(a).transpose()
            b = A.solve_right(vector(c))
            vector2poly = lambda h: sum([ci*x^i for i,ci in enumerate(h)])
            g = vector2poly(b)
            if b[e] != 0 and f % g == 0: 
                print("found non trivial factor {}".format(g))
                print(g.parent())
                return g.numerator()

    return f # f is irreducible

def kronecker_factorization(f):
    """
    INPUT:
        f :- monic element of ZZ[x]
    OUTPUT:
        [f_1,...,f_s] :- monic irreducible factorization of f
    """
    ZZX = ZZ[x]
    factors = [f]; irreducible_factors = []

    while len(factors) > 0:
        g = factors.pop()
        print("splitting {}".format(g))
        h = kronecker_splitting(g)
        # if g was already irreducible
        if h == g: irreducible_factors.append(g)
        else: factors.extend([h, ZZX(g/h)])

    return irreducible_factors


ZZX.<x> = ZZ[x]
f = (x+1)*(x+2)*(x+3)
factors = kronecker_factorization(f)

print("{} factorizes as {}".format(f, factors))
assert(f == reduce(lambda a,b:a*b, factors))
for g in factors:
    assert(g.parent() == ZZX)
    assert(g.is_irreducible())

When I run the snippet above, I often get this output (try running it a few times to reproduce it, bcs of the randomness):

splitting x^3 + 6*x^2 + 11*x + 6
splitting x^2 + 3*x + 2
found non trivial factor 2*x + 2
Univariate Polynomial Ring in x over Rational Field
Error in lines 46-46
Traceback (most recent call last):
  File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute
    flags=compile_flags) in namespace, locals
  File "", line 1, in <module>
  File "", line 15, in kronecker_factorization
  File "sage/structure/parent.pyx", line 921, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9679)
    return mor._call_(x)
  File "sage/categories/map.pyx", line 790, in sage.categories.map.Map._call_ (build/cythonized/sage/categories/map.c:7325)
    cpdef Element _call_(self, x):
  File "/ext/sage/sage-8.4_1804/local/lib/python2.7/site-packages/sage/rings/fraction_field.py", line 1155, in _call_
    raise TypeError("fraction must have unit denominator")
TypeError: fraction must have unit denominator

This confuses me, since if I try reproducing the mistake with the following snippet:

QQX.<x> = QQ[x]
g = 3*x + 3
g.parent()
g.numerator()

I get the output: I expected, no errors in sight:

Univariate Polynomial Ring in x over Rational Field
3*x + 3

What is going on?

Thanks to the kind and patient help of @rburing this finally seems to be working!

vector2poly = lambda h: reduce(lambda a,b:a+b, [ci*x^i for i,ci in enumerate(h)])

def berlekamp_basis(f):
        """
        INPUT: f a monic squarefree polynomial in F_q[x]
        OUTPUT: a list of polynomials hi in F_q[x] that represent
                a basis of the Berlekamp algebra associated to f
        """
        n = f.degree()
        GFX = f.parent()
        x = GFX.gen()
        q = GFX.base().order()
        GFX_quo_f = GFX.quotient(f)
        g = []
        for k in range(n):
            basic_image = GFX_quo_f(x^(q*k)- x^k)
            g.append(vector(basic_image))
        H = matrix(g).transpose().right_kernel().basis()
        H = [h.lift() for h in H]
        H = map(vector2poly, [h.lift() for h in H])
        return H

Let $f\in GF_{q}[x]$ be a squarefree univariate polynomial with coefficients in a finite field, and $\bar{Fr}: \frac{F_q}{(f)}\to \frac{F_q}{(f)}: \bar h \mapsto \bar h^q$ the Frobenius application over $GF_q$ elevated to the cocient $\frac{F_q}{(f)}$.

Then we can see $\frac{F_q}{(f)}$ as a vector space over $GF_q$, and the space of fixed points of $\bar{Fr}$, that is $W = {\bar h \in \frac{F_q}{(f)} : \bar h ^q = q }$ is a linear subspace.

The question: how do I compute the basis of $W$ in sagemath?

EDIT:

Inspired by @rburing I wrote some (incomplete) code:

def berlekamp_basis(f):
    """
    INPUT: f a monic squarefree polynomial in F_q[x]
    OUTPUT: a list of polynomials hi in F_q[x] that represent
            a basis of the Berlekamp algebra associated to f
    """
    n = f.degree()
    GFX = f.parent()
    GFX_quo_f = GFX.quotient(f)
    g = []
    for k in range(n):
        basic_image = GFX(GFX_quo_f(x)^(q*k)) - GFX(x)^k
        g.append(basic_image)
    H = ???
    return H

So I have collected in g the images of a basis of GFX_quo_f under $\bar {Fr} - \bar {id}$. How do I compel sagemath to use this information to compute the kernel of this linear application?

EDIT 2:

Following the technical hints from @rburing below, I scrapped together this code:

def berlekamp_basis(f):
    """
    INPUT: f a monic squarefree polynomial in F_q[x]
    OUTPUT: a list of polynomials hi in F_q[x] that represent
            a basis of the Berlekamp algebra associated to f
    """
    n = f.degree()
    GFX = f.parent()
    x = GFX.gen()
    q = GFX.base().order()
    GFX_quo_f = GFX.quotient(f)
    g = []
    for k in range(n):
        basic_image = GFX_quo_f(x^(q*k)- x^k)
        g.append(vector(basic_image))
    H = matrix(g).right_kernel().basis()
    H = [GFX(h) for h in H]
    return H

When I run the following test:

GF9.<a> = GF(3^2)
GF9X.<x> = GF9[x]
f = x**3 + x**2 + x - a
assert(f.is_squarefree())
berlekamp_basis(f)

I get a big scary error:

Error in lines 42-42
Traceback (most recent call last):
  File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute
    flags=compile_flags) in namespace, locals
  File "", line 1, in <module>
  File "", line 17, in berlekamp_basis
  File "sage/structure/parent.pyx", line 921, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9679)
    return mor._call_(x)
  File "sage/structure/coerce_maps.pyx", line 145, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4574)
    raise
  File "sage/structure/coerce_maps.pyx", line 140, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4442)
    return C._element_constructor(x)
  File "/ext/sage/sage-8.4_1804/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_ring.py", line 460, in _element_constructor_
    return C(self, x, check, is_gen, construct=construct, **kwds)
  File "sage/rings/polynomial/polynomial_zz_pex.pyx", line 145, in sage.rings.polynomial.polynomial_zz_pex.Polynomial_ZZ_pEX.__init__ (build/cythonized/sage/rings/polynomial/polynomial_zz_pex.cpp:15712)
    Polynomial_template.__init__(self, parent, x, check, is_gen, construct)
  File "sage/rings/polynomial/polynomial_template.pxi", line 176, in sage.rings.polynomial.polynomial_zz_pex.Polynomial_template.__init__ (build/cythonized/sage/rings/polynomial/polynomial_zz_pex.cpp:8595)
    x = parent.base_ring()(x)
  File "sage/structure/parent.pyx", line 921, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9679)
    return mor._call_(x)
  File "sage/structure/coerce_maps.pyx", line 145, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4574)
    raise
  File "sage/structure/coerce_maps.pyx", line 140, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4442)
    return C._element_constructor(x)
  File "/ext/sage/sage-8.4_1804/local/lib/python2.7/site-packages/sage/rings/finite_rings/finite_field_givaro.py", line 370, in _element_constructor_
    return self._cache.element_from_data(e)
  File "sage/rings/finite_rings/element_givaro.pyx", line 318, in sage.rings.finite_rings.element_givaro.Cache_givaro.element_from_data (build/cythonized/sage/rings/finite_rings/element_givaro.cpp:7809)
    cpdef FiniteField_givaroElement element_from_data(self, e):
  File "sage/rings/finite_rings/element_givaro.pyx", line 407, in sage.rings.finite_rings.element_givaro.Cache_givaro.element_from_data (build/cythonized/sage/rings/finite_rings/element_givaro.cpp:6166)
    raise TypeError("e.parent must match self.vector_space")
TypeError: e.parent must match self.vector_space

I am writing my own function to solve systems of congruences in arbitrary euclidean domains, but I can't quite make it work for gaussians.

# Extended Euclides Algorithm
def extended_euclides(a,b,quo=lambda a,b:a//b):
    r0 = a; r1 = b
    s0 = 1; s1 = 0
    t0 = 0; t1 = 1

    while r1 != 0:
        q = quo(r0, r1)
        r0, r1 = r1, r0 - q * r1
        s0, s1 = s1, s0 - q * s1
        t0, t1 = t1, t0 - q * t1

    return r0, s0, t0

# Chinese Remainder Theorem
def chinese_remainder(remainders, modules, quo=lambda a,b:a//b):
    """
    INPUT:
        A list of remainders v_i and pairwise coprime modules m_i
        representing a system of congruences {x == v_i mod m_i},
        where v_i, m_i belong to an euclidean domain ED
    OUTPUT:
        A solution of the system of congruences in the domain ED
    """
    m = reduce(lambda x,y:x*y, modules)
    c = 0
    for v_i, m_i in zip(remainders, modules):
        m_div_m_i = quo(m, m_i)
        _, s_i, _ = extended_euclides(m_div_m_i, m_i, quo)
        # c_i = (v_i * s_i) % m_i
        # we need to do it convolutedly bcs gaussian integer quotient is not implemented in SageMath
        a = v_i * s_i; b = m_i
        c_i = a - quo(a,b)*b

        c += c_i * m_div_m_i
    return c

# EXAMPLE ON Z
remainders = [ZZ(2),ZZ(7)]
modules = [ZZ(11), ZZ(13)]
print("System to be solved:")
for r,m in zip(remainders, modules):
    print("x == {} (mod {})".format(r,m))
solution = chinese_remainder(remainders, modules)
print("A solution of the system is {}".format(solution))
#check that the solution is correct
for r, m in zip(remainders, modules):
    assert(solution % m == r)

# EXAMPLE ON Z[i]
ZI = QuadraticField(-1, 'I').ring_of_integers()
gaussian_quo = lambda a,b : ZI(real(a/b).round() + I*imag(a/b).round())
gaussian_rem=lambda a,b: a - gaussian_quo(a,b)*b
remainders = [ZI(2+3*I),ZI(7+5*I)]
modules = [ZI(11), ZI(13-I)]
print("System to be solved:")
for r,m in zip(remainders, modules):
    print("x == {} (mod {})".format(r,m))
solution = chinese_remainder(remainders, modules, gaussian_quo)
print("A solution of the system is {}".format(solution))
#check that the solution is correct
for r, m in zip(remainders, modules):
    print(gaussian_rem(solution, m) / r)
    assert(gaussian_rem(solution, m) == r)

This produces the following ouptut:

System to be solved:
x == 2 (mod 11)
x == 7 (mod 13)
A solution of the system is 46
System to be solved:
x == 3*I + 2 (mod 11)
x == 5*I + 7 (mod -I + 13)
A solution of the system is -41*I - 75
1
36/37*I - 6/37
Error in lines 31-33
Traceback (most recent call last):
  File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute
    flags=compile_flags) in namespace, locals
  File "", line 3, in <module>
AssertionError

So in the integer example we find a correct solution, but in the Gaussian example we find a solution that satisfies the first incongruence relation, but not the second.

I have checked with wolfram alpha, and if I did not mess up it seems like the problem is in the Chinese Remainder Theorem algorithm rather than on the assertion tests.

What is going on?