Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

faster is_prime #16878

Closed
videlec opened this issue Aug 25, 2014 · 88 comments
Closed

faster is_prime #16878

videlec opened this issue Aug 25, 2014 · 88 comments

Comments

@videlec
Copy link
Contributor

videlec commented Aug 25, 2014

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.

Depends on #16997

CC: @nathanncohen

Component: number theory

Author: Vincent Delecroix, Jeroen Demeyer

Branch: 1f8abd9

Reviewer: Jeroen Demeyer, Vincent Delecroix

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

@videlec videlec added this to the sage-6.4 milestone Aug 25, 2014
@videlec
Copy link
Contributor Author

videlec commented Aug 25, 2014

New commits:

dd44a7etrac #16878: faster is_prime

@videlec
Copy link
Contributor Author

videlec commented Aug 25, 2014

Branch: u/vdelecroix/16878

@videlec
Copy link
Contributor Author

videlec commented Aug 25, 2014

Commit: dd44a7e

@videlec

This comment has been minimized.

@jdemeyer
Copy link
Contributor

comment:3

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
Copy link
Contributor

Reviewer: Jeroen Demeyer

@jdemeyer
Copy link
Contributor

comment:4

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
Copy link
Contributor

comment:5

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

cpdef _pari_(self)

accomplish the same thing?

@videlec
Copy link
Contributor Author

videlec commented Aug 25, 2014

comment:6

Thanks for the remark. One question: are we sure that 264 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

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Aug 25, 2014

Changed commit from dd44a7e to d7dc389

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Aug 25, 2014

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

37fc8e8Upgrade to PARI-2.7.1
5db54b6Trac 15767: reviewer patch
1d103daTrac 15767: fix doctest in Ser()
62ca821Trac 15767: explain application of Sturm's theorem
bf435f8Merge tag '6.4.beta1' into ticket/15767
6e9f686trac #16878: faster is_prime
d7dc389trac #16878: review

@videlec
Copy link
Contributor Author

videlec commented Aug 25, 2014

Dependencies: #15767

@jdemeyer
Copy link
Contributor

comment:9

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
Copy link
Contributor

comment:10

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

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Aug 25, 2014

Changed commit from d7dc389 to 175d6b5

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Aug 25, 2014

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

175d6b5trac #16878: review 2

@jdemeyer
Copy link
Contributor

comment:13

This is bikeshedding I know, but wouldn't

mpz_ui_pow_ui(PARI_PSEUDOPRIME_LIMIT, 2, 64)

be easier to understand?

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Feb 26, 2015

Changed commit from 96f1f50 to 1f8abd9

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Feb 26, 2015

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

1f8abd9Doctest fix

@jdemeyer
Copy link
Contributor

comment:60

Now all doctests pass.

@videlec
Copy link
Contributor Author

videlec commented Feb 26, 2015

comment:61
  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
Copy link
Contributor

comment:62

Replying to @videlec:

  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

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

@videlec
Copy link
Contributor Author

videlec commented Feb 26, 2015

comment:63

Replying to @jdemeyer:

Replying to @videlec:

  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.

  1. 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
Copy link
Contributor

comment:64

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.

@jdemeyer

This comment has been minimized.

@videlec
Copy link
Contributor Author

videlec commented Feb 27, 2015

comment:67

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
Copy link
Contributor

comment:68

Replying to @videlec:

The buildbot is complaining about a digit

Fixed by #16997.

@videlec
Copy link
Contributor Author

videlec commented Mar 4, 2015

comment:69

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

Vincent

@vbraun
Copy link
Member

vbraun commented Mar 5, 2015

Changed branch from u/jdemeyer/16878 to 1f8abd9

@fchapoton
Copy link
Contributor

Changed commit from 1f8abd9 to none

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

4 participants