Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Total refactor of any_root() to solve issue in characteristic two and clean up the code #37170

Merged
merged 31 commits into from
Feb 13, 2024

Conversation

GiacomoPope
Copy link
Contributor

@GiacomoPope GiacomoPope commented Jan 26, 2024

This PR was inspired by the issue #37160, which found that computing isomorphisms between fields of even characteristic was taking much much longer than expected (20s vs 100ms).

My first attempt (see commit history) were to simply patch up the any_root() function, but this yielded some results which left the code in an even more confusing state.

Following the advice of the reviewers, I have now implemented any_irreducible_factor() which given a polynomial element computes an irreducible factor of this polynomial with minimal degree. When the optional parameter degree=None is set, then an irreducible factor of degree degree is computed and a ValueError is raise if no such factor exists.

Now, any_root() becomes simply a call to self. any_irreducible_factor(degree=1) and a root is found from this linear polynomial. We also handle all the other cases of optional arguments handled by the old function, so the function should behave as before, but with a cleaner code to read.

Before PR

┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 10.2, Release Date: 2023-12-03                    │
│ Using Python 3.11.4. Type "help()" for help.                       │
└────────────────────────────────────────────────────────────────────┘
sage: set_random_seed(0)
....: K.<z8> = GF(2^8)
....: R.<x> = K[]
....: u = R.irreducible_element(2)
....: K_ext = K.extension(modulus=u, names="a")
....: R_ext.<y_ext> = K_ext[]
....: poly = R_ext.random_element(2).monic()
....: %time poly.roots()
CPU times: user 20.5 s, sys: 19.9 ms, total: 20.5 s
Wall time: 20.5 s
[]

After PR

┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 10.3.beta6, Release Date: 2024-01-21              │
│ Using Python 3.11.4. Type "help()" for help.                       │
└────────────────────────────────────────────────────────────────────┘
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Warning: this is a prerelease version, and it may be unstable.     ┃
┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛
sage: set_random_seed(0)
....: K.<z8> = GF(2^8)
....: R.<x> = K[]
....: u = R.irreducible_element(2)
....: K_ext = K.extension(modulus=u, names="a")
....: R_ext.<y_ext> = K_ext[]
....: %time R_ext.random_element(2).roots()
CPU times: user 110 ms, sys: 9.03 ms, total: 119 ms
Wall time: 150 ms
[]

This fixes #37160 but i think more feedback and thorough doctests need to be included considering this is a very old function which many projects will be using.

@GiacomoPope
Copy link
Contributor Author

Also, any_roots() seems very complicated to me, but maybe it's necessary. I suppose this could be an opportunity to refactor the whole function.

Copy link
Contributor

@grhkm21 grhkm21 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Request a review again after changing, and I'll approve it

src/sage/rings/polynomial/polynomial_element.pyx Outdated Show resolved Hide resolved
@GiacomoPope
Copy link
Contributor Author

I think I will try and refactor any_root as the more I look at it, the less I understand it and think it can be simplified.

@GiacomoPope GiacomoPope changed the title Ugly patch for characteristic two in any_root Total refactor of any_root() to solve issue in characteristic two and clean up the code Jan 27, 2024
@yyyyx4 yyyyx4 self-requested a review January 28, 2024 10:10
src/sage/rings/polynomial/polynomial_element.pyx Outdated Show resolved Hide resolved
src/sage/rings/polynomial/polynomial_element.pyx Outdated Show resolved Hide resolved
src/sage/rings/polynomial/polynomial_element.pyx Outdated Show resolved Hide resolved
src/sage/rings/polynomial/polynomial_element.pyx Outdated Show resolved Hide resolved
src/sage/rings/polynomial/polynomial_element.pyx Outdated Show resolved Hide resolved
src/sage/rings/polynomial/polynomial_element.pyx Outdated Show resolved Hide resolved
src/sage/rings/polynomial/polynomial_element.pyx Outdated Show resolved Hide resolved
src/sage/rings/polynomial/polynomial_element.pyx Outdated Show resolved Hide resolved
src/sage/rings/polynomial/polynomial_element.pyx Outdated Show resolved Hide resolved
@GiacomoPope
Copy link
Contributor Author

@yyyyx4 I'm struggling with the logic of the old version of any_root() and how the ring and degree parameters seem to want to behave.

