`Fraction must have unit denominator`when using polynomial.numerator()

asked 2018-11-20 14:44:51 +0200

Jsevillamol gravatar image

updated 2018-11-20 14:45:25 +0200

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?

edit retag flag offensive close merge delete

Comments

You should learn to debug a little better, at least to pin down where the error comes from. In this instance e.g. search for division / in your code. Here it seems the error comes when you try to convert g/h to a polynomial with integer coefficients. This would work if h was a factor of g, but it seems in your example not to be the case, so you are calculating h in a wrong way.

rburing gravatar imagerburing ( 2018-11-20 17:25:10 +0200 )edit