Ask Your Question
0

Subresultant algorithm taking a lot of time for higher degree univariate polynomials with coefficients from fraction fields

asked 2017-05-23 10:50:37 +0100

mua gravatar image

updated 2017-05-24 10:06:19 +0100

I have to compute the gcd of univariate polynomials over the fraction field of $\mathbb{Z}[x,y]$. I wanted to use the subresultant algorithm already implemented for UFDs. I copied the same function to fraction_field.py. The subresultant algorithm calls the psuedo division algorithm which has the following step : R = dR - cB.shift(diffdeg) - this hangs when we consider random polynomials of degree >6 in $Frac(\mathbb{Z}[x,y])[z]$.

(Note: In the current version of sage it uses the regular Euclidean algorithm implemented in rings.py for computing gcd in this case. It is much slower than the subresultant algorithm (hangs for degrees >4) which is why I thought the subresultant algorithm will improve things.)

Sample input: sage: A.<x,y>=ZZ[]

sage: B= Frac(A)

sage: C.<z> = B[]

sage: p = C.random_element(6)

sage: q = C.random_element(6)

sage: gcd(p,q)

The following function is what I copied into fraction_field.py from unique_factorisation_domain.py.

def _gcd_univariate_polynomial(self, f, g):

    if f.degree() < g.degree():
        A,B = g, f
    else:
        A,B = f, g

    if B.is_zero():
        return A

    a = b = self.zero()
    for c in A.coefficients():
        a = a.gcd(c)
        if a.is_one():
            break
    for c in B.coefficients():
        b = b.gcd(c)
        if b.is_one():
            break

    d = a.gcd(b)

    #d=1
    A = A // a
    B = B // b
    g = h = 1

    delta = A.degree()-B.degree()
    _,R = A.pseudo_quo_rem(B)

    while R.degree() > 0:
        A = B
        B = R // (g*h**delta)
        g = A.leading_coefficient()
        h = h*g**delta // h**delta
        delta = A.degree() - B.degree()
        _, R = A.pseudo_quo_rem(B)
       # print("i am here")


    if R.is_zero():
        b = self.zero()
        for c in B.coefficients():
            b = b.gcd(c)
            if b.is_one():
                break

        return d*B // b

    return d

This calls the following pseudo quo remainder function in polynomial_element.pyx. It is in this function I was able to see that it hangs at R = dR - cB.shift(diffdeg).

def pseudo_quo_rem(self,other):

    if other.is_zero():
        raise ZeroDivisionError("Pseudo-division by zero is not possible")

    # if other is a constant, then R = 0 and Q = self * other^(deg(self))
    if other in self.parent().base_ring():
        return (self * other**(self.degree()), self._parent.zero())

    R = self
    B = other
    Q = self._parent.zero()
    e = self.degree() - other.degree() + 1
    d = B.leading_coefficient()

    while not R.degree() < B.degree():

        c = R.leading_coefficient()

        diffdeg = R.degree() - B.degree()

        Q = d*Q + self.parent()(c).shift(diffdeg)

        R = d*R - c*B.shift(diffdeg)

        e -= 1

    q = d**e
    return (q*Q,q*R)
edit retag flag offensive close merge delete

Comments

1

It is really hard to reconstruct the error. (Without also copying the obvious same function to fraction_field, and guessing the code that builds the fraction fields, and knowing the meaning of d,R, cB, shift, diffdeg after an fgrep of potential modules, and searching / finding an explicit example that reproduces the same problem.) It would be simpler for a potential helper if the question already isolates a minimal example with simple polynomials also leading to the same error, and the code would be welcome...

dan_fulea gravatar imagedan_fulea ( 2017-05-23 18:48:27 +0100 )edit

yes I understand that. I will edit the question with more details.

mua gravatar imagemua ( 2017-05-23 22:54:17 +0100 )edit

I am the author of these lines. I don't have time to have a careful look right now, but I'll do that next week.

B r u n o gravatar imageB r u n o ( 2017-05-24 10:32:54 +0100 )edit

thank you! that will be great!

mua gravatar imagemua ( 2017-05-24 10:38:26 +0100 )edit

1 Answer

Sort by ยป oldest newest most voted
2

answered 2017-06-01 15:39:32 +0100

B r u n o gravatar image

I looked into the issue. You correctly identified the place where it hangs. Actually, in the line R = d*R - c*B.shift(diffdeg), the long time is taken by the two multiplications.

The problem is that d and the coefficients of R are very large rational fractions: I stopped the computation at some point when it hanged, and stored the values of R and d there: The degrees in each variable is around 50 for both d and the coefficients of R, and their number of terms is around 2500 each. Already computing the product of d with one coefficient of R takes around 2s on the machine I use.

Looking more into details, the problem is SageMath's multiplication of bivariate polynomials which is quite slow. To sum up: I think there is (unfortunately) no bug, simply the objects are too large for the current capabilities of SageMath.

As a further comment, the value of d and the coefficients of R are all polynomials in $\mathbb Z[x,y]$ (they are rational fractions, but with denominator $1$). Putting them in the ring A (that represents $\mathbb Z[x,y]$ in your example) does not change anything on the computation time. On the other hand, using as a base ring $\mathbb Z[y][x]$ improves the performance of the multiplication (the computation is approx. 20 times faster). One solution may then be to modify pseudo_quo_rem for your needs, using a recursive univariate polynomial ring (U = ZZ['y']['x']) as base ring instead of a multivariate polynomial ring.

P.S.: Though this website is devoted to SageMath, I think that the (free, of course) Computer Algebra System Nemo, written is Julia, is particularly efficient for such kind of computations. You may have a look to it.

edit flag offensive delete link more

Comments

thank you so much for the reply. We are building some packages for sagemath and this gcd is getting called in them. Can you probably explain what you mean by the limitations of sage math - is it for example because of the limitations of singular (Since in the multiplication function it calls the gcd function in multi_polynomial_libsingular.pyx and always hangs there since the number of monomials is big) ? or is it something else in the architecture? The reason I was asking is because I was wondering if it can be ``rectified" if we spend more time on it? or would it be a complete overhaul of the code?

mua gravatar imagemua ( 2017-06-01 20:48:25 +0100 )edit

Also when i tried adding a dummy variable in the polynomial ring and type casting the polynomials p and q as multivariate polynomials it was computed very fast. Any reason for why it happens differently when it is a multivariate polynomial over a fraction field versus a univariate polynomial over a fraction field?

mua gravatar imagemua ( 2017-06-02 17:12:39 +0100 )edit

Multivariate polynomials are handled by Singular while univariate polynomials are handled by Flint or NTL. This explains that computation times may vary a lot between univariate and multivariate implementations since different libraries are doing the actual work behind the scene. One would need to delve further to fully understand the differences in the computation times between univariate and multivariate polynomials. Note that you do not need to add a dummy variable to use the multivariate implementation. You can define your polynomial ring as R.<x> = PolynomialRing(A, 1) to get a multivariate polynomial ring in one variable.

B r u n o gravatar imageB r u n o ( 2017-06-06 11:04:05 +0100 )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

Stats

Asked: 2017-05-23 10:50:37 +0100

Seen: 1,162 times

Last updated: Jun 01 '17