Opened 12 years ago
Closed 12 years ago
#10836 closed defect (fixed)
primitive root is broken
Reported by: | Karl-Dieter Crisman | Owned by: | William Stein |
---|---|---|---|
Priority: | critical | Milestone: | sage-4.7 |
Component: | number theory | Keywords: | |
Cc: | William Stein, Alyson Deines | Merged in: | sage-4.7.alpha4 |
Authors: | Karl-Dieter Crisman, William Stein, Jeroen Demeyer | Reviewers: | William Stein, Karl-Dieter Crisman, Jeroen Demeyer |
Report Upstream: | Reported upstream. Developers deny it's a bug. | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
Pari has a major bug in its primitive root:
sage: primitive_root(15) 2 sage: mod(2,15).multiplicative_order() 4 sage: euler_phi(15) 8 sage: [mod(2,15)^i for i in [1..8]] [2, 4, 8, 1, 2, 4, 8, 1] sage: b = pari(15) sage: b.znprimroot() Mod(2, 15)
Note that this DOES work correctly with whatever Pari is in Sage 4.4.4.
--------------------------------------------------------------------- | Sage Version 4.4.4, Release Date: 2010-06-23 | | Type notebook() for the GUI, and license() for information. | ---------------------------------------------------------------------- Loading Sage library. Current Mercurial branch is: hackbranch sage: primitive_root(15) <snip> ArithmeticError: There is no primitive root modulo n
Apply trac_10836.2.patch and 10836_additional.patch.
Attachments (4)
Change History (34)
comment:1 Changed 12 years ago by
Cc: | William Stein added |
---|---|
Report Upstream: | Not yet reported upstream; Will do shortly. → Reported upstream. Developers deny it's a bug. |
comment:2 Changed 12 years ago by
Ergo, we will have to catch such cases in an effort to remain mathematically correct.
comment:3 Changed 12 years ago by
Pretty easy to fix but there are some choices regarding speed and whether we now believe in pari's ispower function.
trying: 2 call Integers(n).multiplicative_generator() : 625 loops, best of 3: 3.97 µs per loop pari after Integers(n).multiplicative_group_is_cyclic(): 625 loops, best of 3: 16.1 µs per loop pari after inline Integers(n).mult.._group_is_cyclic() : 625 loops, best of 3: 9.53 µs per loop pari without check: : 625 loops, best of 3: 8.74 µs per loop pari trusting that pari.ispower now works : 625 loops, best of 3: 9.39 µs per loop I.mul_gen() if prime, trusting pari otherwise: : 625 loops, best of 3: 29.3 µs per loop trying: 97 call Integers(n).multiplicative_generator() : 625 loops, best of 3: 3.93 µs per loop pari after Integers(n).multiplicative_group_is_cyclic(): 625 loops, best of 3: 33 µs per loop pari after inline Integers(n).mult.._group_is_cyclic() : 625 loops, best of 3: 27.3 µs per loop pari without check: : 625 loops, best of 3: 10.8 µs per loop pari trusting that pari.ispower now works : 625 loops, best of 3: 19.8 µs per loop I.mul_gen() if prime, trusting pari otherwise: : 625 loops, best of 3: 29.4 µs per loop trying: 10007 call Integers(n).multiplicative_generator() : 625 loops, best of 3: 3.93 µs per loop pari after Integers(n).multiplicative_group_is_cyclic(): 625 loops, best of 3: 34.1 µs per loop pari after inline Integers(n).mult.._group_is_cyclic() : 625 loops, best of 3: 28.4 µs per loop pari without check: : 625 loops, best of 3: 11.4 µs per loop pari trusting that pari.ispower now works : 625 loops, best of 3: 21.3 µs per loop I.mul_gen() if prime, trusting pari otherwise: : 625 loops, best of 3: 30 µs per loop trying: 1000000007 call Integers(n).multiplicative_generator() : 625 loops, best of 3: 3.86 µs per loop pari after Integers(n).multiplicative_group_is_cyclic(): 625 loops, best of 3: 78.6 µs per loop pari after inline Integers(n).mult.._group_is_cyclic() : 625 loops, best of 3: 73.6 µs per loop pari without check: : 625 loops, best of 3: 53.7 µs per loop pari trusting that pari.ispower now works : 625 loops, best of 3: 71.1 µs per loop I.mul_gen() if prime, trusting pari otherwise: : 625 loops, best of 3: 31.6 µs per loop trying: 17^5 call Integers(n).multiplicative_generator() : 625 loops, best of 3: 79.7 µs per loop pari after Integers(n).multiplicative_group_is_cyclic(): 625 loops, best of 3: 57.4 µs per loop pari after inline Integers(n).mult.._group_is_cyclic() : 625 loops, best of 3: 51 µs per loop pari without check: : 625 loops, best of 3: 10.5 µs per loop pari trusting that pari.ispower now works : 625 loops, best of 3: 20.8 µs per loop I.mul_gen() if prime, trusting pari otherwise: : 625 loops, best of 3: 46.7 µs per loop trying: 2 * 109^100 call Integers(n).multiplicative_generator() : 625 loops, best of 3: 857 µs per loop pari after Integers(n).multiplicative_group_is_cyclic(): 625 loops, best of 3: 879 µs per loop pari after inline Integers(n).mult.._group_is_cyclic() : 625 loops, best of 3: 869 µs per loop pari without check: : 625 loops, best of 3: 49.4 µs per loop pari trusting that pari.ispower now works : 625 loops, best of 3: 95.1 µs per loop I.mul_gen() if prime, trusting pari otherwise: : 625 loops, best of 3: 123 µs per loop
Long story short, Integers(n).multiplicative_generator is very fast for primes, but otherwise slow. Trying to mix the two (using multiplicative_generator for primes and pari for everything else) winds up being slower than pari for small primes and faster for large ones; could branch on size but not sure if the speed benefit would be worth the complexity. I don't have a good sense for Python overheads at the microsecond level.
comment:4 Changed 12 years ago by
From some correspondence with myself, William, and the Pari team, where it becomes clear that their implementation has the rather useful virtue of being very very fast:
> Compare : > > (23:07) gp > p = 10^50 + 151; > (23:07) gp > for(i=1,10^2, isprime(p)) > time = 1,324 ms. > (23:07) gp > for(i=1,10^2, ispseudoprime(p)) > time = 32 ms. > (23:07) gp > for(i=1,10^2, znprimroot(p)) > time = 112 ms. > > For simple inputs (where factoring p-1 is not a problem), it's orders of > magnitude faster to find a plausible primitive root than to guarantee > that the question makes sense (primality proof). In fact, already trying > to weed out most composites [ ispseudoprime vs. isprime ] takes a > non-negligible amount of time. > > As p increases, I can produce (contrived but) legitimate prime inputs > where isprime() takes an arbitrary amount of time whereas znprimroot() > remains essentially instantaneous. For random large p, factoring p-1 will > dominate the running time, hiding all this. (In fact, znprimroot won't > even stand a chance of finishing.)
So we will keep using their very fast function and wrap it with something that checks for sensible input, but adding documentation like
sage: b = pari(11) sage: b.znprimroot()
to the documentation for our primitive_root
function for those needing direct access to the speed, warning of 'undefined' input. William, does that sound right?
comment:5 Changed 12 years ago by
Authors: | → Karl-Dieter Crisman |
---|---|
Milestone: | → sage-4.7 |
Status: | new → needs_review |
Well, I have a provisional patch, but because I am very unfamiliar with how gp/GP/Pari/PARI interact with Sage, this might be really kludgy. Also, I have observed the following perplexing behavior in testing (in Sage or Pari):
parisize = 8000000, primelimit = 500509 ? p = 10^80+129 %1 = 100000000000000000000000000000000000000000000000000000000000000000000000000000129 ? znprimroot(p) ^C *** at top-level: znprimroot(p) *** ^------------- *** znprimroot: user interrupt after 29,777 ms. *** Break loop: type <Return> to continue; 'break' to go back to GP break> break ? p = 10^76+133 %2 = 10000000000000000000000000000000000000000000000000000000000000000000000000133 ? znprimroot(p) %3 = Mod(2, 10000000000000000000000000000000000000000000000000000000000000000000000000133)
Should this size difference really make such a huge difference? I would estimate the second computation took about a second. The slowness is all in Pari, too; in Sage I get the same thing. I would think these were not such large numbers.
sage: primitive_root(next_prime(3.32*10^76)) ^C KeyboardInterrupt: sage: primitive_root(next_prime(3.33*10^76)) 2
This also happens in Sage 4.6.2, so it's nothing unusual. Am I just getting unlucky with numbers that have really big primitive roots? Maybe; I just put a problem to that effect on a take-home exam...
On the other hand, when I quit after having done these breaks, I get
*** Warning: I/O: leaked file descriptor (0): /var/folders/Yy/YytEJm5VEB0+pBRD7JNLe++++TQ/-Tmp-/MPQS.503.640ad/LPNEW. *** Warning: I/O: leaked file descriptor (0): /var/folders/Yy/YytEJm5VEB0+pBRD7JNLe++++TQ/-Tmp-/MPQS.503.640ad/FNEW. *** Warning: I/O: leaked file descriptor (0): /var/folders/Yy/YytEJm5VEB0+pBRD7JNLe++++TQ/-Tmp-/MPQS.503.640ac/FNEW.
type stuff, so maybe I'm not crazy.
Anyway, for the purposes here, the problem is finding a doctest where the Sage version takes a long time without the Pari one taking a long time. I tried lots and lots of things.
Any ideas? Anyway, maybe this patch is better than WRONG answers, so here is a patch to look at.
Changed 12 years ago by
Attachment: | trac_10836-primitive.patch added |
---|
comment:6 Changed 12 years ago by
sage: for i in range(10^2): ....: n = randrange(10^50) ....: try: ....: primitive_root(n) ....: except: ....: pass
I've gotten big error rates on this in the old one - 70% or more of the primitive roots were erroneous. Something like this takes 11 seconds in both versions, though.
But there definitely is a slowdown for some numbers:
sage: sage: timeit('primitive_root(10^50+151)') 25 loops, best of 3: 17.6 ms per loop sage: sage: timeit('pari(10^50+151).znprimroot()') 625 loops, best of 3: 1.46 ms per loop
Before, both took about the same amount of time.
Another data point, possibly irrelevant:
sage: timeit('try:\n primitive_root(randrange(10^60))\nexcept:\n pass') 5 loops, best of 3: 97.9 ms per loop sage: timeit('try:\n pari(randrange(10^60)).znprimroot()\nexcept:\n pass') 5 loops, best of 3: 200 ms per loop
The variation is extremely wide on this, as one might expect, and sometimes one is faster, sometimes another. I've gotten from 25 to 300 ms for both commands, irrespective of before/after the patch. I did once get 1.85 seconds for the new primitive_root, so it seems it is sometimes slower, but I also got as long as 1.24 s for the pari version, and neither of these bad times came close to always happening - much more important were the particular random numbers chosen.
comment:7 Changed 12 years ago by
Cc: | Alyson Deines added |
---|
Hopefully someone can review this during Sage Day 29? It's mathematically wrong currently...
comment:8 Changed 12 years ago by
I think maybe there's a little too much duplication with all the branches? ISTM it's simpler to implement the condition exactly from the definition given, with something like
if ((n == 4) or (n % 2 == 1 and n.is_prime_power()) or (n % 2 == 0 and n % 4 != 0 and (n//2).is_prime_power())): return ans
where each condition group matches exactly one of the classes of allowed arguments in the docstring. I find this is also marginally faster on average for successful calls. But I think we should check n for the existence of a primitive root before we call pari, because when testing it I came across this (4.6.2):
sage: factor(1729) 7 * 13 * 19 sage: time primitive_root(1729) [... several minutes go by, so I have no idea if this eventually returns a value or not ...]
Similarly for 2465, 3458, 4930. This is a pari bug not an interface bug, as working directly with the console shows the same issue.
comment:9 Changed 12 years ago by
If you can find a quicker way to write that, this is fine. I was hoping to speed things up by checking certain cases first, such as primes. Branching can be confusing, you are right.
More precisely, it's what Pari claims is user input error. I actually explicitly didn't do this because we didn't want to slow down the checking of primitive roots for legitimate input, which this would do for huge numbers - that is the whole idea behind Pari's thinking.
At this point, I'd appreciate the input of someone who actually needs primitive roots of huge numbers! The question is whether we want the bottleneck with legitimate or non-legitimate input. If it really is that long for a four-digit number, maybe we go with your idea; that's pretty bad.
comment:10 Changed 12 years ago by
Well, it hasn't stopped after half an hour..
I don't quite follow you about the speed, though. Before every "return ans" for n > 5 you either (1) test primality, (2) test for a prime power, or (3) once again test for a prime power, so ISTM you're not avoiding the cost of doing the tests on legitimate input, you're just doing it after the znprimroot call instead of before. (And we get very similar times for large legitimate input, which is what I'd expect if I'm right.)
If we really want a super-fast way to do this when needed -- and like you I suspect, I have zero need for primitive roots of huge numbers so I have no idea if this matters at all in practice :^)
-- I say we add a proof=False option, or whatever the standard is in other functions that don't test their input for speed reasons (I think I've seen exactly this situation before but can't think where right now).
comment:11 Changed 12 years ago by
Status: | needs_review → needs_work |
---|
Yes, you are right about this. I forgot that the number isn't actually returned until then :)
Yeah, I thought about proof=False, but wanted to keep it short to start - adding keywords sometimes turns into a mess.
If you want to post a second patch with these ideas, we can review it jointly very easily. To me, it's half a dozen of one and six of the other, although I am surprised myself about
sage: a = pari(1729) sage: a.znprimroot()
taking so long.
comment:12 Changed 12 years ago by
Wow, regarding
sage: a = pari(1729) sage: a.znprimroot()
That's seriously messed up! This is a case where there is no primitive root. So it's the "undefined" behavior mentioned in the pari docs -- just take forever.
Anyway, looking over this thread, it seems to me that the docstring in the patch is very good as is, but the code should be completely changed to only call pari after explicitly checking the condition for there to be a primitive root (as explained in the docstring).
comment:13 Changed 12 years ago by
Status: | needs_work → needs_review |
---|
comment:14 follow-up: 17 Changed 12 years ago by
Authors: | Karl-Dieter Crisman → Karl-Dieter Crisman, William Stein |
---|---|
Reviewers: | → William Stein, Karl-Dieter Crisman |
Status: | needs_review → needs_info |
Nice use of is_prime_power
to make sure negatives don't get through.
Should we throw a better error message than
AttributeError: 'sage.rings.rational.Rational' object has no attribute 'is_prime_power'
for non-integer input? Or is it enough to say that there isn't a primitive root? I guess same question for negatives. "No primitive root" seems somewhat cryptic. Putting 'needs info'; I can also upload a reviewer patch if you think that's okay. We would check for that after we actually got the primitive root, so it wouldn't slow anything down.
Otherwise the patch looks good, and if William thinks check=True
is good enough, it's fine for me.
comment:15 Changed 12 years ago by
I'm also getting something weird while trying to test one of the files where this occurs.
sage -tp 4 devel/sage/sage/rings/finite_rings Testing that Sage starts... Yes, Sage starts. Global iterations: 1 File iterations: 1 Using cached timings to run longest doctests first. Doctesting 0 files doing 0 jobs in parallel ---------------------------------------------------------------------- All tests passed! Timings have been updated. Total time for all tests: 0.2 seconds
But
ls Downloads/sage-4.7.alpha1/devel/sage/sage/rings/finite_rings/ __init__.py finite_field_base.c all.py finite_field_base.pxd constructor.py finite_field_base.pyx element_base.c finite_field_ext_pari.py element_base.pxd finite_field_givaro.py element_base.pyx finite_field_prime_modn.py element_ext_pari.py homset.py element_givaro.cpp integer_mod.c element_givaro.pxd integer_mod.pxd element_givaro.pyx integer_mod.pyx element_ntl_gf2e.cpp integer_mod_ring.py element_ntl_gf2e.pxd stdint.h element_ntl_gf2e.pyx
What is up with that?
comment:16 Changed 12 years ago by
Okay, that turned out to be a red herring. At least I think so. Sorry for the noise.
Relevant tests pass. I don't have time to make a reviewer patch now, but would appreciate the feedback on whether it's okay to add that check and format. Since it would be after the real work is done, I think that is ok.
comment:17 Changed 12 years ago by
Status: | needs_info → needs_review |
---|
Should we throw a better error message than
AttributeError: 'sage.rings.rational.Rational' object has no attribute 'is_prime_power'for non-integer input?
Actually, we end up with
TypeError: Attempt to coerce non-integral RealNumber to Integer
which I didn't notice but is obvious because of the ZZ(n)
. So that is okay.
Assuming the docs look ok, I now think that adding the other thing is not worth the effort. I think I may need to give a reviewer patch to fix a couple REsT things, though. We'll see :)
comment:18 Changed 12 years ago by
Description: | modified (diff) |
---|---|
Status: | needs_review → positive_review |
Had to fix one tiny thing for the note to look right. Positive review.
Apply only trac_10836.2.patch.
comment:19 Changed 12 years ago by
Authors: | Karl-Dieter Crisman, William Stein → Karl-Dieter Crisman, William Stein, Jeroen Demeyer |
---|---|
Description: | modified (diff) |
Status: | positive_review → needs_work |
comment:20 Changed 12 years ago by
Status: | needs_work → needs_review |
---|
My additional patch needs review.
Changed 12 years ago by
Attachment: | 10836_additional.patch added |
---|
comment:21 follow-up: 22 Changed 12 years ago by
Whoah, Nelly! I see the 'niceification'. But two points.
- PARI doc explicitly says that it does not give the smallest root for non-prime input.
- Do negative integers have primitive roots? I have never heard this. I understand the extension of the definition, but a generator of "Z/(-5)Z" somehow seems different to me than a "primitive root of -5". I was able to crash !W|A with this, incidentally :)
comment:22 Changed 12 years ago by
Replying to kcrisman:
Do negative integers have primitive roots? I have never heard this. I understand the extension of the definition, but a generator of "Z/(-5)Z" somehow seems different to me than a "primitive root of -5". I was able to crash !W|A with this, incidentally :)
Well, I think we should decide what we do with negative integers: either raise an exception or return a sensible value. It seemed artificial to me not to allow negative numbers. Also PARI/GP computes primitive roots of negative numbers.
comment:23 follow-up: 24 Changed 12 years ago by
Status: | needs_review → needs_work |
---|
Well, what we say does refer to just multipl. gp. mod n, so ok.
But should definitely make sure not to say we give smallest prim root. So 'needs work' to fix that. But this is pretty bad not to have fixed, so let's try to get it in asap.
comment:24 follow-up: 25 Changed 12 years ago by
Replying to kcrisman:
But should definitely make sure not to say we give smallest prim root.
Why not? If you things like "this is pretty bad not to have fixed", at least say why.
comment:25 follow-up: 26 Changed 12 years ago by
Replying to jdemeyer:
Replying to kcrisman:
But should definitely make sure not to say we give smallest prim root.
Why not? If you things like "this is pretty bad not to have fixed", at least say why.
Because the PARI documentation says that it no longer gives the smallest prim root, except for primes. Or at least the changeset I saw a while ago says this - please don't ask me to find it. Presumably something about the internal algorithm changed?
Oh, I see what you mean now! "Pretty bad to not have fixed" meant the actual ticket. Wrong mathematics! Not the tiny change I am mentioning above. I'm sorry for the possible misinterpretation.
comment:26 Changed 12 years ago by
Replying to kcrisman:
Because the PARI documentation says that it no longer gives the smallest prim root, except for primes.
But this is consistent with my patch. In this patch, I write "If n
is prime, this is the smallest primitive root.". You mean I should be more explicit about the non-prime case?
comment:27 follow-up: 28 Changed 12 years ago by
Reviewers: | William Stein, Karl-Dieter Crisman → William Stein, Karl-Dieter Crisman, Jeroen Demeyer |
---|---|
Status: | needs_work → positive_review |
My humble apologies. I don't know why I thought it said
"A primitive root of n
. This is the smallest primitive root."
That means we should, in theory, fix
"k
is a nonnegative number."
and
"assumed to be a positive integer possessing a primitive root"
But I agree that unless it's very easy to update the patch for this, this should go in. So I'm inclined to put positive review. Docs look right, still passes tests.
comment:28 Changed 12 years ago by
Replying to kcrisman:
That means we should, in theory, fix
Actually, I think both of these statements are 100% correct in the patch.
"
k
is a nonnegative number."
Here, we are talking about exponents which are allowed to be zero, so this is correct.
"assumed to be a positive integer possessing a primitive root"
This is also correct because PARI's behaviour for negative arguments to znprimroot()
is undocumented, so we cannot rely on it (in practice, it does work though).
comment:29 Changed 12 years ago by
I clearly need more sleep. You're right on the second count, and I copied the wrong thing on the first count:
"A primitive root exists if n=4
or n=p^k
or n=2p^k
, where p
is an odd prime"
Which in theory then should say something about \pm
all those numbers, given that this patch defines it that way. Probably that is lawyering. But wanted to put it out there. Anyway, doesn't really bother me.
comment:30 Changed 12 years ago by
Merged in: | → sage-4.7.alpha4 |
---|---|
Resolution: | → fixed |
Status: | positive_review → closed |
From one of the Pari developers: