Opened 12 years ago

Closed 12 years ago

# primitive root is broken

Reported by: Owned by: Karl-Dieter Crisman William Stein critical sage-4.7 number theory William Stein, Alyson Deines sage-4.7.alpha4 Karl-Dieter Crisman, William Stein, Jeroen Demeyer William Stein, Karl-Dieter Crisman, Jeroen Demeyer Reported upstream. Developers deny it's a bug.

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.        |
----------------------------------------------------------------------
sage: primitive_root(15)
<snip>
ArithmeticError: There is no primitive root modulo n


### comment:1 Changed 12 years ago by Karl-Dieter Crisman

Cc: William Stein added Not yet reported upstream; Will do shortly. → Reported upstream. Developers deny it's a bug.

From one of the Pari developers:

Please read the documentation of znprimroot in PARI 2.4.3:

? ??znprimroot
znprimroot(n):

returns a primitive root (generator) of (Z/nZ)^*, whenever this latter group is cyclic (n = 4
or n = 2p^k or n = p^k,  where p is an odd prime and k >= 0).   If the group is not cyclic,  the
result  is undefined.   If n is a prime,  then the smallest positive primitive root is returned.
This is no longer true for composites.

(Z/15Z)* is not cyclic so znprimroot(15) is undefined.


### comment:2 Changed 12 years ago by Karl-Dieter Crisman

Ergo, we will have to catch such cases in an effort to remain mathematically correct.

### comment:3 Changed 12 years ago by D.S. McNeil

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 Karl-Dieter Crisman

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 Karl-Dieter Crisman

Authors: → Karl-Dieter Crisman → sage-4.7 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.

### comment:6 Changed 12 years ago by Karl-Dieter Crisman

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 Karl-Dieter Crisman

Hopefully someone can review this during Sage Day 29? It's mathematically wrong currently...

### comment:8 Changed 12 years ago by D.S. McNeil

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 Karl-Dieter Crisman

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 D.S. McNeil

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 Karl-Dieter Crisman

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 William Stein

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).

### Changed 12 years ago by William Stein

apply only this patch

### comment:13 Changed 12 years ago by William Stein

Status: needs_work → needs_review

### comment:14 follow-up:  17 Changed 12 years ago by Karl-Dieter Crisman

Authors: Karl-Dieter Crisman → Karl-Dieter Crisman, William Stein → William Stein, Karl-Dieter Crisman 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 Karl-Dieter Crisman

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 Karl-Dieter Crisman

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 in reply to:  14 Changed 12 years ago by Karl-Dieter Crisman

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 :)

### Changed 12 years ago by Karl-Dieter Crisman

Apply only this patch

### comment:18 Changed 12 years ago by Karl-Dieter Crisman

Description: modified (diff) 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 Jeroen Demeyer

Authors: Karl-Dieter Crisman, William Stein → Karl-Dieter Crisman, William Stein, Jeroen Demeyer modified (diff) positive_review → needs_work

### comment:20 Changed 12 years ago by Jeroen Demeyer

Status: needs_work → needs_review

### comment:21 follow-up:  22 Changed 12 years ago by Karl-Dieter Crisman

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 in reply to:  21 Changed 12 years ago by Jeroen Demeyer

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 Karl-Dieter Crisman

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 in reply to:  23 ; follow-up:  25 Changed 12 years ago by Jeroen Demeyer

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 in reply to:  24 ; follow-up:  26 Changed 12 years ago by Karl-Dieter Crisman

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 in reply to:  25 Changed 12 years ago by Jeroen Demeyer

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 Karl-Dieter Crisman

Reviewers: William Stein, Karl-Dieter Crisman → William Stein, Karl-Dieter Crisman, Jeroen Demeyer 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 in reply to:  27 Changed 12 years ago by Jeroen Demeyer

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 Karl-Dieter Crisman

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 Jeroen Demeyer

Merged in: → sage-4.7.alpha4 → fixed positive_review → closed
Note: See TracTickets for help on using tickets.