sage: R.<x> = GF(11)[]
....: g = (x-1)*(x^2 + 3*x + 9) * (x^5 + 5*x^4 + 8*x^3 + 5*x^2 + 3*x + 5)
....: g.any_root(ring=GF(11^10, 'b'), degree=1)
1
sage: g.any_root(ring=GF(11^10, 'b'), degree=2)
5*b^9 + 4*b^7 + 4*b^6 + 8*b^5 + 10*b^2 + 10*b + 5
sage: g.any_root(ring=GF(11^10, 'b'), degree=5)
5*b^9 + b^8 + 3*b^7 + 2*b^6 + b^5 + 4*b^4 + 3*b^3 + 7*b^2 + 10*b
sage: 
sage: g_ext = g.change_ring(GF(11^10, 'b'))
sage: g_ext.any_root(degree=1)
4*b^9 + 8*b^8 + 8*b^7 + 8*b^6 + 3*b^5 + b^4 + 4*b^3 + 4*b^2 + 2*b
sage: g_ext.any_root(degree=2)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [6], line 1
----> 1 g_ext.any_root(degree=Integer(2))

File ~/sage/sage-10.2/src/sage/rings/polynomial/polynomial_element.pyx:2257, in sage.rings.polynomial.polynomial_element.Polynomial.any_root()
   2255 if 2*d > e:
   2256     if degree != e:
-> 2257         raise ValueError("no roots D %s" % self)
   2258     break
   2259 d = d+1

ValueError: no roots D 1
sage: g.any_root(ring=GF(11^10, 'b'), degree=2).parent()
Finite Field in b of size 11^10
sage: g.factor()
(x + 10) * (x^2 + 3*x + 9) * (x^5 + 5*x^4 + 8*x^3 + 5*x^2 + 3*x + 5)
sage: g_ext.factor()
(x + 10) * (x + 6*b^9 + 10*b^8 + 8*b^7 + 9*b^6 + 10*b^5 + 7*b^4 + 8*b^3 + 4*b^2 + b) * (x + 7*b^9 + 3*b^8 + 3*b^7 + 3*b^6 + 8*b^5 + 10*b^4 + 7*b^3 + 7*b^2 + 9*b) * (x + 6*b^9 + 7*b^7 + 7*b^6 + 3*b^5 + b^2 + b + 6) * (x + 5*b^9 + 4*b^7 + 4*b^6 + 8*b^5 + 10*b^2 + 10*b + 8) * (x + 4*b^9 + 9*b^8 + 5*b^7 + 8*b^6 + 9*b^5 + 4*b^3 + 4*b^2 + 9) * (x + 9*b^9 + 8*b^8 + 4*b^7 + 9*b^6 + 8*b^5 + 7*b^4 + 5*b^3 + 3*b + 9) * (x + 7*b^9 + 3*b^8 + 2*b^7 + 4*b^6 + 9*b^5 + 9*b^4 + 9*b^3 + 7*b^2 + 9*b + 9)

So for my new implementation, if someone adds a ring param I assume they want to find factors in this extension, but when the ring and degree params are together, it seems the old one seems to try and find a root in the base field before extension and then casts this the bigger field at the end. This is what allows g.any_root(degree=1) to always return 1, despite the fact that GF(11^10) is the splitting field so any of these roots are "Good".

Should I try and do more work to align these ideas? I dont really see their logic and why what they have here is useful but my doctests fail because over the splitting field there are no degree two irreducible factors so my code breaks.

@GiacomoPope
Copy link
Contributor Author

With some small tweaks to the doctests (root finding is not deterministic anymore, which i guess it never was due to using CZ, but the doctests assumed it was for certain cases?)

The any_root() now seems to behave as it used to. Although making this happen made the code uglier, so I think it would be good to deprecate certain API. I think the case with degree and ring both not None is the weirdest of the choices and could be avoided, but i need opinions.

@GiacomoPope
Copy link
Contributor Author

Added new doctests, some tests fail elsewhere in the codebase mainly as a result of different roots being computed now. I can fix these with the --fixdoctests but I'll await more feedback on my code itself before running this.

@yyyyx4 or @grhkm21 if you can't be bothered to read again, this is roughly ready for review.

@GiacomoPope
Copy link
Contributor Author

Thanks for all the reviews and advice on refactoring. (I think) I even feel happy with how this function works haha!

@vbraun
Copy link
Member

vbraun commented Feb 4, 2024

I'm getting

sage -t --warn-long 40.4 --random-seed=123 src/sage/rings/polynomial/polynomial_quotient_ring_element.py
**********************************************************************
File "src/sage/rings/polynomial/polynomial_quotient_ring_element.py", line 717, in sage.rings.polynomial.polynomial_quotient_ring_element.PolynomialQuotientRingElement.minpoly
Failed example:
    ext(u + 1).minpoly()  # indirect doctest                              # needs sage.modules
Expected:
    x^3 + (396*i + 428)*x^2 + (80*i + 39)*x + 9*i + 178
Got:
    x^3 + (35*i + 428)*x^2 + (351*i + 39)*x + 422*i + 178
**********************************************************************
1 item had failures:
   1 of  19 in sage.rings.polynomial.polynomial_quotient_ring_element.PolynomialQuotientRingElement.minpoly
    [160 tests, 1 failure, 0.51 s]
