faster is_prime

[vdelecroix]

Right now to test if a Sage integer is prime it is faster to call prime_range rather than .is_prime()...

sage: timeit("bool(prime_range(121,122))", number=10000)
10000 loops, best of 3: 1.09 µs per loop
sage:  timeit("bool(prime_range(1009,1010))", number=10000)
10000 loops, best of 3: 1.2 µs per loop

versus

sage: timeit("121.is_prime()", number=10000)
10000 loops, best of 3: 5.3 µs per loop
sage: timeit("1009.is_prime()", number=10000)
10000 loops, best of 3: 4.19 µs per loop

The patch does some tiny modifications in integer.pyx and we get

sage: timeit("121.is_prime()", number=10000)
10000 loops, best of 3: 812 ns per loop
sage: timeit("1009.is_prime()", number=10000)
10000 loops, best of 3: 730 ns per loop

We also modify is_prime_power() to return False for 1 and raise an error for non-integral rationals like 1/2.

See also ​this sage-devel discussion.

[vdelecroix]

New commits:

​dd44a7e

trac #16878: faster is_prime

[jdemeyer]

I would prefer not to include stuff from pari/pari.h directly, but instead use the declarations from src/sage/libs/pari/decl.pxi. And please rebase this over #15767.

[jdemeyer]

Can you also change the value of PARI_PSEUDOPRIME_LIMIT (see PARI docs) to 2^64 and use mpz_fits_ui() before calling mpz_get_ui()

[jdemeyer]

Instead of changing _pari_() to _pari_c() everywhere, wouldn't

cpdef _pari_(self)

accomplish the same thing?

[vdelecroix]

Thanks for the remark. One question: are we sure that 2⁶⁴ fits into an unsigned long on any platform? Otherwise we can not safely call uisprime and uisprimepower. If not I will use a double check

if mpz_cmp(PARI_PSEUDOPRIME_LIMIT, v) > 0 and mpz_cmp_ui(v, ULONG_MAX) < 0:
     # use uisprime or uisprimepower

Vincent

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​37fc8e8	Upgrade to PARI-2.7.1
​5db54b6	Trac 15767: reviewer patch
​1d103da	Trac 15767: fix doctest in Ser()
​62ca821	Trac 15767: explain application of Sturm's theorem
​bf435f8	Merge tag '6.4.beta1' into ticket/15767
​6e9f686	trac #16878: faster is_prime
​d7dc389	trac #16878: review

[jdemeyer]

10000000000000000 != 2^64 :-)

I don't think we should lower PARI_PSEUDOPRIME_LIMIT because it's also used in places where it's not bounded by ULONG_MAX, see next_prime() for example.

Personally, I would assume that ULONG_MAX < 2^64 and just use mpz_fits_uint_p().

[jdemeyer]

Also, you wrote "set set" instead of "set".

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​175d6b5

trac #16878: review 2

[jdemeyer]

This is bikeshedding I know, but wouldn't

mpz_ui_pow_ui(PARI_PSEUDOPRIME_LIMIT, 2, 64)

be easier to understand?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​d55dc67

trac #16878: bikeshedding

[vdelecroix]

I don't mind...

---

New commits:

​d55dc67

trac #16878: bikeshedding

[vdelecroix]

ping! (applies cleanly on 6.4.beta2)

[jdemeyer]

Big reviewer patch coming up...

[jdemeyer]

As you can see, I made quite a lot of additional changes, mostly trying to make the code cleaner. I am currently running doctests.

Can you review this extra commit?

---

New commits:

​7f61321

Reorganize is_prime_power() and next_prime() logic

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​5937af5

Sloane considers 1 to be a prime power

[vdelecroix]

Great! I especially like that now isprimepower is accessible.

In my last commit, I added the keyword certificate to is_prime_power in order to get back the pair (n,p) such that the integer is pⁿ (following pari convention). I think it is worth it!

Vincent

---

New commits:

​f0d716f

trac #16878: certificate flag in is_prime_power

[vdelecroix]

Note the treatment of 1 is different inis_prime_powers and prime_powers... either we fix it in that ticket or it can wait #16880.

Vincent

[jdemeyer]

Replying to vdelecroix:

Great! I especially like that now isprimepower is accessible.

In my last commit, I added the keyword certificate to is_prime_power in order to get back the pair (n,p) such that the integer is pⁿ (following pari convention). I think it is worth it!

I like the idea, but I'm not sure that I like the interface with the certificate keyword.

