Allow "generic" PolynomialRing implementation

[jdemeyer]

It is possible to specify a specific backend like PolynomialRing(QQ, 'x', implementation="singular") but the opposite (asking for a generic Sage implementation) is not. Fix this by allowing implementation="generic".

This requires various fixes involving the implementation keyword being ignored or breaking caching. In particular, we fix this caching bug:

sage: p = 2^64 + 13
sage: A = GF(p^2)
sage: B = GF(p^3)
sage: R = A.modulus().parent()
sage: S = B.modulus().parent()
sage: R is S
False

We also introduce a new implementation "GF2X" for the GF2X package. This used be called "NTL" which we still allow as alias.

CC: @embray @tscrim

Component: basic arithmetic

Author: Jeroen Demeyer

Branch/Commit: 0efda1b

Reviewer: Travis Scrimshaw

Issue created by migration from https://trac.sagemath.org/ticket/24264

[embray]

comment:1

+1

[jdemeyer]

Description changed:

--- 
+++ 
@@ -1 +1,3 @@
 It is possible to specify a specific backend like `PolynomialRing(QQ, 'x', implementation="singular")` but the opposite (asking for a generic Sage implementation) is not. Fix this by allowing `implementation="generic"`.
+
+This requires minor fixes in `src/sage/rings/polynomial/polynomial_ring_constructor.py`

[jdemeyer]

Author: Jeroen Demeyer

[jdemeyer]

comment:5

This is harder than expected due to caching. Currently, the caching does not take into account the implementation keyword. That is an existing bug, but it makes it pointless to add more support for the implementation keyword.

[embray]

comment:6

Wouldn't it make sense to fix that? It doesn't need to use the "implementation" keyword per se, but it would make sense if the type a method was called on was also cached?

I'm likely missing something though since I don't know exactly what you mean in this case by "caching".

[jdemeyer]

Description changed:

--- 
+++ 
@@ -1,3 +1,13 @@
 It is possible to specify a specific backend like `PolynomialRing(QQ, 'x', implementation="singular")` but the opposite (asking for a generic Sage implementation) is not. Fix this by allowing `implementation="generic"`.

-This requires minor fixes in `src/sage/rings/polynomial/polynomial_ring_constructor.py`
+This requires various fixes involving the `implementation` and caching. In particular, we fix this caching bug:
+
+```
+sage: p = 2^64 + 13
+sage: A = GF(p^2)
+sage: B = GF(p^3)
+sage: R = A.modulus().parent()
+sage: S = B.modulus().parent()
+sage: R is S
+False
+```

[jdemeyer]

comment:8

Replying to @embray:

Wouldn't it make sense to fix that?

Of course. I'm just saying that this makes it harder than I initially expected.

[jdemeyer]

Branch: u/jdemeyer/allowgenericpolynomialring_implementation

[jdemeyer]

Commit: f8449e9

[jdemeyer]

New commits:

f8449e9

Allow "generic" PolynomialRing implementation

[jdemeyer]

Description changed:

--- 
+++ 
@@ -1,6 +1,6 @@
 It is possible to specify a specific backend like `PolynomialRing(QQ, 'x', implementation="singular")` but the opposite (asking for a generic Sage implementation) is not. Fix this by allowing `implementation="generic"`.

-This requires various fixes involving the `implementation` and caching. In particular, we fix this caching bug:
+This requires various fixes involving the `implementation` keyword being ignored or breaking caching. In particular, we fix this caching bug:

sage: p = 2^64 + 13

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

b8990cb

Allow "generic" PolynomialRing implementation

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from f8449e9 to b8990cb

[tscrim]

comment:15

This change:

diff --git a/src/sage/algebras/free_algebra.py b/src/sage/algebras/free_algebra.py
index e6c66e9..09db35d 100644
--- a/src/sage/algebras/free_algebra.py
+++ b/src/sage/algebras/free_algebra.py
@@ -269,7 +269,7 @@ class FreeAlgebraFactory(UniqueFactory):
         # test if we can use libSingular/letterplace
         if implementation == "letterplace":
             args = [arg for arg in (arg1, arg2) if arg is not None]
-            kwds = dict(sparse=sparse, order=order, implementation="singular")
+            kwds = dict(order=order, implementation="singular")
             if name is not None:
                 kwds["name"] = name
             if names is not None:

seems to cause multiple failures with the free algebra:

sage -t --long src/sage/algebras/free_algebra.py  # 12 doctests failed
sage -t --long src/sage/algebras/letterplace/letterplace_ideal.pyx  # Bad exit: 14
sage -t --long src/sage/algebras/letterplace/free_algebra_letterplace.pyx  # 19 doctests failed
sage -t --long src/sage/algebras/letterplace/free_algebra_element_letterplace.pyx  # 31 doctests failed

