Opened 5 years ago

Closed 4 years ago

#15790 closed defect (fixed)

GCD of univariate polynomials: generic implementation for UFD, and sparse case

Reported by: bruno Owned by:
Priority: minor Milestone: sage-6.9
Component: basic arithmetic Keywords: polynomial, gcd
Cc: vdelecroix Merged in:
Authors: Bruno Grenet Reviewers: Julian Rueth
Report Upstream: N/A Work issues:
Branch: d77ee90 (Commits) Commit: d77ee902e808082dcc44daf71bb0ea2756c621bb
Dependencies: Stopgaps:

Description (last modified by bruno)

[Update] Due to #13442, the purpose of this ticket changed slightly. It provides two functionalities:

  1. A generic implementation for the GCD of two polynomials over a UFD;
  2. A gcd method for sparse polynomials which allows to choose between the generic implementation, or using the implementation for dense polynomials. Default is generally generic, but for polynomial over ZZ it is dense. This is due to the very fast implementations available in Flint and NTL.

These implementations solve the problem mentioned in the first version of this ticket, that remains below for information.


[Old description, not relevant anymore]

The GCD of two sparse polynomials over ZZ is not computed. If the polynomials are in dense representation or if they are in sparse representation but defined over QQ, there is no problem.

sage: R1.<x>=PolynomialRing(ZZ,'x',sparse=true)
sage: R2.<x>=PolynomialRing(ZZ,'x',sparse=false)
sage: R3.<x>=PolynomialRing(QQ,'x',sparse=true)
sage: p=x^2-3*x+2
sage: q=x^2-1
sage: gcd(R2(p),R2(q))
x - 1
sage: gcd(R3(p),R3(q))
x - 1
sage: gcd(R1(p),R1(q))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-6-0d7380207fef> in <module>()
----> 1 gcd(R1(p),R1(q))

/home/bgrenet/bin/sage-6.0-x86_64-Linux/local/lib/python2.7/site-packages/sage/rings/arith.pyc in gcd(a, b, **kwargs)
   1605                 return ZZ(a).gcd(ZZ(b))
   1606             except TypeError:
-> 1607                 raise TypeError("unable to find gcd")
   1608         return GCD(b, **kwargs)
   1609 

TypeError: unable to find gcd

Change History (29)

comment:1 Changed 5 years ago by bruno

For PolynomialRing(ZZ,'x',sparse=false), one of FLINT or NTL is used (FLINT by default).

For PolynomialRing(QQ,'x',sparse=true), the gcd computation is based on quo_rem. This is not possible over the integer since quo_rem works for polynomials over a field.

There are several possible solutions:

  1. Coerce the polynomials to PolynomialRing(QQ,'x',sparse=true); compute their GCD g, the LCM l of the denominators of the coefficients of g, and return l*g.
  2. Coerce the polynomials to PolynomialRing(ZZ,'x',sparse=false); compute their GCD and return the result in sparse representation with a second coercion.
  3. Implement a pseudo-division algorithm (cf. Cohen GTM 138 Algorithm 3.1.2) and a GCD algorithm based on the sub-resultant (cf. Cohen GTM 138 Algorithm 3.3.1).

Solution 3 seems cleaner but does not give especially efficient results since the algorithms described in Cohen GTM 138 have complexity polynomial in the degree of the input polynomial. Actually, computing the GCD of two sparse polynomial is NP-hard ![1] so no better solution than converting to the dense representation is likely to exist. To my mind, Solution 2 is then the best one.

Yet, I am new to Sage so comments are welcome!

![1] D. Plaisted, Sparse complex polynomials and polynomial reducibility. J. Comput. Syst. Sci. 14 (2), 210–221, 1977.

comment:2 Changed 5 years ago by tscrim

I'd implement all 3 with an optional algorithm argument since each has different strengths/weaknesses. For example, 2 is good as long as the polynomials are relatively dense, but for something like x^100000 - 2, I wouldn't want to convert that to a dense polynomial because of the memory usage (and take the time to construct it). From that, I think 3 should be the default. Possible names:

  1. "fraction_field"
  2. "dense"
  3. "pseudo-division"

Nevertheless, I think we should have some gcd since this would be a surprise for the more casual user.

comment:3 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:4 Changed 5 years ago by vdelecroix

  • Cc vdelecroix added

comment:5 Changed 5 years ago by bruno

  • Branch set to u/bruno/gcd_of_sparse_univariate_polynomials_over_zz

