avoid unnecessary divisions and calls to gcd

[mantepse]

In an effort to make #38108 more usable, we optimize the computation of gcd's in generic polynomial rings.

[github-actions]

Documentation preview for this PR (built with commit 46eddc8; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[mantepse]

This reduces the time for the example from #38108 and some others significantly. The old timing is

sage: h = SymmetricFunctions(QQ).h()
sage: L.<t, u> = LazyPowerSeriesRing(h.fraction_field())
sage: D = L.undefined()
sage: s1 = L.sum(lambda n: h[n]*t^(n+1)*u^(n-1), 1)
sage: L.define_implicitly([D], [u*D - u - u*s1*D - t*(D - D(t, 0))])
sage: D[0]
h[]
sage: %time repr(D)
CPU times: user 28.6 s, sys: 8.55 ms, total: 28.6 s
Wall time: 28.7 s
'h[] + h[1]*t^2 + ((h[1,1]+h[2])*t^4+h[2]*t^3*u) + ((h[1,1,1]+3*h[2,1]+h[3])*t^6+(2*h[2,1]+h[3])*t^5*u+h[3]*t^4*u^2) + O(t,u)^7'                             

The new timing is about 5 seconds. With a more aggressive patch we could reduce it to 3.6 seconds.

sage: P.<x,y>=ProjectiveSpace(QQbar, 1)
sage: E=EllipticCurve([1, 2])
sage: f=P.Lattes_map(E, 2)
sage: %time f.Lattes_to_curve(check_lattes=true)
CPU times: user 24.5 s, sys: 156 ms, total: 24.7 s
Wall time: 24.9 s
Elliptic Curve defined by y^2 = x^3 + x + 2 over Rational Field

The new timing is about 15 seconds.

The drawback is that some computations give less beautiful (but of course, equivalent) results. For example:

sage: R.<a,b> = NumberField(x^2 - 3, 'g').extension(x^2 - 7, 'h')[]
sage: h = R.base_ring().gen()
sage: S.<y> = R.fraction_field()[]
sage: xgcd(y^2, a*h*y + b)
(1, 7*a^2/b^2, (((-7)*a)/(h*b^2))*y + 7/(7*b))

instead of (1, 7*a^2/b^2, (((-h)*a)/b^2)*y + 1/b).

[mantepse]

The failures are a bit puzzling. In arith/misc.py:

        sage: # needs sage.rings.number_field
        sage: R.<a,b> = NumberField(x^2 - 3, 'g').extension(x^2 - 7, 'h')[]
        sage: h = R.base_ring().gen()
        sage: S.<y> = R.fraction_field()[]
        sage: xgcd(y^2, a*h*y + b)
        (1, 7*a^2/b^2, (((-h)*a)/b^2)*y + 1/b)

which gives with this patch

(1, 7*a^2/b^2, (((-7)*a)/(h*b^2))*y + 7/(7*b))

These are the same elements, just represented differently.

In schemes/affine/affine_morphism.py:

            sage: A.<z> = AffineSpace(QQbar, 1)
            sage: H = End(A)
            sage: f = H([2*z / (z^2 + 2*z + 3)])
            sage: f.homogenize(1)
            Scheme endomorphism of Projective Space of dimension 1
             over Algebraic Field
              Defn: Defined on coordinates by sending (x0 : x1) to
                    (x0*x1 : 1/2*x0^2 + x0*x1 + 3/2*x1^2)

we now get

Scheme endomorphism of Projective Space of dimension 1 over Algebraic Field
  Defn: Defined on coordinates by sending (x0 : x1) to
        (2*x0*x1 : x0^2 + 2*x0*x1 + 3*x1^2)

I don't know, whether this is really the same thing.

I guess they come from the fact that gcd(0, a) may differ from a by a unit.

[tscrim]

Yes, that is the same morphism as it is projective space (and thus you can scale the coordinates freely).

[tscrim]

I am pretty sure that b divides d in the last part because d is the gcd of the coefficients of the initial A and B, which then we reduce B by some common factor with A and then get the gcd of the resulting coefficients. So it should be safe to do return (d // b) * B. Although perhaps we should cc someone who is more of an expert to triple-check this.

[mantepse]

It needs work, because I want to incorporate the much faster version I now have, and I think that the current version has some weak spots. For example, gcd is not required to return 1 for units - at least this is not done in several cases, so testing whether the result equals one is a bit weak.

While my implementation is very much faster on some examples, there are others where it leads to problems, and I have to investigate this.

[tscrim]

It needs work, because I want to incorporate the much faster version I now have, and I think that the current version has some weak spots. For example, gcd is not required to return 1 for units - at least this is not done in several cases, so testing whether the result equals one is a bit weak.

Indeed, the gcd is only defined up to a unit.

While my implementation is very much faster on some examples, there are others where it leads to problems, and I have to investigate this.

It would be good with your push to add some tests for these cases (at least, given the bot results, it seems like they are not already included).

[enriqueartal]

It is maybe unrelated but it is possible to compute the gcd for a list of rational functions. I encountered a case where the computation takes long (I do not know how long since I stopped it) but working with the gcd's (and maybe lcm's) of numerators and denominators, it becomes very fast.

[enriqueartal]

Sorry for the noise, lcm is slow.

[mantepse]

I think the patch as is is quite good, so setting to needs-review, I will also advertise it on sage-devel. The results are not as beautiful as with the current version, but it is significantly faster.

I have an additional improvement, which I am not so sure about. It avoids constructing a new polynomial ring for each variable. Instead, it constructs a univariate ring over the polynomial ring once, and reuses this until all coefficients are constants. I am sure this could be further improved and streamlined, but it is unclear whether it is worth the effort.

[mantepse]

I am pretty sure that b divides d in the last part because d is the gcd of the coefficients of the initial A and B, which then we reduce B by some common factor with A and then get the gcd of the resulting coefficients. So it should be safe to do return (d // b) * B. Although perhaps we should cc someone who is more of an expert to triple-check this.

This turns out not to be the case:

sage: B.<X, Y> = QQ[]
sage: C.<f> = B[]
sage: p = X*f^3 + X^2*Y*f
sage: q = X*Y*f^2
sage: C._gcd_univariate_polynomial(p, q)

At the end of the algorithm, B = X*Y*f, b = X*Y and d = X, so b does not divide d.

I did not yet check whether this is to be expected.

[mantepse]

I actually do not quite understand what happens concerning the beauty of results. For example, also in develop we have

sage: R.<x,y> = QQbar[]
sage: 7*x / (7*y)
7*x/(7*y)

[tscrim]

That's a good example. Can you add it as a test and a note in the code?

I don't know why that isn't reducing itself properly. That would be a bug IMO.

[mantepse]

That's a good example. Can you add it as a test and a note in the code?

I'm not sure where, but yes.

I don't know why that isn't reducing itself properly. That would be a bug IMO.

I don't think it's a bug: the gcd is one:

sage: a = 7*x; b = 7*y
sage: gcd(a, b)
1

We could require fractions of polynomials to have monic leading terms, but I don't think we always want that:

sage: R.<r,s> = QQ[]
sage: (3*r+s)/(2*r+s)
(3*r + s)/(2*r + s)

[mantepse]

In fact:

sage: R.<r,s> = QQ[]
sage: (3*r+6)/(3*r+3)
(3*r + 6)/(3*r + 3)

[tscrim]

I don't know why that isn't reducing itself properly. That would be a bug IMO.

I don't think it's a bug: the gcd is one:

sage: a = 7*x; b = 7*y
sage: gcd(a, b)
1

Okay, but only because 7 is a unit. (Although if you do the gcd computation in QQbar, you get 7.) It is kind of a bug since it isn't detecting that the unit in the rational function is the leading coefficients of both...

[mantepse]

I don't know why that isn't reducing itself properly. That would be a bug IMO.

I don't think it's a bug: the gcd is one:

sage: a = 7*x; b = 7*y
sage: gcd(a, b)
1

Okay, but only because 7 is a unit. (Although if you do the gcd computation in QQbar, you get 7.) It is kind of a bug since it isn't detecting that the unit in the rational function is the leading coefficients of both...

Yes, but then

sage: R.<r,s> = QQ[]
sage: (3*r+6)/(3*r+3)
(3*r + 6)/(3*r + 3)

is also a bug. Put differently, it is a decision to take in multivariate polynomial normalization, I think.

[tscrim]

This is an interesting inconsistency (tested on 10.5.beta3):

sage: R.<x,y> = QQbar[]
sage: gcd(7*x, 7*y)
7
sage: (7*x).gcd(7*y)
7
sage: R.<x,y> = QQ[]
sage: gcd(7*x, 7*y)
1

Personally, I feel that putting at least the leading coefficient of the denominator when working over a field to 1 is desirable to normalize the elements. Well, I guess this is off topic to this PR at the end of the day.

[mantepse]

Personally, I feel that putting at least the leading coefficient of the denominator when working over a field to 1 is desirable to normalize the elements. Well, I guess this is off topic to this PR at the end of the day.

OK, I'll adapt the doctests then. I really dislike

            sage: P.<x,y> = ProjectiveSpace(Qp(3), 1)
            sage: f = DynamicalSystem_Berkovich([2*x^2, 2*y^2])
            sage: f.normalize_coordinates(); f
            Dynamical system of Projective Berkovich line over Cp(3) of precision 20
             induced by the map
              Defn: Defined on coordinates by sending (x : y)
               to ((2 + O(3^20))*x^2 : (2 + O(3^20))*y^2)

but I got no reaction on sage-devel.

[mantepse]

So, I just checked how univariate polynomial rings are treated specially. There is an implementation of the method fraction_field in PolynomialRing_field (which is for univariate polynomial rings only), which returns FractionField_1poly_field. On the other hand, there is apparently no special class for multivariate polynomial rings over fields.

[tscrim]

@bhutz Do you have any comments on the change to the dynamics doctest?

[bhutz]

I'm fine with those doc test changes for dynamics

[tscrim]

@bhutz Thanks for checking. It is uglier now since it isn't canceling the leading coefficients. Of course mathematically this is the same, but computationally this might make a difference...

[tscrim]

LGTM.

My only guess for why the coefficients are no longer so nice is that you immediately stop the content with a unit instead of 1. Perhaps we do want to add that into into the gcd since the coefficient rings want to do so for better normalization properties (if indeed that is the cause)?

[mantepse]

My only guess for why the coefficients are no longer so nice is that you immediately stop the content with a unit instead of 1.

I'm not sure. Indeed, one way the proposal differs from the old one is that the old version computes the content starting with the zero element, and stops when it reaches the one. But to make the differences precise seems non-obvious to me. For example, in the new version we have

sage: A.<a> = QQ[]
sage: B.<b> = A[]
sage: A._gcd_univariate_polynomial(B(b), B(7))
7
sage: A._gcd_univariate_polynomial(B(b), B(7*a))
1

In the old version, both give 7.

Perhaps we do want to add that into into the gcd since the coefficient rings want to do so for better normalization properties (if indeed that is the cause)?

I don't really know what you mean here. Could you be more explicit?

In preparing this comment, I discovered another potential optimisation: we currently have

            b = content(B)
            if b is not None:
                B //= b
...
            _, R = A.pseudo_quo_rem(B)
            while R.degree() > 0:
...
            if R.is_zero():
                b = content(B)

If we have R == 0 after computing pseudo_quo_rem, the while loop is not executed at all, and therefore B has not changed. In this case, it should not be necessary to compute the content of B again. However, at least in the examples from #38924 (comment) this does not seem to make any difference. Should I push it anyway?

diff --cc src/sage/categories/unique_factorization_domains.py
index d12ff868eb3,d12ff868eb3..f544bf5d098
--- a/src/sage/categories/unique_factorization_domains.py
+++ b/src/sage/categories/unique_factorization_domains.py
@@@ -210,9 -210,9 +210,14 @@@ class UniqueFactorizationDomains(Catego
                  d = a.gcd(b)
              else:
                  d = one
--            g = h = one
--            delta = A.degree() - B.degree()
++
              _, R = A.pseudo_quo_rem(B)
++            if R.is_zero():
++                return d * B
++            if not R.degree():
++                return f.parent()(d)
++            delta = A.degree() - B.degree()
++            g = h = one
              while R.degree() > 0:
                  A = B
                  h_delta = h**delta

[mantepse]

I tried to change is_unit to is_one. This does not change the result of the doctest in arith.misc.py.

I tried to find a ring where is_unit takes significantly longer than is_one, and where UniqueFactorizationDomains.ParentMethods._gcd_univariate_polynomial is called. However, so far I failed. A natural candidate would be the various quotient rings, but they don't come with a gcd method.

[tscrim]

Perhaps working over Zmod(n) where n is composite with large-ish prime factors?

Hmm..interesting, I am somewhat surprised by this. I thought that would have been the difference. Well, in that case, we might want to leave it as is_unit()…

[mantepse]

Perhaps working over Zmod(n) where n is composite with large-ish prime factors?

I don't think that I can create a UFD from that.

[tscrim]

Ah, right, we need a UFD. Then perhaps $Z[x]$ mod a reasonably high degree irreducible polynomial with large coefficients?

[mantepse]

I don't know how to do that:

sage: R.<x> = ZZ[]
sage: S = R.quotient(x^2 + 1)
sage: S in UniqueFactorizationDomains()
False
sage: S._refine_category_(UniqueFactorizationDomains())
sage: S.an_element().gcd(S.an_element())
...
AttributeError: 'PolynomialQuotientRing_domain_with_category.element_class' object has no attribute 'gcd'

See also #38937

[tscrim]

Hmmm...I don't know what a good example, much less one that is also easy to construct in Sage. Perhaps we should ask for an example on sage-devel.

[mantepse]

Since nobody was able to come up with an example, could we please deal with that problem once it arises?