polynomial arithmetic

edit

Short solution:

p = 207661
K.<a> = GF(p^2)
print "GF(%s^2) has modulus %s" % ( p, K.modulus()  )

R.<z> = PolynomialRing( K, sparse=True )

f =     (88998*a + 148331)*z^207661 \
    - ( (29666*a + 118664)*z^207661 + (177995*a +  88998)*z )^3 \
    + ( (59332*a +  29666)*z^207661 + (148329*a + 177995)*z )^2 \
    +   (118663*a + 59333)*z - 1

g = (98319*a + 142902)*z^207661 + (109342*a + 44583)*z + 29041


def getAdHocRemainder( f, g ):
    """f,g being univariate polynomials over a common field
    compute f modulo g, the result is a polynomial of degree < g.degree()
    Note that we do *not* compute the quotient.
    """
    A = f
    gdeg = g.degree()
    G = g.monic()    # = g / g.leading_coefficient()
    z = g.parent().gens()[0]

    while A.degree() >= gdeg:
        A -= A.leading_coefficient() * G * z**(A.degree()-gdeg)

    if A:
        return A.monic()
    return A


K._gcd_univariate_polynomial = ( lambda a, b:    getAdHocRemainder( a, b ) )    # overwrite method
gcd( f, g )

Result:

z^3 + (87501*a + 22266)*z^2 + (185003*a + 55703)*z + 70580*a + 171526

Long discussion:

Let $f,g$ be the polynomials as the two posted above, here renamed. First of all, we have a quick computation in sage of the reminder of $f$ at the division with $g$, it is easy to see also mathematically with bare eyes, that we may replace modulo $g$ the big power $z^p$ in $f$ with a linear term, so the reminder must be of a low degree. In code:

sage: f % g
(82741*a + 48061)*z^3 + (165888*a + 202480)*z^2 + (83905*a + 148072)*z + 42850*a + 158093
sage: quo, rem = f.quo_rem(g)
sage: rem
(82741*a + 48061)*z^3 + (165888*a + 202480)*z^2 + (83905*a + 148072)*z + 42850*a + 158093

The quotient quo is of course also sparse, omitted here. We will meet again the method quo_rem soon.

Let us see in detail, which is the path for the computation done by sage. Calling gcd(f, g) is the same as f.gcd(g) . We can of course type f.gcd(g) and it also takes too long for the "few operations" needed, one may think. So we still have the problem, it must be there. Let us look into the code...

f.gcd??

so, first, in our case, the algorithm = "generic" path is taken. In this case, the code goes to Polynomial.gcd, so we type again Polynomial.gcd( f, g ) to convince ourselves that the problem persists. Yes, so evaluation was stopped by my Control+C in the interpreter. The code of Polynomial.gcd can be obtained in the sage interpreter via

sage: ??Polynomial.gcd

It delivers among many other rows:

    if hasattr(self.base_ring(), '_gcd_univariate_polynomial'):
        return self.base_ring()._gcd_univariate_polynomial(self, other)

And it also contains in the doc string the

    """
    .. NOTE::

        The actual algorithm for computing greatest common divisors depends
        on the base ring underlying the polynomial ring. If the base ring
        defines a method ``_gcd_univariate_polynomial``, then this method
        will be called (see examples below).
    """

So we further go to

sage: f.base_ring()
Finite Field in a of size 207661^2
sage: K._gcd_univariate_polynomial??

and among the many lines we get the place that hurts:

q, r = a.quo_rem(b)

Yes, the quotient q is explicitly stored in the second we call the method.

Note: There is a similar problem over R, the polynomial ring over K, and the implemented method obtained via R._gcd_univariate_polynomial?? has an other information in the doc string:

        """
        ALGORITHM:

        Algorithm 3.3.1 in [Coh1993]_, based on pseudo-division.
        """

where we see the calls of the method pseudo_quo_rem . A cousin of quo_rem. We also go into this method, and see that it computes and returns also the (unneeded) quotient. In our case, the first step works in...

sage: _, rem = f.quo_rem(g)
sage: rem
(82741*a + 48061)*z^3 + (165888*a + 202480)*z^2 + (83905*a + 148072)*z + 42850*a + 158093

but

sage: _, rem = g.quo_rem(rem)    # stopped after some time, Control+C

Also (and much slower):

sage: _, rem = f.pseudo_quo_rem(g)
sage: rem
(18510*a + 199577)*z^3 + (191808*a + 60662)*z^2 + (167956*a + 170038)*z + 60490*a + 63076

but

sage: _, rem = g.pseudo_quo_rem(rem)    # stopped after some time, Control+C

So computing the reminder of the division of $g$, a polynomial of degree $p=207661$, with the degree three polynomial rem, takes a loooong time already. Yes, it is not sparse, so only trying to store it in between is a difficult task.

So this is the reason. This is of course not satisfying the computational needs. Can we implement this better, for this one purpose, in our own class? Yes. This would be the gentle way.

The quick answer goes the dirty way, it uses the cheap function, getAdHocRemainder, which is not optimized for speed: But it works. With it and with $f,g$ as above we have after a not so long time:

sage: g
(98319*a + 142902)*z^207661 + (109342*a + 44583)*z + 29041
sage: rem = getAdHocRemainder( f,g )
sage: rem
z^3 + (87501*a + 22266)*z^2 + (185003*a + 55703)*z + 70580*a + 171526
sage: ( f % g ).monic()
z^3 + (87501*a + 22266)*z^2 + (185003*a + 55703)*z + 70580*a + 171526
sage: rem2 = getAdHocRemainder( g, rem )
sage: rem2
0

So the above rem is the g.c.d of $f,g$. We could stop here, but...

How can we test this quickly (in an other way)? For instance by getting the extension field of K = GF(p^2) by a root $w$ of rem, then testing if $f$ and $g$ vanish in $w$ or not.

# p = 207661
# K.<a> = GF(p^2)
r = rem    # use it as modulus for L
L.<w> = K.extension( r )
RL.<Z> = PolynomialRing( L )
rL = RL(r)
print "rL    = %s" % rL
print "rL(w) = %s" % rL(w)
gL = RL(g)
print "gL    = %s" % gL
print "gL(w) = %s" % gL(w)

And we get, after some time:

rL    = Z^3 + (87501*a + 22266)*Z^2 + (185003*a + 55703)*Z + 70580*a + 171526
rL(w) = 0
gL    = (98319*a + 142902)*Z^207661 + (109342*a + 44583)*Z + 29041
gL(w) = 0

I should maybe stop here. If any questions, please do not hesitate!