Improved nth root for finite fields and integer_mods

roed314 commented 14 years ago

Implements an algorithm described in

Johnston, Anna M. A generalized qth root algorithm. 
Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms. 
Baltimore, 1999: pp 929-930.

This means we can take nth roots with large n, since we no longer need to create the polynomial xⁿ - a.

CC: @robertwb

Component: algebra

Keywords: finite fields, nth root

Author: David Roe

Reviewer: Paul Zimmermann, David Loeffler, Bill Hart

Merged: sage-4.6.2.alpha4

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

roed314 commented 14 years ago

Attachment: 7931.patch.gz

roed314 commented 14 years ago

Description changed:

--- 
+++ 
@@ -1,8 +1,9 @@
 Implements an algorithm described in 
-{{
+
+```
 Johnston, Anna M. A generalized qth root algorithm. 
 Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms. 
 Baltimore, 1999: pp 929-930.
-}}}
+```

-This means we can take nth roots with large n, since we no longer need to create the polynomial x^n-a.
+This means we can take nth roots with large n, since we no longer need to create the polynomial x<sup>n</sup> - a.

b12b007b-ced8-4aef-b4e9-5b123b0bb8b3 commented 14 years ago

comment:2

Is there something I am missing?

sage: K = GF(31)
sage: b=K(22)^200
sage: b.nth_root(200)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)

/virtual/scratch/rishi/sage-4.3.1.alpha2-sage.math.washington.edu-x86_64-Linux/<ipython console> in <module>()

/virtual/scratch/rishi/sage-4.3.1.alpha2-sage.math.washington.edu-x86_64-Linux/local/lib/python2.6/site-packages/sage/rings/integer_mod.so in sage.rings.integer_mod.IntegerMod_abstract.nth_root (sage/rings/integer_mod.c:11446)()

ValueError: no nth root
sage:

JohnCremona commented 14 years ago

comment:3

It looks to me as if the original code (now deleted) worked with finite field elements, while thie new code works with integer-mods. That must surely mean that for non-prime finite fields there was something which worked, but now there isn't? If I am right, it would be a good idea to allow the original code as a fall-back.

JohnCremona commented 14 years ago

comment:4

Patch applies fine to 4.3.1 and tests pass, but I reproduced the example from rishi, so -- "needs work".

roed314 commented 14 years ago

Attachment: 7931_nth_root.patch.gz

Should fix the problem observed earlier

zimmermann6 commented 14 years ago

comment:6

David, should the 2nd patch be applied after the 1st one or alone?

roed314 commented 14 years ago

comment:7

Alone. Sorry about that.

roed314 commented 14 years ago

comment:8

Oops. I made it depend on a bunch of other changes

8218 -> 8332 -> 7880 -> 7883 -> 8333 -> 8334 -> this patch

The actual change is fairly small; I'll try to extract it out and get something based only on 4.3.3, but I won't be able to do that tonight.

zimmermann6 commented 14 years ago

comment:9

The actual change is fairly small; I'll try to extract it out and get something based only on 4.3.3, but I won't be able to do that tonight.

that would be nice. Otherwise we can wait for the other patches to be reviewed. Paul

aghitza commented 14 years ago

comment:10

Replying to @roed314:

Oops. I made it depend on a bunch of other changes
8218 -> 8332 -> 7880 -> 7883 -> 8333 -> 8334 -> this patch
The actual change is fairly small; I'll try to extract it out and get something based only on 4.3.3, but I won't be able to do that tonight.

Just so we have trac's help to see which of the dependencies have already been merged: this depends on #8218, #8332, #7880, #7883, #8333, #8334.

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

comment:11

Patch 7931_nth_root.patch applies and builds fine and all doctests in finite_rings pass (with the dependencies installed); but I'm not terribly happy about the amount of code duplication that's going on -- the code added to element_base and integer_mod has about 30 lines that are virtually identical. Could it not be refactored to avoid this?

roed314 commented 13 years ago

comment:12

I completely agree. In the long run I want to change FiniteFieldElement to FiniteRingElement and make IntegerMod_abstract inherit from that. The nth_root function on FiniteRingElement can then be modified to handle both cases. An additional benefit of having such a common parent is that the finite field elements can then no longer require p prime (or modulus irreducible for that matter), which will be useful for p-adic extensions. In fact, I want to modify the finite field element classes to use templates (a la sage.rings.polynomial.polynomial_template): then we can share common code between polynomials, finite fields and rings, and p-adic extension fields. Coercions between them will be much easier and faster and it will make implementing extensions of extensions easier (which are necessary for p-adic extensions that are neither unramified nor totally ramified).

The first step though, of making IntegerMod_abstract inherit from FiniteRingElement caused compile-time bugs that I couldn't figure out and I gave up for the moment.

So currently their closest common superclass is CommutativeRingElement. It didn't seem like a good idea to put the common nth_root code there, and there wasn't another obvious way to do it. Any suggestions?

roed314 commented 13 years ago

comment:13

Another problem with the current patch is that for large n, finding the factorization of the size of the unit group becomes a bottleneck. There's a change in #8335 that helps with this problem for finite fields of small characteristic by caching the factorization of pⁿ-1 and using the Cunningham package at #7240 to speed up the computation in the first place. Thematically it would make sense to backport it to this patch (which I've done already locally), but it makes this patch depend on #7240 (otherwise there are warnings printed when _factor_cunningham is used). What do you think?

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