See patchbots.

Also, we (and singular?) only have sparse implementations for multivariate polynomials, right?

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

a299cb4

Allow "generic" PolynomialRing implementation

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from b8990cb to a299cb4

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

21b4952

Fix -Wsign-compare compiler warning

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from a299cb4 to 21b4952

[jdemeyer]

comment:20

Replying to @tscrim:

seems to cause multiple failures with the free algebra:

Fixed by allowing sparse=None to mean sparse=True in the multi-variate case.

[jdemeyer]

comment:21

Patchbot is green now.

[tscrim]

comment:22

What are your thoughts about changing _implementation_names to _implementation, which only takes a string? It is only used to differentiate between FLINT and NTL implementations in coercion. I think it could simplify the code overall.

Furthermore, how do you feel about instead of storing a key for all implementation keywords possible, we instead standardize the implementation keyword to always give the actual implementation. I know it is a little more brittle, but I think it will make the code more clean and would work well with the above suggestion. I haven't yet fully attempted this, but I will be happy to do this if you do not object.

One thing that definitely needs fixing: By allowing generic coercions, we are introducing a memory leak via cyclic coercion between the generic and FLINT implementations:

sage: RF.<x> = PolynomialRing(ZZ, implementation='FLINT')
sage: RG.<x> = PolynomialRing(ZZ, implementation='generic')
sage: RF.coerce_map_from(RG)
sage: RF.has_coerce_map_from(RG)
True
sage: RG.has_coerce_map_from(RF)
True

I've been told that this circular coercion gives a memory leak because the co(?)domain is cached by the coercion model. Although I cannot reproduce it:

sage: import gc
sage: for i in range(10000):
....:     if i % 100 == 0:
....:         get_memory_usage()
....:     R = PolynomialRing(ZZ, 'x%s'%i, implementation='FLINT')
....:     S = PolynomialRing(ZZ, 'x%s'%i, implementation='generic')
....:     assert R is not S
....:     assert R.coerce_map_from(S) is not None
....:     assert S.coerce_map_from(R) is not None
....:     del R
....:     del S
....:     _ = gc.collect()
5872.29296875
5872.29296875
5872.29296875
5872.29296875
...

However, it does introduce a very subtle reason why code could suddenly become slow by multiplying things in the opposite order. We will have to update things accordingly in that spot in PolynomialRing._coerce_map_from_.

[jdemeyer]

comment:23

Replying to @tscrim:

Furthermore, how do you feel about instead of storing a key for all implementation keywords possible, we instead standardize the implementation keyword to always give the actual implementation.

Well, this will slow down creating polynomial rings. Every time that you create a polynomial ring with implementation=None (which is the default, so the most common case), you will need to run some code to figure out the implementation.

Let me give a counter-proposal: instead of making _implementation_names an attribute, we create a static method on the polynomial ring class:

# Assuming that the implementation for a given polynomial ring class
# depends only on the base_ring, which is currently the case.
@staticmethod
def _implementation_names(base_ring, implementation):
    # Default implementation for generic polynomial rings
    if implementation is None or implementation == "generic":
        return ["generic", None]
    raise ValueError("unknown implementation {!r}".format(implementation))

I think that this will lead to cleaner code without the performance loss of normalizing the implementation every time.

[tscrim]

comment:24

Replying to @jdemeyer:

Replying to @tscrim:

Furthermore, how do you feel about instead of storing a key for all implementation keywords possible, we instead standardize the implementation keyword to always give the actual implementation.

Well, this will slow down creating polynomial rings. Every time that you create a polynomial ring with implementation=None (which is the default, so the most common case), you will need to run some code to figure out the implementation.

Ah, that is a good point about having to process the input.

Let me give a counter-proposal: instead of making _implementation_names an attribute, we create a static method on the polynomial ring class:

# Assuming that the implementation for a given polynomial ring class
# depends only on the base_ring, which is currently the case.
@staticmethod
def _implementation_names(base_ring, implementation):
    # Default implementation for generic polynomial rings
    if implementation is None or implementation == "generic":
        return ["generic", None]
    raise ValueError("unknown implementation {!r}".format(implementation))

I think that this will lead to cleaner code without the performance loss of normalizing the implementation every time.

+1 Do you want me to do this or will you?

[jdemeyer]

comment:25

Let me give this a try.

[jdemeyer]

Description changed:

--- 
+++ 
@@ -11,3 +11,5 @@
 sage: R is S
 False

+ +We also introduce a new implementation "GF2X" for the GF2X package. This used be called "NTL" which we still allow as alias.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

d679f13	Fix -Wsign-compare compiler warning
bfd388f	Allow "generic" PolynomialRing implementation

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 21b4952 to bfd388f

