primitive root is broken

[kcrisman]

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 attachment: trac_10836.2.patch and attachment: 10836_additional.patch.

Upstream: Reported upstream. Developers deny it's a bug.

CC: @williamstein @adeines

Component: number theory

Author: Karl-Dieter Crisman, William Stein, Jeroen Demeyer

Reviewer: William Stein, Karl-Dieter Crisman, Jeroen Demeyer

Merged: sage-4.7.alpha4

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

[kcrisman]

comment:1

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.

[kcrisman]

Changed upstream from Not yet reported upstream; Will do shortly. to Reported upstream. Developers deny it's a bug.

[kcrisman]

comment:2

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

[56048686-9665-4c5e-973e-6c3add3aa805]

comment:3

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.

[kcrisman]

comment:4

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?

[kcrisman]

Author: Karl-Dieter Crisman

[kcrisman]

comment:5

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.

[kcrisman]

Attachment: trac_10836-primitive.patch.gz

[kcrisman]

comment:6

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.

[kcrisman]

comment:7

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

[56048686-9665-4c5e-973e-6c3add3aa805]

comment:8

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.

[kcrisman]

comment:9

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.

[56048686-9665-4c5e-973e-6c3add3aa805]

comment:10

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

[kcrisman]

comment:11

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.

[williamstein]

comment:12

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

[williamstein]

Attachment: trac_10836.patch.gz

apply only this patch

[kcrisman]

Changed author from Karl-Dieter Crisman to Karl-Dieter Crisman, William Stein

[kcrisman]

comment:14

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.

[kcrisman]

Reviewer: William Stein, Karl-Dieter Crisman

[kcrisman]

comment:15

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?

[kcrisman]

comment:16

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.

[kcrisman]

comment:17

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

[kcrisman]

Apply only this patch

[kcrisman]

Description changed:

--- 
+++ 
@@ -25,3 +25,6 @@
 <snip> 
 ArithmeticError: There is no primitive root modulo n 

+ +--- +Apply only attachment: trac_10836.2.patch.

[kcrisman]

comment:18

Attachment: trac_10836.2.patch.gz

Had to fix one tiny thing for the note to look right. Positive review.

Apply only attachment: trac_10836.2.patch.

[jdemeyer]

Changed author from Karl-Dieter Crisman, William Stein to Karl-Dieter Crisman, William Stein, Jeroen Demeyer

[jdemeyer]

Description changed:

--- 
+++ 
@@ -27,4 +27,4 @@

---

-Apply only attachment: trac_10836.2.patch. +Apply attachment: trac_10836.2.patch and attachment: 10836_additional.patch.

[jdemeyer]

comment:20

My additional patch needs review.

[jdemeyer]

Attachment: 10836_additional.patch.gz

[kcrisman]

comment:21

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

[jdemeyer]

comment:22

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.

[kcrisman]

comment:23

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.

[jdemeyer]

comment:24

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.

[kcrisman]

comment:25

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.

[jdemeyer]

comment:26

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?

[kcrisman]

Changed reviewer from William Stein, Karl-Dieter Crisman to William Stein, Karl-Dieter Crisman, Jeroen Demeyer

[kcrisman]

comment:27

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.

[jdemeyer]

comment:28

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

[kcrisman]

comment:29

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.

[jdemeyer]

Merged: sage-4.7.alpha4