comment:6 Changed 5 years ago by git

  • Commit set to 8ddbb9efad8d53552fa26727fb484720bda66025

Branch pushed to git repo; I updated commit sha1. New commits:

8ddbb9eAdd gcd method to class Polynomial_generic_sparse

comment:7 Changed 5 years ago by bruno

I've implemented the gcd method in the class Polynomial using the pseudo-division algorithm, as well as a method gcd in Polynomial_generic_sparse which allows to choose the algorithm between "pseudo-division", "dense" or "fraction_field".

Needs review!

One remaining limitation is the absence of quo_rem method in Polynomial_generic_sparse. I created a specific ticket #16544 for this.

comment:8 Changed 5 years ago by bruno

  • Keywords polynomial gcd added; GCD removed
  • Status changed from new to needs_review

comment:9 Changed 5 years ago by bruno

  • Authors set to Bruno Grenet

comment:10 Changed 5 years ago by git

  • Commit changed from 8ddbb9efad8d53552fa26727fb484720bda66025 to 65943e1a4ae6da9d2c7579ba14317aeada6eb102

Branch pushed to git repo; I updated commit sha1. New commits:

96f11a1Add special case for p.pseudo_quo_rem(q) when q is a constant
9000fb8Add quo_rem method in class Polynomial_generic_sparse
84833c9Move inversion of unit
fb1d21cAdd doc tests
65943e1Merge branch 't/16544/fb1d21cbee05d6be577bf66120aa20d3bd498ace' into t/15790/gcd_of_sparse_univariate_polynomials_over_zz

comment:11 Changed 5 years ago by bruno

  • Status changed from needs_review to needs_work

comment:12 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:13 Changed 4 years ago by git

  • Commit changed from 65943e1a4ae6da9d2c7579ba14317aeada6eb102 to 5c84d7115a1e5406487b5d59746554b44be995e4

Branch pushed to git repo; I updated commit sha1. New commits:

bf5a846Correct formatting of doctests + raise error for unknown algorithm
5c84d71Change Error for pseudo-division by zero

comment:14 Changed 4 years ago by bruno

  • Status changed from needs_work to needs_review

The ticket is ready for review. Note that as default, I have the "dense" algorithm for polynomials over ZZ and "pseudo-division" for the other ones. The idea is that arithmetic on ZZ[X] is very efficient.

Of course, gcd(x^100000-1,x^1000-1) is much faster (by a factor approx. 400) with "pseudo-division" algorithm because of memory allocation for the "dense" algorithm, and thanks to a very good behavior of the "pseudo-division" algorithm on this example. Yet, if you consider gcd(x^100000-1,x^2-1) (which should be easy btw), the "dense" algorithm is faster and I even stopped the "pseudo-division" algorithm after a quite long time.

Note: It may be possible to reduce the huge allocation time since with NTL implementation, this is much faster: https://groups.google.com/forum/#!topic/sage-devel/6qhW90dgd1k.

Last edited 4 years ago by bruno (previous) (diff)

comment:15 follow-up: Changed 4 years ago by tscrim

Do you have a heuristic for determining when one algorithm is faster than another (when the base ring is ZZ)? It doesn't have to be a particularly good one, but considering the time differences you're getting, it seems like we should have some automatic choice of algorithm.

comment:16 Changed 4 years ago by git

  • Commit changed from 5c84d7115a1e5406487b5d59746554b44be995e4 to ec2ef11336040e95c8cdf1bed7523006747b5ab1

Branch pushed to git repo; I updated commit sha1. New commits:

ec2ef11Choice between FLINT and NTL

comment:17 in reply to: ↑ 15 Changed 4 years ago by bruno

Replying to tscrim:

Do you have a heuristic for determining when one algorithm is faster than another (when the base ring is ZZ)? It doesn't have to be a particularly good one, but considering the time differences you're getting, it seems like we should have some automatic choice of algorithm.

I don't have one. Based on the comments on https://groups.google.com/forum/#!topic/sage-devel/6qhW90dgd1k, I add the possibility¹ for the "dense" algorithm to use NTL instead of FLINT. After some tests, I chose a heuristic for using NTL rather than FLINT (basically, FLINT is faster for low-degree polynomials and not-too-sparse ones). Then the "dense" algorithm is almost always faster than "pseudo-division" but for some very specific examples (I haven't been able to build other ones): gcd(x^N-1,x^M-1) with both N and M pretty large. An the ratio is not hundreds anymore, but approx. 10.

