gcd can be very slow in AA[x]

[vdelecroix]

Here is a x50 slower example

sage: x,y = polygens(QQ,"x,y")
sage: p1 = x^5 + 6*x^4 - 42*x^3 - 142*x^2 + 467*x + 422
sage: p2 = p1(x=(x-1)^2)
sage: p3 = p2(x=x*y).resultant(p2,x).univariate_polynomial()
sage: p4, = [f[0] for f in p3.factor() if f[0].degree() == 80]
sage: %time _ = p4.squarefree_decomposition()
CPU times: user 807 µs, sys: 0 ns, total: 807 µs
Wall time: 883 µs
sage: %time _ = p4.change_ring(AA).squarefree_decomposition()
CPU times: user 40.1 s, sys: 3.21 ms, total: 40.1 s
Wall time: 40.1 s

This problem originally appeared in #17895 (where a better workaround has been found).

[git]

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

​62db08c

28199: Try to use subresultants algorithm

[bruno]

I've made a very quick try: Instead of using the plain old Euclidean algorithm, one can use the so-called subresultant algorithm. This algorithm is already provided for UFDs so I tried to branch this algorithms for fields also. The resulting timings are still bad, but slightly less so:

sage: x,y = polygens(QQ,"x,y")
sage: p1 = x^5 + 6*x^4 - 42*x^3 - 142*x^2 + 467*x + 422
sage: p2 = p1(x=(x-1)^2)
sage: p3 = p2(x=x*y).resultant(p2,x).univariate_polynomial()
sage: p4, = [f[0] for f in p3.factor() if f[0].degree() == 80]
sage: %time _ = p4.squarefree_decomposition()
CPU times: user 1.89 ms, sys: 0 ns, total: 1.89 ms
Wall time: 1.91 ms
sage: %time _ = p4.change_ring(AA).squarefree_decomposition()
CPU times: user 3.14 s, sys: 15.3 ms, total: 3.16 s
Wall time: 3.16 s

[vdelecroix]

You were quick :-) This is great: 10x improvement! There is no hope to beat the C-implementation from flint anyway.

Do you know why the subresultant is not the default method over fields? Could it be always better?

This reminds me the problem for computing determinant of matrices with rational coefficients. Going via Gauss elimination (which is a sort of "default" for fields) will create a terrible coefficient blowup.

[bruno]

Replying to vdelecroix:

You were quick :-) This is great: 10x improvement! There is no hope to beat the C-implementation from flint anyway.

Of course, but I would hope for less than a second...

Do you know why the subresultant is not the default method over fields? Could it be always better?

That's why my first try. The problem is that it requires a _floordiv_ and not all fields have it. For instance:

sage: R.<x> = RR[]
sage: p = (x^2-1)^2
sage: p.is_squarefree()
Traceback (most recent call last):
...
TypeError: unsupported operand parent(s) for //: 'Real Field with 53 bits of precision' and 'Real Field with 53 bits of precision'

or

sage: S.<y> = GF(5)[]
sage: Q = S.quotient_ring(y^2+y+1)
sage: T.<z> = Q[]
sage: gcd(y*z+1, z + y)
Traceback (most recent call last):
...
TypeError: unsupported operand parent(s) for //: 'Univariate Quotient Polynomial Ring in ybar over Finite Field of size 5 with modulus y^2 + y + 1' and 'Univariate Quotient Polynomial Ring in ybar over Finite Field of size 5 with modulus y^2 + y + 1'

One possibility to check is whether adding a _floordiv_ method that simply calls self._div_(other) in the Fields category is sufficient and yields better timings.

[bruno]

Actually I checked: It is not faster in general (you can for instance test in the ring T of my example of the previous comment). This is reasonable: The only cases where it is faster is when there is coefficient growth. So there is a difficulty to decide when to use each algorithm. Maybe I should add a test as: if not self.base_ring().is_finite(). Note that in that case, RR would use subresultants while it is slightly slower. (But not that computing (non-approximate) gcds for polynomials over RR is mainly nonsense because of approximation issues.)

[vdelecroix]

I think that exact algorithm should raise an error with non-approximate base ring.

[embray]

Ticket retargeted after milestone closed

[mkoeppe]

Batch modifying tickets that will likely not be ready for 9.1, based on a review of the ticket title, branch/review status, and last modification date.

[mkoeppe]

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.