Ask Your Question
1

polynomial arithmetic

asked 2017-10-08 01:41:41 -0600

this post is marked as community wiki

This post is a wiki. Anyone with karma >750 is welcome to improve it.

i have two polynomial over polynomial ring z . I am facing problem in computation of gcd of two polynomial over polynomial ring z.

p=207661
K.<a>=GF(p^2);K.modulus() 
R.<z> = PolynomialRing( K, sparse=True )
fE=(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;fE
fL= (98319*a + 142902)*z^207661 + (109342*a + 44583)*z + 29041;fL
f1=gcd(fE,fL)

required lot of time in the computation of gcd.

edit retag flag offensive close merge delete

Comments

Please provide the code correctly indented, with multiplication explicitly written if possible, and with matching parenthesis.

Else it is hard to start.

Use that button with 101 and 010 to shift 4 spaces in each row the marked code fragment. This makes a markdown, so that the message is correctly displayed on the page.

dan_fulea gravatar imagedan_fulea ( 2017-10-08 03:01:13 -0600 )edit

1 answer

Sort by ยป oldest newest most voted
1

answered 2017-10-09 16:17:05 -0600

dan_fulea gravatar image

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!

edit flag offensive delete link more

Comments

thanks a lot

santoshi gravatar imagesantoshi ( 2017-10-10 04:17:24 -0600 )edit

@dan_fulea great, perhaps should you consider implement your solution within Sage source code (computing the reminder without the quotient).

tmonteil gravatar imagetmonteil ( 2017-10-10 04:58:20 -0600 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2017-10-08 01:41:41 -0600

Seen: 108 times

Last updated: Oct 09