Implement Square Roots to FiniteFields Category

[Brian-Heckel]

Fixes #39688

Added the following methods to the finite field category

• quadratic_non_residue()
• square_root() and sqrt()
• _tonelli() and _cipolla(): the algorithms for finding a square root
• is_square() Moved is_square() and sqrt() methods from element.pyx to their proper categories.

Most profiling tests showed that Tonelli's algorithm outperforms Cipolla. However a profile with a Finite Field with order of a Cullen prime resulted in Cipolla outperforming Tonelli.

I'm running the meson build and currently cannot find the way to build the html documentation.

📝 Checklist

• The title is concise and informative.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.
• I have updated the documentation and checked the documentation preview.

⌛ Dependencies

[vincentmacri]

Brian did this work as part of an undergrad mentorship program at the University of Calgary, and I was the grad student mentoring him. I'll review this for mathematical correctness/style/documentation/etc., but since I worked closely with him on this I'd like a second person to review it as well (especially for the refactoring that was done to move the default square root implementation out of element.pyx and into the appropriate category).

[vincentmacri]

First pass, mostly style things and wording in the documentation.

[vincentmacri]

Docstring fixes

[vincentmacri]

Docstring formatting

[user202729]

Not completely related to this change:

@tobiasdiez another ruff failure here, even though the code is just moved from one file to another without any whitespace change. What's going on? (different rules are applied on .py and .pyx files, I guess? Is there a linter for pyx?)

[user202729]

do fix the lints.

Also, build the HTML files locally (the CI for HTML preview is currently broken) and check that they're rendered correctly. They're incorrect in many places.

Check some existing docstring for the formatting. I think first word should be in infinitive form, first sentence has full stop, and no bullet point in output if only one output. (Do verify that.)

Refer to #39798 .

First character of error message should not be capitalized.

What's the algorithm for computing square root in e.g. GF(5^2)? If it's also tonelli or cipolla, consider giving a different name for the sage implementation e.g. sage_tonelli to distinguish it.

[tobiasdiez]

Not completely related to this change:

@tobiasdiez another ruff failure here, even though the code is just moved from one file to another without any whitespace change. What's going on? (different rules are applied on .py and .pyx files, I guess? Is there a linter for pyx?)

Yes, ruff is not supporting cython.

[vincentmacri]

Docstrings, cleanup.

[vincentmacri]

do fix the lints.

Curious, did you approve the CI run? I don't have permission to approve CI runs was wondering who did (so I can ping them to do it again).

Also, build the HTML files locally (the CI for HTML preview is currently broken) and check that they're rendered correctly. They're incorrect in many places.

That's annoying that the CI is broken. The process to build the docs with Meson is not very well-documented (I should get around to filing a bug report about that eventually...), but I eventually figured it out and sent the instructions to Brian.

What's the algorithm for computing square root in e.g. GF(5^2)? If it's also tonelli or cipolla, consider giving a different name for the sage implementation e.g. sage_tonelli to distinguish it.