How about a separate method write_as_prime_power() (raising ArithmeticError if it's not a prime power) or something?

[jdemeyer]

Replying to vdelecroix:

Note the treatment of 1 is different inis_prime_powers and prime_powers... either we fix it in that ticket or it can wait #16880.

Let's ignore that for now and try to merge this ticket first.

[vdelecroix]

Replying to jdemeyer:

Replying to vdelecroix:

Great! I especially like that now isprimepower is accessible.

In my last commit, I added the keyword certificate to is_prime_power in order to get back the pair (n,p) such that the integer is pⁿ (following pari convention). I think it is worth it!

I like the idea, but I'm not sure that I like the interface with the certificate keyword.

How about a separate method write_as_prime_power() (raising ArithmeticError if it's not a prime power) or something?

Note that there is already factor that does the job. Moreover, having two methods (is_prime_power and write_as_prime_power) with the exact same code looks bad as well. From my implementation, it is clear that the certificate keyword only modifies the output.

My design choice also comes from the fact that this certifcate keyword is already in use (see sage.graphs or sage.combinat.designs). I am not particularly happy with it but I find it better than having two functions.

On the other hand, treating a Python exception takes time, and this is_prime_power is the kind of function that I might use in time critical code.

Vincent

[jdemeyer]

One more reviewer patch...

---

New commits:

​bfd69d3

is_prime_power() reviewer fixes

[vdelecroix]

I am happy with the current version (at commit bfd69d3). Do you agree to set to positive review?

Vincent

[vbraun]

This takes about a second before this ticket and really long after this ticket:

sage: with proof.WithProof('arithmetic', False):  
....:    sage.rings.arith.is_prime_power((10^1000 + 453)^2)

In fact, it get so slow that the ARM buildbot times out in

sage -t --long src/sage/rings/finite_rings/constructor.py
    Timed out
**********************************************************************
...
sage: k = FiniteField(10^1000 + 453, proof=False) ## line 230 ##
sage: k = FiniteField((10^1000 + 453)^2, 'a', proof=False)      # long time -- about 5 seconds ## line 231 ##

**********************************************************************

[vdelecroix]

very strange indeed

sage: sage.rings.arith.is_prime_power((10^999 + 453)^2)   # very fast
False
sage: sage.rings.arith.is_prime_power((10^1001 + 453)^2)  # very fast
False
sage: sage.rings.arith.is_prime_power((10^1000 + 453)^2)
... <still waiting> ...

[vdelecroix]

I see... seems to be because we do not care about the flag proof=True or proof=False anymore in is_prime_power. We should still rely on is_pseudoprime_small_power when proof=False.

And by the way, in is_pseudoprime_small_power the choice of the keyword to return (n,p) is get_data. I guess we should change the certificate in is_prime_power to get_data as well.

Vincent

[vdelecroix]

New commits:

​6e9f686	trac #16878: faster is_prime
​d7dc389	trac #16878: review
​175d6b5	trac #16878: review 2
​d55dc67	trac #16878: bikeshedding
​7f61321	Reorganize is_prime_power() and next_prime() logic
​5937af5	Sloane considers 1 to be a prime power
​f0d716f	trac #16878: certificate flag in is_prime_power
​bfd69d3	is_prime_power() reviewer fixes
​e6f88c5	trac #16878: test is_pseudoprime_small_power + clean

[jdemeyer]

In is_pseudoprime_small_power, the fixed exponent bound of 1024 is clearly not acceptable.

The PARI implementation of isprimepower() is as follows: first do trial division by small primes up to 100. If we do not find a divisor, then find out if the number is a perfect power. This is now a lot cheaper since the exponent can be at most log(n)/log(p), where n is the number to check and p the smallest prime that we didn't check in the trial division.

But overall, I would prefer PARI to add a new function ispseudoprimepower() such that Sage can interface it.

[vdelecroix]

Replying to jdemeyer:

In is_pseudoprime_small_power, the fixed exponent bound of 1024 is clearly not acceptable.

The PARI implementation of isprimepower() is as follows: first do trial division by small primes up to 100. If we do not find a divisor, then find out if the number is a perfect power. This is now a lot cheaper since the exponent can be at most log(n)/log(p), where n is the number to check and p the smallest prime that we didn't check in the trial division.

But overall, I would prefer PARI to add a new function ispseudoprimepower() such that Sage can interface it.

Me too. I saw that you already submit the question ​http://pari.math.u-bordeaux.fr/archives/pari-dev-1409/msg00011.html...

On the other hand, we need a reasonable fix in order to get this ticket in.

Vincent

[vbraun]

Jeroen, did you end up writing a pari patch? Maybe we can use the current is_pseudoprime_small_power until there is a better implementation?

[jdemeyer]

Replying to vbraun:

Jeroen, did you end up writing a pari patch?

Not yet, but it would be trivial to do that. I will add it to #16997.

[jdemeyer]

On the other hand, perhaps better add the patch here...

[jdemeyer]

...or a new ticket.

[jdemeyer]

Turns out adding the PARI patch to #16997 was the easiest (since there would have been merge conflicts otherwise). It's now done.

[jdemeyer]

I have now a version of #16997 which is ready for review. For this ticket, I would wait until that ticket is settled and then have a look again at this.

[vdelecroix]

Hi,

Good news: Karim recently committed a ispseudoprimepower in pari (​commit 0f20474).

Vincent

[jdemeyer]

New commits:

​ade4f75	Upgrade PARI to git master
​9751ee8	Trac #16878: faster is_prime

[vdelecroix]

(I guess you did it because of conflicts with #17852)

[jdemeyer]

I'm doing various fixes, like remove is_pseudoprime_small_power.

I think it would be better to rebase #17852 on this ticket, since this ticket removes lots of things that #17852 fixes.

[vdelecroix]

Replying to jdemeyer:

I'm doing various fixes, like remove is_pseudoprime_small_power.

I think it would be better to rebase #17852 on this ticket, since this ticket removes lots of things that #17852 fixes.

Fine for me.

[jdemeyer]

I would also like to switch around the get_data result: when given p^n as input, I would find it more natural to return (p, n) instead of (n, p). What do you think?

[vdelecroix]

gp is doing the other way around: it returns k as the main argument since it is compatible with the boolean answer. But doing (p,n) will be consistent with Integer.perfect_power so I would go for it.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​96f1f50

Trac #16878: further fixes to is_prime() and friends

[jdemeyer]

As far as I'm concerned, this is the final version. If you agree, then you can set this to positive_review.

I haven't yet run all doctests, doing that right now...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1f8abd9

Doctest fix

[jdemeyer]

Now all doctests pass.

[vdelecroix]

1. Now arith.is_prime_power and arith.is_pseudoprime_power do not handle negative powers anymore (i.e. rational numbers p^-n). This is different from previous versions.

2. In Integer.is_prime_power why did you use a long for n and p instead of an unsigned long. Is pari uisprimepower not able to handle up to ULONG_MAX?

[jdemeyer]

Replying to vdelecroix:

1. Now arith.is_prime_power and arith.is_pseudoprime_power do not handle negative powers anymore (i.e. rational numbers p^-n). This is different from previous versions.

Yes, this was discussed on sage-devel and the former behaviour was considered to be a bug => no deprecation needed

2. In Integer.is_prime_power why did you use a long for n and p instead of an unsigned long. Is pari uisprimepower not able to handle up to ULONG_MAX?

I assume that pari can handle it, but smallInteger() for the result takes a long.

[vdelecroix]

Replying to jdemeyer:

Replying to vdelecroix:

1. Now arith.is_prime_power and arith.is_pseudoprime_power do not handle negative powers anymore (i.e. rational numbers p^-n). This is different from previous versions.

Yes, this was discussed on sage-devel and the former behaviour was considered to be a bug => no deprecation needed

True, I forgot about here and it is ​this thread.

2. In Integer.is_prime_power why did you use a long for n and p instead of an unsigned long. Is pari uisprimepower not able to handle up to ULONG_MAX?

I assume that pari can handle it, but smallInteger() for the result takes a long.

I see... I do not like the fact that we avoid pari uisprimepower for the range between LONG_MAX and ULONG_MAX. Especially because otherwise we import the proof module which takes some non negligible amount of time. We can either get rid of smallInteger (that would be bad if either p or n would be below 256) or add a smallunsignedInteger

cdef inline Integer smallunsignedInteger(unsigned long value):
    """
    This is the fastest way to create a (likely) small positive Integer.
    """
    cdef Integer z
    if value <= small_pool_max:
        return <Integer>small_pool[value - small_pool_min]
    else:
        z = PY_NEW(Integer)
        mpz_set_ui(z.value, value)
        return z

What do you think?

[jdemeyer]

To be honest, I think introducing smallunsignedInteger is just overkill. I don't think it makes a big difference if the limit for a fast is_prime_power is 2^63 or 2^64.

[vdelecroix]

The buildbot is complaining about a digit

sage -t --long --warn-long 35.8 src/sage/libs/pari/gen.pyx
**********************************************************************
File "src/sage/libs/pari/gen.pyx", line 6237, in sage.libs.pari.gen.gen.ellheight
Failed example:
    e.ellheight([1,0], precision=128).sage()
Expected:
    0.47671165934373953737948605888465305945902294217            
Got:
    0.47671165934373953737948605888465305945902294218
**********************************************************************

Might be related to the doctest failure in pari update in #16997.

[jdemeyer]

Replying to vdelecroix:

The buildbot is complaining about a digit

Fixed by #16997.

[vdelecroix]

Merges cleanly with #16997 and all tests pass! Cool.

Vincent