----------------------------------------------------------------------
sage -t --warn-long 40.4 --random-seed=123 src/sage/rings/polynomial/polynomial_quotient_ring_element.py  # 1 doctest failed
----------------------------------------------------------------------

@@ -105,7 +105,7 @@ cdef class SectionFiniteFieldHomomorphism_givaro(SectionFiniteFieldHomomorphism_
sage: K.<T> = GF(3^4)
sage: f = FiniteFieldHomomorphism_givaro(Hom(k, K))
sage: g = f.section()
sage: g(f(t+1))
sage: g(f(t+1)) # random
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are a great number of added # random tags, are those really random? Isn't this test checking an important identity?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I am uncomfortable with this, but all these doctests seem to depend on the old, undocumented behaviour of .any_root() being deterministic in certain places. Before, the end of the recursion would be a call of f.roots()[0], which fully factors what remains and always returns the same root.

Now, we recursively call f.any_root() and this uses CZ for the splitting and so is not deterministic.

In my opinion "any_root()" should not be deterministic, but if we want these to be fixed then we need to have that the method in hom calls f.roots()[0] instead of the current any_root().

I would love to have more advice on this though!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A lot of these tests really are # random since field embeddings aren't unique: They were deterministic due to a particular implementation choice, but that wasn't really guaranteed anywhere and @GiacomoPope's new algorithm happens to be (1) better, and (2) randomized.
This particular test, however, seems like it should still be deterministic: Here g is supposed to be a left inverse of f by construction, and the test tests that by evaluating the composition on a particular input value.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for my confusion above, I did not notice that I had over-zealously added random tags to these special tests. I have removed this in the latest commit.

@GiacomoPope
Copy link
Contributor Author

GiacomoPope commented Feb 4, 2024

Thank you for reporting the error -- I'll see whether I can find the root of this, I assume it'll be again that any_root() is no longer deterministic

Edit: update to the failure above, this is again due to the randomisation of "any_root" used when computing the extension before computing the minimal polynomial. I don't really know how to proceed... See the above comment re: doctests.

@GiacomoPope
Copy link
Contributor Author

Removed the random tags from those which never needed them. Added a random tag for the min poly over an extension field, as this embedding is not unique and so the result can change.

Fundamental Question

Are we happy with these embeddings being "randomised"? They were not before, but by the accident of how any_root() used to be implemented. This PR fixed a bugged function, but it also now randomised (you really get any root from those from which you ask). This has consequences in many other places (and doctests)

@vbraun
Copy link
Member

vbraun commented Feb 5, 2024

I can live with it being nondeterministic

Ideally we'd not just slap # random on tests but also check something that is true (like the g(f(t+1)) == t+1, or that any root is actually a root), both as clarification for the docs and to guard against regressions.

@GiacomoPope
Copy link
Contributor Author

GiacomoPope commented Feb 5, 2024

Right, I could look at including checks like these in more places, at least for the hom_* files where the # random has appeared the most. If it helps this PR I am happy to take some time to rewrite some of the hom doctests to be better "tests" in light of the change to non-determinism.

Would changes like this be needed in other files too, or simply hom_finite_field*? For the PCL.polynomial I am afraid I have no experience here and do not know "good" doctests

@vbraun
Copy link
Member

vbraun commented Feb 5, 2024

If its easy enough to add then go for it, but if not then not.

@GiacomoPope
Copy link
Contributor Author

Modified the tests to check things like f(a + b) == f(a) + f(b) and also made the section tests more explicit by picking random a and then ensuring g(f(a)) == a.

Did not touch the PCL.polynomial code again, as I'm just not sure what to do.

Copy link

github-actions bot commented Feb 6, 2024

Documentation preview for this PR (built with commit 490ffcb; changes) is ready! 🎉

vbraun pushed a commit to vbraun/sage that referenced this pull request Feb 7, 2024
…aracteristic two and clean up the code

    
This PR was inspired by the issue sagemath#37160, which found that computing
isomorphisms between fields of even characteristic was taking much much
longer than expected (20s vs 100ms).

My first attempt (see commit history) were to simply patch up the
`any_root()` function, but this yielded some results which left the code
in an even more confusing state.

Following the advice of the reviewers, I have now implemented
`any_irreducible_factor()` which given a polynomial element computes an
irreducible factor of this polynomial with minimal degree. When the
optional parameter `degree=None` is set, then an irreducible factor of
degree `degree` is computed and a `ValueError` is raise if no such
factor exists.

Now, `any_root()` becomes simply a call to `self.
any_irreducible_factor(degree=1)` and a root is found from this linear
polynomial. We also handle all the other cases of optional arguments
handled by the old function, so the function *should* behave as before,
but with a cleaner code to read.

### Before PR

```python
┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 10.2, Release Date: 2023-12-03                    │
│ Using Python 3.11.4. Type "help()" for help.                       │
└────────────────────────────────────────────────────────────────────┘
sage: set_random_seed(0)
....: K.<z8> = GF(2^8)
....: R.<x> = K[]
....: u = R.irreducible_element(2)
....: K_ext = K.extension(modulus=u, names="a")
....: R_ext.<y_ext> = K_ext[]
....: poly = R_ext.random_element(2).monic()
....: %time poly.roots()
CPU times: user 20.5 s, sys: 19.9 ms, total: 20.5 s
Wall time: 20.5 s
[]
```

### After PR

```py
┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 10.3.beta6, Release Date: 2024-01-21              │
│ Using Python 3.11.4. Type "help()" for help.                       │
└────────────────────────────────────────────────────────────────────┘
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Warning: this is a prerelease version, and it may be unstable.     ┃
┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛
sage: set_random_seed(0)
....: K.<z8> = GF(2^8)
....: R.<x> = K[]
....: u = R.irreducible_element(2)
....: K_ext = K.extension(modulus=u, names="a")
....: R_ext.<y_ext> = K_ext[]
....: %time R_ext.random_element(2).roots()
CPU times: user 110 ms, sys: 9.03 ms, total: 119 ms
Wall time: 150 ms
[]
```
This fixes sagemath#37160 but i think more feedback and thorough doctests need
to be included considering this is a very old function which many
projects will be using.
    
URL: sagemath#37170
Reported by: Giacomo Pope
Reviewer(s): Giacomo Pope, grhkm21, Lorenz Panny, Volker Braun
@vbraun vbraun merged commit 7c0a5ec into sagemath:develop Feb 13, 2024
17 of 18 checks passed
vbraun pushed a commit to vbraun/sage that referenced this pull request Feb 28, 2024
    
This pull request aims to fix some issues which came with the
refactoring of `any_root()` in
sagemath#37170.

I am afraid that apart from a little more cleaning up and verbose
comments, this "fix" really relies on `f.roots()` rather than
`f.any_root()` when roots over extension fields are computed.

The core problem currently is that the proper workflow should be the
following:

- Given a polynomial $f$, find an irreducible factor of smallest degree.
When `ring` is `None` and `degree` is `None`, only linear factors are
searched for. This code is unchanged and works fast.
- When `degree` is `None` and `ring` is not `None`, the old code used to
perform `f.change_ring(ring).any_root()` and look for a linear factor in
`ring`. The issue is when `ring` is a field with a large degree,
arithmetic in this ring is very slow and root finding is very slow.
- When `degree` is not `None` then an irreducible of degree `degree` is
searched for, $f$, and then a root is found in an extension ring. This
is either the user supplied `ring`, or it is in the extension
`self.base_ring().extension(degree)`. Again, as this extension could be
large, this can be slow.

Now, the reason that there's a regression in speed, as explained in
sagemath#37359, is that the method before
refactoring simply called `f.roots(ring)[0]` when working for these
extensions. As the root finding is always performed by a C library this
is much faster than the python code which we have for `any_root()` and
so the more "pure" method I wrote simply is slower.

The real issue holding this method back is that if we are in some field
$GF(p^k)$ and `ring` is some extension $GF(p^{nk})$ then we are not
guaranteed a coercion between these rings with the current coercion
system. Furthermore, if we extend the base to some $GF(p^{dk})$ with
$dk$ dividing $nk$ we also have no certainty of a coercion from this
minimal extension to the user supplied ring.

As a result, a "good" fix is delayed until we figure out finite field
coercion. **Issue to come**.

Because of all of this, this PR goes back to the old behaviour, while
keeping the refactored and easier to read code. I hope one day in the
future to fix this.

Fixes sagemath#37359
Fixes sagemath#37442
Fixes sagemath#37445
Fixes sagemath#37471

**EDIT**

I have labeled this as critical as the slow down in the most extreme
cases causes code to complete hang when looking for a root. See for
example: sagemath#37442. I do not want
sagemath#37170 to be merged into 10.3
without this patch.

Additionally, a typo meant that a bug was introduced in the CZ
splitting. This has also been fixed.

Also, also in the refactoring the function behaved as intended (rather
than the strange behaviour from before) which exposed an issue
sagemath#37471 which I have fixed.
    
URL: sagemath#37443
Reported by: Giacomo Pope
Reviewer(s): Giacomo Pope, grhkm21, Lorenz Panny, Travis Scrimshaw
@GiacomoPope GiacomoPope deleted the fix_any_root_even_char branch March 6, 2024 16:25
@mkoeppe mkoeppe added this to the sage-10.3 milestone Mar 7, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Slow root finding over GF(2^k)[] quotient rings
5 participants