I'm not too sure what you mean by wanting to distinguish with the naming, but I don't think we need to worry about it. GF uses optimized implementations written in Cython, so the sqrt added here will not be called by an element of GF(q), only by elements of quotient rings that are known to be fields by Sage (or any other way one manages to construct something belonging to the FiniteFields category that isn't an instance of GF, but the only such way I know of is the polynomial quotient ring). Notably, the GF implementations do not have an algorithm option for their sqrt methods.

In terms of what algorithm GF uses, for GF(5^2) we use the Givaro library. For GF(p) we use Tonelli (might be a variant of Tonelli, not 100% sure). For GF(p^n) where p^n is large and n > 1 we use Pari. There might be other situations in which we do different things or use other libraries.

That said, I do think that at least for GF(p) when p is large and p - 1 is highly divisible by 2, the _cipolla function can outperform the Cython integer_mod Tonelli implementation. One such case we found was the Cullen prime p = 141 * 2^141 + 1. For future work (i.e. I consider out of scope for this PR) we might want to do some more detailed profiling on exactly when Cipolla outperforms Tonelli and either implement that in Cython and use it when p is large enough, or just call the Python _tonelli implementation added by this PR.

[user202729]

I'm not too sure what you mean by wanting to distinguish with the naming, but I don't think we need to worry about it. GF uses optimized implementations written in Cython, so the sqrt added here will not be called by an element of GF(q), only by elements of quotient rings that are known to be fields by Sage (or any other way one manages to construct something belonging to the FiniteFields category that isn't an instance of GF, but the only such way I know of is the polynomial quotient ring). Notably, the GF implementations do not have an algorithm option for their sqrt methods.

Yes, I know. But that is a breakage of Liskov substitution principle — in general, elements of parents in FiniteFields() category supports .sqrt(algorithm=), except elements of GF(p^e), even though GF(p^e) in FiniteFields().

In practice, yes, the user probably wouldn't care and just not pass algorithm= either way, I agree. But changing the external interface requires deprecation period, I'd prefer being extra careful.

or just call the Python _tonelli implementation added by this PR.

Bad idea. Methods with name start with underscore should not be called by the user.

[vincentmacri]

Bad idea. Methods with name start with underscore should not be called by the user.

I meant Sage itself could call that, not the user.

[vincentmacri]

Yes, I know. But that is a breakage of Liskov substitution principle — in general, elements of parents in FiniteFields() category supports .sqrt(algorithm=), except elements of GF(p^e), even though GF(p^e) in FiniteFields().

In practice, yes, the user probably wouldn't care and just not pass algorithm= either way, I agree. But changing the external interface requires deprecation period, I'd prefer being extra careful.

Sage's current behaviour doesn't respect this principle. The implementation of is_square in CommutativeRingElement has an optional root parameter. The implementation of is_square for elements of GF(p^e) doesn't have an optional root parameter. And yet:

K = GF(5)
a = K(3)
from sage.structure.element import CommutativeRingElement
isinstance(a, CommutativeRingElement)  # True

[user202729]

I consider that a bug. While the behavior remains, we could at least not introduce more of it.

[Brian-Heckel]

So the options to deal with this are

1. Add the algorithm parameter to the sqrt methods in GF implementations
2. Remove algorithm parameter to sqrt that I'm adding and provide a heuristic for which algorithm to use (or add sqrt_cipolla)
3. Leave it as is. Maybe I can add a warning to the documentation that GF implementations don't have algorithm so that users are aware.

[vincentmacri]

I consider that a bug. While the behavior remains, we could at least not introduce more of it.

Fair enough. For this PR I think not breaking the Liskov substitution principle with respect to the GF implementations is all that's realistic to aim for. There's too much inconsistent behaviour for square roots between the rings code and finite fields code to address in this PR.

So the options to deal with this are

1. Add the `algorithm` parameter to the `sqrt` methods in `GF` implementations

2. Remove `algorithm` parameter to `sqrt` that I'm adding and provide a heuristic for which algorithm to use (or add `sqrt_cipolla`)

3. Leave it as is. Maybe I can add a warning to the documentation that GF implementations don't have `algorithm` so that users are aware.

Using a heuristic to decide what algorithm to use is only a good option if we have an accurate heuristic for which algorithm is better for all cases, and I don't think we do. If there are any cases where Cipolla is meaningfully faster than Tonelli but the heuristic still decides to use Tonelli then we'd still want to have an optional algorithm parameter.

I don't like the idea of adding a public sqrt_cipolla method. As far as I know every other function in Sage that supports multiple algorithms has an optional algorithm parameter, and not public methods with different names.

In theory my preference would be to add the algorithm parameter to sqrt in the GF implementations and allow the user to choose Cipolla. Then I looked at the GF implementations of sqrt and there are a lot of them all with slightly different behaviour, especially with regards to extend. This would also cause a lot of code duplication. I think for now perhaps the best option if @user202729 doesn't want to merge this with the Liskov substitution principle violation is to add an optional algorithm parameter that defaults to None to all the sqrt implementations for the various finite fields in src/sage/rings/finite_rings and don't do anything with the parameter, just ignore it (as many of those implementations do for extend currently). Perhaps the biggest issue is that Zmod(6) and GF(7) use the same sqrt function, and Cipolla relies on the fact that we are working over a field.

Given those issues, I think the only realistic solution is to add the algorithm parameter to GF but don't use it in the implementations. Leave refactoring the sqrt code for finite fields to all have the same interface and and ability to use Cipolla for a future PR.

@user202729 What's your preference for fixing this? Any other ideas?

[user202729]

I was thinking of something like

class FiniteFieldGeneric:
    _default_sqrt_algorithm = "sage_tonelli"

    def sqrt(algorithm=None):
        if algorithm is None:
            algorithm = self._default_sqrt_algorithm
        if algorithm == "sage_tonelli":
            return self._sqrt_sage_tonelli()
        elif algorithm == "pari":
            if hasattr(self, '_sqrt_pari'):  # alternatively can also provide a default implementation of `_sqrt_pari` that raise NotImplementedError or ValueError
                return self._sqrt_pari()  # else fallthrough
        elif ...:
            ...
        raise NotImplementedError(f"algorithm {algorithm} not supported for elements of {self.parent()}")

    def _sqrt_sage_tonelli():
        ...


class FiniteFieldPari(FiniteFieldGeneric):
    def _sqrt_pari(self):
        ...

    _default_sqrt_algorithm = "pari"

class FiniteFieldGivaro(FiniteFieldGeneric):
    def _sqrt_givaro(self):
        ...

    _default_sqrt_algorithm = "givaro"

but then there's extend= which I didn't account for.

While it ought to be possible to deduplicate/clean up things, depends on how difficult it is to do so, it may be reasonable to just get the current pull request in in the current situation anyway.

In that case, my only remaining complaint is the naming of the algorithm. Assume, in the future we add .sqrt(algorithm="tonelli") and .sqrt(algorithm="givaro") to GF(p^n) (let's say it uses givaro), then

• it would be confusing if givaro also uses tonelli internally,
• it isn't at all clear that .sqrt(algorithm="tonelli") is slower than .sqrt(algorithm="givaro") (which is why I suggest naming it sage_tonelli, then it's clear that it's an implementation in sage)

[vincentmacri]

While it ought to be possible to deduplicate/clean up things, depends on how difficult it is to do so, it may be reasonable to just get the current pull request in in the current situation anyway.

I think it is definitely possible but it would be large enough to warrant its own PR. I would be in favour of getting this PR in in its current state and opening a separate issue to clean up the finite field square root code so its less likely to be forgotten about.

So I think we're in agreement that this is ready to merge, pending the CI run. I do disagree with the suggestion to call the algorithm sage_tonelli though as I'll explain below.

In that case, my only remaining complaint is the naming of the algorithm. Assume, in the future we add .sqrt(algorithm="tonelli") and .sqrt(algorithm="givaro") to GF(p^n) (let's say it uses givaro), then it would be confusing if givaro also uses tonelli internally,

Sage doesn't make any guarantees on what algorithm external libraries use internally, so I don't think it would be confusing if Givaro uses Tonelli internally (it doesn't) because nowhere would Sage document that fact or guarantee it to be true if your Givaro version changes.

Additionally, I don't see a reason to make the Tonelli implementation added in this PR available for the Givaro/NTL/Flint implementations of finite fields, whatever they have is likely faster as it's written in a lower-level language like C. So that's another reason why there shouldn't be risk of confusion. Indeed, this PR does not allow one to use either the Tonelli or Cipolla implementations with any of the existing finite field implementations (at least not through public methods). The only way to run either of the square roots algorithms implemented in this PR is irreducible polynomial quotient rings over a finite field (or call the private _tonelli or _cipolla methods directly, but that's not supported and users shouldn't do that and expect things not to break between versions).

it isn't at all clear that .sqrt(algorithm="tonelli") is slower than .sqrt(algorithm="givaro")

For what it's worth, Givaro represents elements as g^k where g is a generator of the multiplicative group. So Givaro takes square roots by checking if k is even and if so the square root is g^(k // 2). This algorithm is very fast if you're representing elements in this form, and doesn't make sense otherwise. I don't see a reason for Givaro to support other square root algorithms than this. It also doesn't make sense to use algorithm=givaro unless you're using the Givaro finite field implementation.

If Givaro is slower then that would indicate to me that it needs to be optimized (our Givaro implementation does some arithmetic mod 2 and division by 2 which could probably be replaced by bit operations). If you've measured them and don't see a substantial difference in speed then that's probably explained by Givaro only allowing finite fields of size up to 2^16. I checked and there aren't any prime powers q where q - 1 is highly divisible by two (which is the case where Tonelli is slow) for any 0 <= q <= 2^16.

In fact, Tonelli (very broadly speaking, skipping a lot of details here) computes sqrt(a) by finding an element g and an even integer k such that a = g^k and then returns sqrt(a) = g^(k // 2). So if you already have a = g^k like we do with Givaro you wouldn't use Tonelli (or Cipolla).

* it isn't at all clear that `.sqrt(algorithm="tonelli")` is slower than `.sqrt(algorithm="givaro")` (which is why I suggest naming it `sage_tonelli`, then it's clear that it's an implementation in sage)

I think prefixing the algorithm name with sage_ would be inconsistent with how we name algorithm options in Sage. We usually name the algorithm after either what the algorithm is called by mathematicians (so 'tonelli' or 'cipolla') or in the case where we are calling an external library we name it by the name of the external library ('giavaro', 'pari', etc.). We do have a few instances where 'sage' is an algorithm option, but I don't think we ever have 'sage_[name of an existing algorithm]'.

While obviously this is subjective, I'd actually find it more confusing if it were called sage_tonelli. If I were reading the documentation and saw an algorithm option called sage_tonelli I would assume it's a variant of the Tonelli algorithm invented by a Sage contributor (whereas @Brian-Heckel's implementation is pretty much a line-by-line implementation of the algorithm in the textbook he cited).

Basically, Tonelli(-Shanks) seems to be the fastest algorithm for most cases, although one can certainly construct finite fields where Cipolla is faster (as @Brian-Heckel shows in one of the examples in the sqrt docstring). So the only sensible options for a finite field square root are whatever the default square root is for that implementation, and possibly also Cipolla. The Tonelli implementation in this PR should only be used by polynomial quotient rings where the quotient is an irreducible, as every other way to construct a finite field in Sage has a more tailored implementation.

In terms of your code snippet, this is how I think it should work (after refactoring, not for this PR). This is just in terms of the logic, the categories framework and the fact that the GF implementations don't get most of their behaviours from it complicates what the real code would look like of course.

class FiniteFieldElementGeneric:

    def sqrt(algorithm=None):
        if algorithm == 'cipolla':
            return self._cipolla()
        return self._sqrt()  # By defaulting to a sensible algorithm rather than throwing an error, it's easier for us to add/remove square root algorithms without breaking changes.

    def _tonelli(self):
        ...

    _sqrt = _tonelli

    def _cipola(self):
        ...


class FiniteFieldElementPari(FiniteFieldElementGeneric):
    def _sqrt_pari(self):
        ...

    # I haven't  checked what PARI does, but for this example let's assume it uses Tonelli, in which case there is no reason to offer the Sage Tonelli implementation when we have a C implementation available
    _sqrt = _sqrt_pari

class FiniteFieldElementGivaro(FiniteFieldElementGeneric):
    def _sqrt_givaro(self):
        ...

    def sqrt(self, algorithm=None):
        return self._sqrt_givaro()  # The Givaro implementation is optimized enough that it doesn't make sense to use other algorithms, although it may be hard to see that in timings since Givaro only works for small finite fields and all algorithms are fast for small finite fields.

To make this a lot shorter: I don't think every finite field implementation needs to support every square root algorithm, only the ones that someone would reasonably want to use. For the default category implementation this is Tonelli and Cipolla. For other implementations it would be whatever they currently use and maybe also Cipolla (with the exception of Givaro which I think should only use the algorithm it has now).

If you still find the tonelli algorithm name confusing, for the sake of consistency with all the other algorithm options in Sage, I would suggest that rather than add the sage_ prefix to the algorithm name that we simply explain in the docstring that sqrt(algorithm='tonelli') uses an implementation of Tonelli written in Sage (and if so, we might as well for Cipolla as well).