[jdemeyer]

comment:28

Replying to @jdemeyer:

Let me give this a try.

Done. This was more work than I initially thought, but I am happy with the result. The code looks cleaner than before.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

a0aefd8

Allow "generic" PolynomialRing implementation

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from bfd388f to a0aefd8

[tscrim]

comment:30

It definitely is cleaner to me as well. Some comments:

• This seems fragile:

      implementation = self._implementation_names(implementation, base_ring)[0]
      if implementation == "NTL":

  I would instead do

      implementation = self._implementation_names(implementation, base_ring)[0]
      if "NTL" in implementation:

• The cyclic coercion (comment:22) should get a doctest for only the generic to FLINT (if I am reading the code correctly).

• This feels a bit fragile to me

  if self.element_class.__name__ != 'Polynomial_integer_dense_flint':

  I'm guessing you are doing this to avoid a circular import? What about doing something similar to what is currently there and checking the implementation name?

• Some _implementation_names_impl have doctests and others do not (including some with mildly non-trivial behavior).

Can you also explain why you have the @classmethod and the @staticmethod? I'm not quite sure I understand why we need the redirection.

[jdemeyer]

comment:31

Replying to @tscrim:

• This seems fragile:

      implementation = self._implementation_names(implementation, base_ring)[0]
      if implementation == "NTL":

  I would instead do

      implementation = self._implementation_names(implementation, base_ring)[0]
      if "NTL" in implementation:

Why? I don't understand your objection.

• Some _implementation_names_impl have doctests and others do not (including some with mildly non-trivial behavior).

Well, all of them are tested somewhere. Some are tested explicitly by calling _implementation_names_impl; some are tested indirectly by constructing an instance of a polynomial ring class; still others are testing using the polynomial ring constructor.

Can you also explain why you have the @classmethod and the @staticmethod? I'm not quite sure I understand why we need the redirection.

It was just simpler that way. First of all, the actual implementations can just return NotImplemented instead of raising an exception. This also simplifies one check in the polynomial constructor where I'm using list comprehension to check for NotImplemented. You cannot use list comprehension with a condition "raises ValueError".

Second, the wrapper method _implementation_names() does a sanity check of the output of _implementation_names_impl(): it checks that the given implementation is in the output list.

[jdemeyer]

comment:32

Replying to @tscrim:

One thing that definitely needs fixing: By allowing generic coercions, we are introducing a memory leak via cyclic coercion between the generic and FLINT implementations:

That is a very deep rabbit hole.

First of all, note that cyclic coercion is already possible in other cases:

sage: R = PolynomialRing(GF(3), 't', implementation='FLINT')
sage: S = PolynomialRing(GF(3), 't', implementation='NTL')
sage: R
Univariate Polynomial Ring in t over Finite Field of size 3
sage: S
Univariate Polynomial Ring in t over Finite Field of size 3 (using NTL)
sage: R.has_coerce_map_from(S)
True
sage: S.has_coerce_map_from(R)
True

There are two related issues when trying to fix this:

(A) The polynomial ring itself does not store its implementation or any other information which could help to determine in which way a coercion should go.

(B) Pickling doesn't work correctly for non-default implementations:

# Continuing the example above
sage: loads(dumps(S)) is S
False
sage: loads(dumps(S)) is R
True

Currently, one could fix (A) by not allowing coercions between different implementations. But then (B) gives trouble because then there isn't even a coercion between loads(dumps(S)) and S.

Suggestions??

[tscrim]

comment:33

Replying to @jdemeyer:

Replying to @tscrim:

• This seems fragile:

      implementation = self._implementation_names(implementation, base_ring)[0]
      if implementation == "NTL":

  I would instead do

      implementation = self._implementation_names(implementation, base_ring)[0]
      if "NTL" in implementation:

Why? I don't understand your objection.

I either missed, didn't remember, or didn't quite understand the implication of the "The first element in the list is the canonical name." part of the OUTPUT. So I thought you were using an undocumented specification.

• Some _implementation_names_impl have doctests and others do not (including some with mildly non-trivial behavior).

Well, all of them are tested somewhere. Some are tested explicitly by calling _implementation_names_impl; some are tested indirectly by constructing an instance of a polynomial ring class; still others are testing using the polynomial ring constructor.

Would it be possible to have at least one of those tests local to each of these functions? It would make it easier to debug if something broke later on.

Can you also explain why you have the @classmethod and the @staticmethod? I'm not quite sure I understand why we need the redirection.

It was just simpler that way. First of all, the actual implementations can just return NotImplemented instead of raising an exception. This also simplifies one check in the polynomial constructor where I'm using list comprehension to check for NotImplemented. You cannot use list comprehension with a condition "raises ValueError".