¹ Actually, no choice for the user but only an automatic switch.

Last edited 4 years ago by bruno (previous) (diff)

comment:18 Changed 4 years ago by chapoton

  • Status changed from needs_review to needs_work

needs work, does not build, see patchbot report:

File "sage/structure/element.pyx", line 3433, in sage.structure.element.NamedBinopMethod.__init__ (build/cythonized/sage/structure/element.c:26905)
AttributeError: 'sage.structure.element.NamedBinopMethod' object has no attribute '__name__'

comment:19 Changed 4 years ago by git

  • Commit changed from ec2ef11336040e95c8cdf1bed7523006747b5ab1 to 3015456844cb54f37cf4bbae8423521cdd33624a

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

b272960Add pseudo-division and gcd for Polynomials (as generic as possible)
c88b600Add gcd method to class Polynomial_generic_sparse
6f6a167Add special case for p.pseudo_quo_rem(q) when q is a constant
9f797efCorrect formatting of doctests + raise error for unknown algorithm
6c31b67Change Error for pseudo-division by zero
f623b88Choice between FLINT and NTL
3015456Adapted to 6.9beta1 => Merged two gcd methods

comment:20 Changed 4 years ago by bruno

  • Description modified (diff)
  • Milestone changed from sage-6.4 to sage-6.9
  • Status changed from needs_work to needs_review
  • Summary changed from GCD of sparse univariate polynomials over ZZ to GCD of [sparse] univariate polynomials [over ZZ] [updated]

This is rebased, and now works (as far as I can tell).

comment:21 Changed 4 years ago by chapoton

  • Status changed from needs_review to needs_work

Doc does not build

OSError: [polynomia] docstring of sage.rings.polynomial.polynomial_element.Polynomial.pseudo_quo_rem:38: WARNING: Duplicate explicit target name: "gtm138".

comment:22 Changed 4 years ago by git

  • Commit changed from 3015456844cb54f37cf4bbae8423521cdd33624a to 47badc1b7322077ffbf2f89ab6ad2acc2dce18b4

Branch pushed to git repo; I updated commit sha1. New commits:

47badc1Move fallback algorithm to category UFD

comment:23 Changed 4 years ago by git

  • Commit changed from 47badc1b7322077ffbf2f89ab6ad2acc2dce18b4 to f9089369482cc589372b0f3a7461d633d9540964

Branch pushed to git repo; I updated commit sha1. New commits:

f908936Simplify GCD for sparse polynomials

comment:24 Changed 4 years ago by bruno

  • Description modified (diff)
  • Keywords sparse removed
  • Status changed from needs_work to needs_review
  • Summary changed from GCD of [sparse] univariate polynomials [over ZZ] [updated] to GCD of univariate polynomials: generic implementation for UFD, and sparse case

In the two latest commits, I propose the following changes:

  1. Instead of having a fallback algorithm in the method Polynomial.gcd which tests whether the base ring in a unique factorization domain, this algorithm is moved to the category UniqueFactorizationDomains.
  2. I simplified the method Polynomial_generic_sparse.gcd: There remain only two possibilities (out of three), namely the one previously called pseudo-division (renamed generic, since it uses the generic code) and the dense version, which may in some cases be faster, for instance over ZZ. I removed the fraction_field option which does not make much sense I think.

At the same time, I also take Frédéric's comment into account.

comment:25 Changed 4 years ago by saraedum

  • Branch changed from u/bruno/gcd_of_sparse_univariate_polynomials_over_zz to u/saraedum/gcd_of_sparse_univariate_polynomials_over_zz

comment:26 Changed 4 years ago by git

  • Commit changed from f9089369482cc589372b0f3a7461d633d9540964 to d77ee902e808082dcc44daf71bb0ea2756c621bb

Branch pushed to git repo; I updated commit sha1. New commits:

d77ee90Improved docstring of PolynomialElement.pseudo_quo_rem()

comment:27 Changed 4 years ago by saraedum

  • Status changed from needs_review to positive_review

Made some minor docstring changes. Everything else looks fine.

comment:28 Changed 4 years ago by bruno

  • Reviewers set to Julian Rueth

Thanks for improving the docstring, and for the review! I add your name as a reviewer.

comment:29 Changed 4 years ago by vbraun

  • Branch changed from u/saraedum/gcd_of_sparse_univariate_polynomials_over_zz to d77ee902e808082dcc44daf71bb0ea2756c621bb
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.