Opened 10 years ago
Closed 9 years ago
#7931 closed enhancement (fixed)
Improved nth root for finite fields and integer_mods
Reported by: | roed | Owned by: | roed |
---|---|---|---|
Priority: | major | Milestone: | sage-4.6.2 |
Component: | algebra | Keywords: | finite fields, nth root |
Cc: | robertwb | Merged in: | sage-4.6.2.alpha4 |
Authors: | David Roe | Reviewers: | Paul Zimmermann, David Loeffler, Bill Hart |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
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.
Attachments (10)
Change History (51)
Changed 10 years ago by
comment:1 Changed 10 years ago by
- Description modified (diff)
- Status changed from new to needs_review
comment:2 Changed 10 years ago by
comment:3 Changed 10 years ago by
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.
comment:4 Changed 10 years ago by
- Status changed from needs_review to needs_work
Patch applies fine to 4.3.1 and tests pass, but I reproduced the example from rishi, so -- "needs work".
comment:5 Changed 9 years ago by
- Status changed from needs_work to needs_review
comment:6 Changed 9 years ago by
- Status changed from needs_review to needs_info
David, should the 2nd patch be applied after the 1st one or alone?
comment:7 Changed 9 years ago by
- Status changed from needs_info to needs_review
Alone. Sorry about that.
comment:8 follow-up: ↓ 10 Changed 9 years ago by
- Owner changed from AlexGhitza to roed
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.
comment:9 Changed 9 years ago by
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
comment:10 in reply to: ↑ 8 Changed 9 years ago by
Replying to roed:
Oops. I made it depend on a bunch of other changes
8218 -> 8332 -> 7880 -> 7883 -> 8333 -> 8334 -> this patchThe 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.
comment:11 Changed 9 years ago by
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?
comment:12 Changed 9 years ago by
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?
comment:13 Changed 9 years ago by
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^{n}-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?
comment:14 Changed 9 years ago by
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.
comment:15 Changed 9 years ago by
Changed 9 years ago by
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
comment:16 Changed 9 years ago by
- Reviewers set to Paul Zimmermann
- Status changed from needs_review to needs_info
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
comment:17 Changed 9 years ago by
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
comment:18 Changed 9 years ago by
- Status changed from needs_info to needs_review
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.
comment:19 Changed 9 years ago by
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
comment:20 Changed 9 years ago by
That's bizarre, because the patches apply and build fine for me.
comment:21 Changed 9 years ago by
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
comment:22 Changed 9 years ago by
- Status changed from needs_review to needs_work
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
comment:23 Changed 9 years ago by
- Status changed from needs_work to needs_review
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
.
comment:24 Changed 9 years ago by
- Status changed from needs_review to needs_work
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
comment:25 Changed 9 years ago by
- Status changed from needs_work to needs_review
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
.
comment:26 Changed 9 years ago by
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?
comment:27 Changed 9 years ago by
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
comment:28 Changed 9 years ago by
- Status changed from needs_review to needs_work
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
comment:29 follow-up: ↓ 30 Changed 9 years ago by
- Status changed from needs_work to needs_review
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).
comment:30 in reply to: ↑ 29 Changed 9 years ago by
- Status changed from needs_review to needs_work
Replying to roed:
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
whenn<=0
(it either returns the empty list or raises aValueError
, depending on the value ofall
;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
comment:31 Changed 9 years ago by
- Status changed from needs_work to needs_review
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?
comment:32 Changed 9 years ago by
- Cc robertwb added
- Status changed from needs_review to needs_info
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
comment:33 Changed 9 years ago by
- Status changed from needs_info to needs_review
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
.
comment:34 Changed 9 years ago by
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)
comment:35 Changed 9 years ago by
This is now #10047.
comment:36 Changed 9 years ago by
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.
comment:37 follow-up: ↓ 38 Changed 9 years ago by
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?
comment:38 in reply to: ↑ 37 Changed 9 years ago by
Replying to davidloeffler:
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
comment:39 Changed 9 years ago by
- Reviewers changed from Paul Zimmermann to Paul Zimmermann, David Loeffler, Bill Hart
- Status changed from needs_review to positive_review
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.
comment:40 Changed 9 years ago by
Great! Thanks to all of the reviewers for looking at this.
comment:41 Changed 9 years ago by
- Merged in set to sage-4.6.2.alpha4
- Resolution set to fixed
- Status changed from positive_review to closed
Is there something I am missing?