Second, the wrapper method _implementation_names() does a sanity check of the output of _implementation_names_impl(): it checks that the given implementation is in the output list.

Thank you for the explanations. Makes sense to me, and it made it easy to see where these features are being used.

[tscrim]

comment:34

Replying to @jdemeyer:

Replying to @tscrim:

One thing that definitely needs fixing: By allowing generic coercions, we are introducing a memory leak via cyclic coercion between the generic and FLINT implementations:

That is a very deep rabbit hole.

First of all, note that cyclic coercion is already possible in other cases:

sage: R = PolynomialRing(GF(3), 't', implementation='FLINT')
sage: S = PolynomialRing(GF(3), 't', implementation='NTL')
sage: R
Univariate Polynomial Ring in t over Finite Field of size 3
sage: S
Univariate Polynomial Ring in t over Finite Field of size 3 (using NTL)
sage: R.has_coerce_map_from(S)
True
sage: S.has_coerce_map_from(R)
True

I know some people would consider that a bug, and even if is not, it could create strange behavior and slowdowns because of moving between implementations.

There are two related issues when trying to fix this:

(A) The polynomial ring itself does not store its implementation or any other information which could help to determine in which way a coercion should go.

It has to store something in order to get the different string representations. As for determining which way to make the coercion is a bit more arbitrary, but something we can do in the _coerce_map_from_.

(B) Pickling doesn't work correctly for non-default implementations:

# Continuing the example above
sage: loads(dumps(S)) is S
False
sage: loads(dumps(S)) is R
True

Currently, one could fix (A) by not allowing coercions between different implementations. But then (B) gives trouble because then there isn't even a coercion between loads(dumps(S)) and S.

That's a bug (and somewhat horrifying).

Suggestions??

(B) deserves a ticket all on its own. Although (A) is currently addressed for ZZ rings:

sage: R = PolynomialRing(ZZ, 't', implementation="FLINT")
sage: S = PolynomialRing(ZZ, 't', implementation="NTL")
sage: R.has_coerce_map_from(S)
True
sage: S.has_coerce_map_from(R)
False

I can't think of a simple way with the current framework that we can have an ordering on the implementations that would correspond to the coercion order. I feel like we would need to add another @staticmethod specifying this (with a default of [None] when there is only 1 implementation, the default).

I am also happy addressing (A) more thoroughly on a followup if we just add a test for

sage: T = PolynomialRing(ZZ, 't', implementation="generic")

and the non-cyclic coercion with R and S above.

[jdemeyer]

comment:35

Replying to @tscrim:

I am also happy addressing (A) more thoroughly on a followup if we just add a test for

sage: T = PolynomialRing(ZZ, 't', implementation="generic")

and the non-cyclic coercion with R and S above.

Yes, I would rather do that. This ticket is already getting quite big...

[jdemeyer]

comment:36

Follow-up tickets: #24319, #24320 (Note that I do not plan to work on these)

[jdemeyer]

comment:37

Part of the mess in polynomial rings is that the logic for constructing them is separated in different places:

1. The constructor PolynomialRing

2. The _implementation_names_impl() method of the polynomial ring.

3. The __init__ method of the polynomial ring.

There are also two non-equivalent ways to pass an implementation: using the implementation keyword or using the element_class keyword.

[jdemeyer]

comment:38

I feel like only PolynomialRing should take an implementation and not __init__. Instead, the implementation should be translated somewhere to an element_class which is passed to __init__.

But that will be for a follow-up.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

0020291	Fix -Wsign-compare compiler warning
ee237be	Allow "generic" PolynomialRing implementation

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from a0aefd8 to ee237be

[jdemeyer]

comment:40

I hope that this addresses all your concerns.

[tscrim]

Reviewer: Travis Scrimshaw

[tscrim]

comment:41

Yes, it does. Thank you.

[jdemeyer]

comment:42

Finally... I initially started with this ticket thinking "this should only take 1 hour" :-)

[tscrim]

comment:43

Ah, I know that feeling. Although in my case, that has always been a sign I am in for a whole day of programming. :P

[embray]

comment:44

Replying to @jdemeyer:

Finally... I initially started with this ticket thinking "this should only take 1 hour" :-)

I "love" when that happens. I haven't really followed this closely since it quickly went beyond any understanding I have of the coercion system, but I'm glad it will be fixed.

The next thing to do might be to fix cases like the one I pointed out in #24263 to ensure that the polynomial implementation that's supposed to be being tested is actually the one being tested.

[vbraun]

comment:45

See patchbot

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from ee237be to fff5d7e

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

f33025e	Fix coercion of polynomial rings
fff5d7e	Better error message for mismatched constructions