comment:14

Could one perhaps have an auxilliary helper function (not a class method) that both of the nth root methods call? This is not terribly elegant, but it solves the duplication issue.

I'd argue for not making this patch depend on #7240, since that patch seems to be stalled at the moment.

roed314 commented 13 years ago

Attachment: 7931_nth_root.2.patch.gz

Apply just this patch

roed314 commented 13 years ago

comment:15

I factored out most of the common code, and backported a few changes from #8335. It still doesn't depend on #7240.

roed314 commented 13 years ago

Apply on top of 7931_nth_root.2.patch; moves _nth_root_common to FiniteRingElement, a new superclass of IntegerMod_abstract and the finite field elements

zimmermann6 commented 13 years ago

comment:16

Attachment: 7931_common_superclass.patch.gz

I tried applying both patches in sage-4.6.alpha1, but got a failure when running sage -br:

python `which cython` --embed-positions --directive cdivision=True -I/tmp/sage-4.6.alpha1/devel/sage-7931 -o sage/rings/finite_rings/element_ntl_gf2e.cpp sage/rings/finite_rings/element_ntl_gf2e.pyx

Error converting Pyrex file to C:
------------------------------------------------------------
...
        if PY_TYPE_CHECK(e, int) \
               or PY_TYPE_CHECK(e, long) or PY_TYPE_CHECK(e, Integer):
            GF2E_conv_long(res.x,int(e))
            return res

        if PY_TYPE_CHECK(e, FiniteFieldElement) or \
                                             ^
------------------------------------------------------------

/tmp/sage-4.6.alpha1/devel/sage-7931/sage/rings/finite_rings/element_ntl_gf2e.pyx:515:46: undeclared name not builtin: FiniteFieldElement
Error running command, failed with status 256.
sage: There was an error installing modified sage library code.

Which version was used to produce the patches?

Paul

zimmermann6 commented 13 years ago

Reviewer: Paul Zimmermann

zimmermann6 commented 13 years ago

comment:17

I realized that I should also apply the patches #7883 and #8334, which were included in sage-4.6.alpha2, but were not yet in sage-4.6.alpha1. However after importing #7883 successfully in sage-4.6.alpha1, importing #8334 gives:

sage: hg_sage.import_patch("/tmp/8333_8334_ALL-better_commit_string.patch")
cd "/tmp/sage-4.6.alpha1/devel/sage" && hg status
cd "/tmp/sage-4.6.alpha1/devel/sage" && hg status
cd "/tmp/sage-4.6.alpha1/devel/sage" && hg import   "/tmp/8333_8334_ALL-better_commit_string.patch"
applying /tmp/8333_8334_ALL-better_commit_string.patch
patching file sage/rings/ideal_monoid.py
Hunk #1 succeeded at 90 with fuzz 2 (offset 0 lines).
patching file sage/rings/residue_field.pyx
Hunk #3 succeeded at 74 with fuzz 2 (offset 0 lines).
Hunk #15 FAILED at 624
1 out of 36 hunks FAILED -- saving rejects to file sage/rings/residue_field.pyx.rej
abort: patch failed to apply

and I'm stuck there. Maybe I should wait that sage-4.6.alpha2 is out.

Paul

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

comment:18

Ticket #8334 has some other dependencies -- they are listed on that ticket (Jeroen's comment 24). With those installed the patch applies fine, and all doctests in sage/rings pass (I'm running a full ptestlong cycle now), so I'm putting this back to "needs_review".

