Opened 5 years ago

Closed 4 years ago

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

Reported by: Owned by: bruno minor sage-6.9 basic arithmetic polynomial, gcd vdelecroix Bruno Grenet Julian Rueth N/A d77ee90 (Commits) d77ee902e808082dcc44daf71bb0ea2756c621bb

[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
```

### 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: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

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

 ​8ddbb9e `Add 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

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

 ​96f11a1 `Add special case for p.pseudo_quo_rem(q) when q is a constant` ​9000fb8 `Add quo_rem method in class Polynomial_generic_sparse` ​84833c9 `Move inversion of unit` ​fb1d21c `Add doc tests` ​65943e1 `Merge 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:

 ​bf5a846 `Correct formatting of doctests + raise error for unknown algorithm` ​5c84d71 `Change 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: ↓ 17 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:

 ​ec2ef11 `Choice between FLINT and NTL`

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

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:

 ​b272960 `Add pseudo-division and gcd for Polynomials (as generic as possible)` ​c88b600 `Add gcd method to class Polynomial_generic_sparse` ​6f6a167 `Add special case for p.pseudo_quo_rem(q) when q is a constant` ​9f797ef `Correct formatting of doctests + raise error for unknown algorithm` ​6c31b67 `Change Error for pseudo-division by zero` ​f623b88 `Choice between FLINT and NTL` ​3015456 `Adapted 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

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

 ​47badc1 `Move fallback algorithm to category UFD`

### comment:23 Changed 4 years ago by git

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

 ​f908936 `Simplify 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:

 ​d77ee90 `Improved 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.