By the way, Mitesh has a lovely script merger-4.6.alpha2, which you can run on a clean 4.6.alpha1 build and it will download and apply all of the patches and spkgs so far merged.

zimmermann6 commented 13 years ago

comment:19

I tried Mitesh's script on my sage-4.6.alpha1 build (after sage -b main), then I applied the two patches from that ticket, but I still get the same error as in comment 16.

Paul

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

comment:20

That's bizarre, because the patches apply and build fine for me.

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

comment:21

Actually I hadn't myself tried Mitesh's script when I posted above; I just did, and I couldn't get it to work either. But it should work if you install:

9898_pari_decl.patch
9753.patch
9764_ideal_repr_new.patch
trac_7883-ideals-folded.patch
8333_8334_ALL-rebased_for_9764.patch
7931_nth_root.2.patch
7931_common_superclass.patch

zimmermann6 commented 13 years ago

comment:22

I managed to apply the patches following comment 21, however the following seems incorrect to me:

sage: b=Integers(3)(2)
sage: b.nth_root(2)
1

whereas in say Sage 4.4.4 we got ValueError: no nth root.

PS: I used the following reviewer program. Feel free to add it to the doctests.

for n in range(2,100):
   K=Integers(n)
   for e in range(1,100):
      for a in range(1,n):
         b = K(a)
         r = b.nth_root(e)
         if r^e <> b:
            print n, e, a
            raise ValueError

roed314 commented 13 years ago

comment:23

Thanks for catching that. I've fixed the problem and added a _nth_root_naive and check that the output of b.nth_root(e, all=True) matches b._nth_root_naive(e) for all b in Zmod(n) for 2 <= n < 100 and 1 <= e < 100. A shortened version of this test is left in as a doctest for _nth_root_naive.

roed314 commented 13 years ago

Apply on top of 7931_nth_root.2.patch and 7931_common_superclass.patch

zimmermann6 commented 13 years ago

comment:24

Attachment: 7931_fixes.patch.gz

with the new patch, I still get unexpected results:

sage: K=Integers(6)
sage: b=K(3)
sage: b.nth_root(0,all=True)
[3]

Paul

roed314 commented 13 years ago

Attachment: 7931_fix2.patch.gz

Apply on top of previous patches

roed314 commented 13 years ago

comment:25

Unexpected indeed. Thanks for checking all these corner cases that I neglected (though some of the earlier ones were bugs I should have found if I'd tested sufficiently).

In this case I'm not quite sure what the right output is: an empty list or a ZeroDivisionError? I've gone for the empty list, on the theory that b.nth_root(e) should return a list of all elements a with a^e=b, and this makes sense even if e=0.

JohnCremona commented 13 years ago

comment:26

e=0 is a hard case since there are too many solutions -- do we want to return all a in the ring? or all invertible a?

Similarly, do we allow negative e, where b must be invertible and b.nth_root(e) == (1/b).nth_root(-e), with an error if b is not invertible?

Or do we want to specify in the documentation that e must be positive, and raise an error if not?

zimmermann6 commented 13 years ago

comment:27

My 2 cents:

for e=0 and b<>1, b.nth_root(e) should return the empty list
for e=0 and b=1, b.nth_root(e) should return only one solution, for example 1 (with all=True, return all a in the ring)

Paul

zimmermann6 commented 13 years ago

comment:28

sorry, another problem:

sage: n=13777831336991389951
sage: b=3798677550250515336
sage: e=10608321776141318019
sage: K = Integers(n)
sage: b=K(b)
sage: b.nth_root(e)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)

The documentation does not say when ZeroDivisionError can be obtained. This was found with the following program:

while True:
   n = randint(0,2^64)
   K = Integers(n)
   b = K.random_element()
   e = randint(0,2^64)
   print n, b, e, a
   sys.stdout.flush()
   try:
      a = b.nth_root(e)
      if a^e <> b:
         print n, b, e, a
         raise NotImplementedError
   except ValueError:
      n = 0

roed314 commented 13 years ago

comment:29

You've found a bug in crt: it claims to work as long as the moduli are relatively prime, but in fact would often fail if one of them was 1. I fixed that and clarified the behavior of nth_root when n<=0 (it either returns the empty list or raises a ValueError, depending on the value of all; mod(1,n).nth_root(0) returns the list of all nonzero elements modulo n).

zimmermann6 commented 13 years ago

comment:30

Replying to @roed314:

You've found a bug in crt: it claims to work as long as the moduli are relatively prime, but in fact would often fail if one of them was 1. I fixed that and clarified the behavior of nth_root when n<=0 (it either returns the empty list or raises a ValueError, depending on the value of all; mod(1,n).nth_root(0) returns the list of all nonzero elements modulo n).

Hi David,

I'm not very happy with that answer. If a bug in crt was found, I would expect you to show a concrete example, to report it as a new ticket, and to mention in your 3rd patch that it depends on the new ticket, and can be removed once the new ticket is fixed.

Paul

roed314 commented 13 years ago

Attachment: 7931_fix3.patch.gz

Apply on top of previous patches

roed314 commented 13 years ago

comment:31

You're right that it needs doctests; I've added some. I also changed 7931_fix3.patch by updating the nth_root function in element_base to match the behavior of the one in integer_mod for non-positive inputs.

Personally, I don't think that the crt bug is a major enough problem to need it's own ticket. It only applies if one of the arguments is mod(0,1). But I've separated the fix into its own patch; if you think it needs it's own ticket would you be willing to make the ticket and move the patch over there?

zimmermann6 commented 13 years ago

comment:32

I don't understand the crt patch: the examples given already worked in 4.6.alpha1 (and similarly in 4.4.4):

----------------------------------------------------------------------
| Sage Version 4.6.alpha1, Release Date: 2010-09-18                  |
| Type notebook() for the GUI, and license() for information.        |
----------------------------------------------------------------------
**********************************************************************
*                                                                    *
* Warning: this is a prerelease version, and it may be unstable.     *
*                                                                    *
**********************************************************************
sage: mod(0,1).crt(mod(4,17)) 
4
sage: mod(0,1).crt(mod(0,1)) 
0
sage: mod(21,22).crt(mod(0,1))
21

What is exactly the crt bug (if any)?

Paul

roed314 commented 13 years ago

comment:33

Ah. Now I understand why this problem only appeared with your latest test program that used large moduli. The issue only occurs for mod(0,1).crt(a) when a has type IntegerMod_gmp. So we didn't see it for smaller moduli.

I've updated 7931_crt.patch, and the doctest there does fail on my unpatched 4.6.alpha1.

zimmermann6 commented 13 years ago

comment:34

I understand now, this is really a bug, and I think it should be considered in a separate ticket:

sage: mod(0,1).crt(mod(4,2^31-2)) 
4
sage: mod(0,1).crt(mod(4,2^31-1)) 
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)

roed314 commented 13 years ago

comment:35

This is now #10047.

roed314 commented 13 years ago

Attachment: 7931_crt.patch.gz

This is now #10047.

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

Attachment: 7931_nth_root-folded.patch.gz

Qfolded patch. Apply only this patch.

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

comment:36

It took me some while to understand what patches to apply in what order, so having done so I thought it might be helpful to qfold everything into one patch. (I'm pleasantly surprised that patchbot understood immediately.)

Now I'll have a look at the code.

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

comment:37

The code looks plausible, and I'm pleased to report that we aren't going to have the same problem as #9304 -- old pickled objects seem to unpickle just fine. I'm not qualified to comment on the details of the algorithm though. Paul?

zimmermann6 commented 13 years ago

comment:38

Replying to @loefflerd:

The code looks plausible, and I'm pleased to report that we aren't going to have the same problem as #9304 -- old pickled objects seem to unpickle just fine. I'm not qualified to comment on the details of the algorithm though. Paul?

I won't have much time in the near future to look closely at the algorithm. Unless someone else has time before me (John?) I'll look at this later. Don't hesitate to ping me.

Paul

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

Fixed ReST formatting

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

Changed reviewer from Paul Zimmermann to Paul Zimmermann, David Loeffler, Bill Hart

11d1fc49-71a1-44e1-869f-76be013245a0 commented 13 years ago

comment:39

Attachment: 7931_nth_root-folded.2.patch.gz

I checked this one with Bill Hart, who knows much more about these things than I, and he reckons that the implementation looks correct; I've run the tests and it seems to be fine. My only contribution has been to qfold everything and adjust some of the reference manual formatting, so I guess that doesn't need a separate review.

Apply only the last patch.

sagemath / sage

Improved nth root for finite fields and integer_mods #7931