Re: [Python-ideas] Add the imath module
On 2018-07-12 20:11, Steven D'Aprano wrote:
about 4.7 seconds to test 2**800 + 1;
In SageMath: sage: n = 2**800+1; timeit('is_prime(n)') 625 loops, best of 3: 303 µs per loop That's 4 orders of magnitude faster...
On Fri, Jul 13, 2018 at 10:10:00AM +0200, Jeroen Demeyer wrote:
Go on, rub it in why don't you... I'm pretty sure that Sage is using numpy, which would almost certainly have a C-based implementation, not pure Python like mine. And it's probably making far less conservative choices about the rate of false positives it is willing to accept. And your computer is probably newer and faster than mine :-( My point was that for casual use, even the bad cases are still fast enough for interactive use. If you think 5 seconds is ridiculously slow, you've never tried a naive trial division algorithm that took literally days to crawl through the millions of divisions needed. *wink* -- Steve
[Steven D'Aprano] about 4.7 seconds to test 2**800 + 1;
[Jeroen Demeyer]
More like 6, yes? My Python version is more like 3, but is way fast enough for most things I do: $ py -m timeit -s "from temp3 import pp; n=2**800+1" "pp(n)" 200 loops, best of 5: 1.71 msec per loop $ py -m timeit -s "from temp3 import pp; n=2**900+1" "pp(n)" 100 loops, best of 5: 2.32 msec per loop $ py -m timeit -s "from temp3 import pp; n=2**1000+1" "pp(n)" 100 loops, best of 5: 3.15 msec per loop Steven's numbers are pretty baffling to me, since these are all composite and so iterating Miller-Rabin "should get out" pretty fast: about 4.7 seconds to test 2**800 + 1;
less than a tenth of a millisecond to test 2**900 + 1; and about 8.6 seconds to test 2**1000 + 1.
Here's the Python code I used: def pp(n, k=10): """ Return True iff n is a probable prime, using k Miller-Rabin tests. When False, n is definitely composite. """ from random import randrange if n < 4: return 2 <= n <= 3 d = nm1 = n-1 s = 0 while d & 1 == 0: s += 1 d >>= 1 if s == 0: return False # n >= 4 and even ss = range(1, s) for _ in range(k): a = randrange(2, nm1) x = pow(a, d, n) if x == 1 or x == nm1: continue for _ in ss: x = x*x % n if x == nm1: break if x == 1: return False else: return False return True
If we have a fast and exact primality test up to 2**64, a signature like this would help: def is_prime(n, probabilistic=False): When probabilistic is False and the number is beyond 2**64, it'll return OverflowError. When probabilistic is True, it's unbounded, but probabilistic. The docs should be explicit about that, though. Instead of an OverflowError, a number beyond 2**64 might simply use a slower strategy. But to make it an explicit option, the signature could be: def is_prime(n, exact=True, unbounded=False): By tellling it should be both exact and unbounded, the user is aware that it might be slow. The alternatives are: - exact, unbounded: Slower test - not exact, unbounded: Probabilistic test - exact, bounded: Default, raises OverflowError beyond 2**64 - not exact, bounded: Invalid, but it can just ignore exact input -- Danilo J. S. Bellini --------------- "*It is not our business to set up prohibitions, but to arrive at conventions.*" (R. Carnap)
[Tim]
More like 6, yes?
Heh - sorry about that. A speck of dirt on my monitor made me read '303" as "3.03".
My Python version is more like 3, but is way fast enough for most things I do:
So my Python Miller-Rabin code (already shared) was about one order of magnitude slower. ...
On Fri, Jul 13, 2018 at 10:52:38AM -0500, Tim Peters wrote:
Steven's numbers are pretty baffling to me, since these are all composite and so iterating Miller-Rabin "should get out" pretty fast:
That's because you haven't seen my code :-) It's over-engineered, class-based, and written as a learning exercise. There's a compatibility layer so it will work back to Python 2.4 and an instrumentation layer which I inserted in a fit of enthusiasm to gather staticistics, which I have never once looked at since. And *even so* it is still fast enough for casual use at the interactive interpreter, compared to more naive algorithms with worse Big O performance characteristics. You might think 5 seconds is slow, but I'm serious when I say some of the other algorithms I played with would take *days* to generate, or check, largish primes. -- Steve
Looking for some information on how long it has taken to generate large primes I stumbled across https://arxiv.org/ftp/arxiv/papers/1709/1709.09963.pdf which makes interesting reading. It concentrates on giving no false negatives (never saying n is not a prime when it is) but giving an increasing probability that n is a prime as testing goes on. I did find an article from 2016 that mentioned M77232917 (a 23,249,425 digit Mersenne prime) and M74207281 (22,338,618 digit Mersenne prime) with the latter taking a month to check with the Lucas-Lehmer test. On 14/07/2018 05:54, Steven D'Aprano wrote:
-- Steve (Gadget) Barnes Any opinions in this message are my personal opinions and do not reflect those of my employer. --- This email has been checked for viruses by AVG. https://www.avg.com
[Steve Barnes]
That's a very nice, gentle introduction! Thanks for sharing it. The `pp()` function I posted before is a Python implementation of the Miller-Rabin test they ended with. The latter is very widely used, because the code is relatively simple. short, and fast, and has no "bad cases". For _generating_ primes, the article uses other tricks to weed out sure-to-be-composite candidates before bothering with probable-prime testing. But the tricks they use in the paper are specific to a base 10 representation. In binary, it's common to skip candidates such that gcd(candidate, 2 * 3 * 5 * ...) > 1 where the product-of-small-primes second argument is picked to be the largest such that bigint division by it is exceptionally fast (usually one "digit" in whatever internal power-of-2 base is used by the bigint implementation). I did find an article from 2016 that mentioned M77232917 (a 23,249,425
digit Mersenne prime) and M74207281 (22,338,618 digit Mersenne prime) with the latter taking a month to check with the Lucas-Lehmer test.
More, it's a special version of Lucas-Lehmer (LL) that only works for testing candidates of the form 2**p-1 where p is itself an odd prime. In that context, it's not probabilistic - it definitively answers whether 2**p-1 is prime. Fast as it is, Miller-Rabin is impractical for guessing about numbers of this size. And as enormously faster as _it_ is, the special version of LL used for Mersenne testing is largely hand-coded in assembler to squeeze out every possible machine cycle. I don't want any of that in Python's standard library either ;-) Here! You should enjoy this account of the even larger Mersenne prime discovered this year: https://www.mersenne.org/primes/press/M77232917.html Note that the hyper-optimized testing code they use is so historically error-prone (& so good at provoking HW errors on overheated machines!) that they now verify each new Mersenne prime by 4 entirely different LL implementations, on 4 different machines.
[Tim]
Steven's numbers are pretty baffling to me, since these are all composite and so iterating Miller-Rabin "should get out" pretty fast:
[Steven]
That's because you haven't seen my code :-)
I'd be baffled anyway ;-) about 4.7 seconds to test 2**800 + 1;
I got 1.71 msec, over a thousand times faster.
less than a tenth of a millisecond to test 2**900 + 1;
I got 2.32 msec, over 20 times _slower_(!). and about 8.6 seconds to test 2**1000 + 1. And I got 3.14 msec, again over a thousand times faster. It's over-engineered, class-based, and written as a learning exercise.
While my code had no fluff at all. It's hard to believe you could add enough cruft to slow it down by a factor of 1000, but pretty much impossible to believe adding cruft could make it way faster in the middle case above. Or do you first try division by small primes? 2**900+1 could get out cheap in that case, because it happens to be divisible by 17.
No argument from me. There's a reason I wrote my `pp()` to begin with too -) State-of-the-art code for these things requires planet-sized brains and commitment, but even a little expert knowledge can yield simple code that's an enormous improvement over any "dead obvious" approach.
On Sat, Jul 14, 2018 at 03:39:25PM -0500, Tim Peters wrote:
While my code had no fluff at all. It's hard to believe you could add enough cruft to slow it down by a factor of 1000,
There's a factor of 20 because my computer is older and slower than yours. (I ran your version, and it was ~20 times slower on my computer than the numbers you give.) So the cruft factor is only 50 times :)
Yes indeed. My approach is: 1. trial division by small primes; 2. Miller-Rabin by a minimal set of provably deterministic witnesses (if such a set exists, which it does for n < 2**64); 3. or by the first twelve primes if no such known set exists; 4. followed by 30 random witnesses. 12+30 = 42 Miller-Rabin tests may be overkill :-) -- Steve
On Fri, Jul 13, 2018 at 10:10:00AM +0200, Jeroen Demeyer wrote:
Go on, rub it in why don't you... I'm pretty sure that Sage is using numpy, which would almost certainly have a C-based implementation, not pure Python like mine. And it's probably making far less conservative choices about the rate of false positives it is willing to accept. And your computer is probably newer and faster than mine :-( My point was that for casual use, even the bad cases are still fast enough for interactive use. If you think 5 seconds is ridiculously slow, you've never tried a naive trial division algorithm that took literally days to crawl through the millions of divisions needed. *wink* -- Steve
[Steven D'Aprano] about 4.7 seconds to test 2**800 + 1;
[Jeroen Demeyer]
More like 6, yes? My Python version is more like 3, but is way fast enough for most things I do: $ py -m timeit -s "from temp3 import pp; n=2**800+1" "pp(n)" 200 loops, best of 5: 1.71 msec per loop $ py -m timeit -s "from temp3 import pp; n=2**900+1" "pp(n)" 100 loops, best of 5: 2.32 msec per loop $ py -m timeit -s "from temp3 import pp; n=2**1000+1" "pp(n)" 100 loops, best of 5: 3.15 msec per loop Steven's numbers are pretty baffling to me, since these are all composite and so iterating Miller-Rabin "should get out" pretty fast: about 4.7 seconds to test 2**800 + 1;
less than a tenth of a millisecond to test 2**900 + 1; and about 8.6 seconds to test 2**1000 + 1.
Here's the Python code I used: def pp(n, k=10): """ Return True iff n is a probable prime, using k Miller-Rabin tests. When False, n is definitely composite. """ from random import randrange if n < 4: return 2 <= n <= 3 d = nm1 = n-1 s = 0 while d & 1 == 0: s += 1 d >>= 1 if s == 0: return False # n >= 4 and even ss = range(1, s) for _ in range(k): a = randrange(2, nm1) x = pow(a, d, n) if x == 1 or x == nm1: continue for _ in ss: x = x*x % n if x == nm1: break if x == 1: return False else: return False return True
If we have a fast and exact primality test up to 2**64, a signature like this would help: def is_prime(n, probabilistic=False): When probabilistic is False and the number is beyond 2**64, it'll return OverflowError. When probabilistic is True, it's unbounded, but probabilistic. The docs should be explicit about that, though. Instead of an OverflowError, a number beyond 2**64 might simply use a slower strategy. But to make it an explicit option, the signature could be: def is_prime(n, exact=True, unbounded=False): By tellling it should be both exact and unbounded, the user is aware that it might be slow. The alternatives are: - exact, unbounded: Slower test - not exact, unbounded: Probabilistic test - exact, bounded: Default, raises OverflowError beyond 2**64 - not exact, bounded: Invalid, but it can just ignore exact input -- Danilo J. S. Bellini --------------- "*It is not our business to set up prohibitions, but to arrive at conventions.*" (R. Carnap)
[Tim]
More like 6, yes?
Heh - sorry about that. A speck of dirt on my monitor made me read '303" as "3.03".
My Python version is more like 3, but is way fast enough for most things I do:
So my Python Miller-Rabin code (already shared) was about one order of magnitude slower. ...
On Fri, Jul 13, 2018 at 10:52:38AM -0500, Tim Peters wrote:
Steven's numbers are pretty baffling to me, since these are all composite and so iterating Miller-Rabin "should get out" pretty fast:
That's because you haven't seen my code :-) It's over-engineered, class-based, and written as a learning exercise. There's a compatibility layer so it will work back to Python 2.4 and an instrumentation layer which I inserted in a fit of enthusiasm to gather staticistics, which I have never once looked at since. And *even so* it is still fast enough for casual use at the interactive interpreter, compared to more naive algorithms with worse Big O performance characteristics. You might think 5 seconds is slow, but I'm serious when I say some of the other algorithms I played with would take *days* to generate, or check, largish primes. -- Steve
Looking for some information on how long it has taken to generate large primes I stumbled across https://arxiv.org/ftp/arxiv/papers/1709/1709.09963.pdf which makes interesting reading. It concentrates on giving no false negatives (never saying n is not a prime when it is) but giving an increasing probability that n is a prime as testing goes on. I did find an article from 2016 that mentioned M77232917 (a 23,249,425 digit Mersenne prime) and M74207281 (22,338,618 digit Mersenne prime) with the latter taking a month to check with the Lucas-Lehmer test. On 14/07/2018 05:54, Steven D'Aprano wrote:
-- Steve (Gadget) Barnes Any opinions in this message are my personal opinions and do not reflect those of my employer. --- This email has been checked for viruses by AVG. https://www.avg.com
[Steve Barnes]
That's a very nice, gentle introduction! Thanks for sharing it. The `pp()` function I posted before is a Python implementation of the Miller-Rabin test they ended with. The latter is very widely used, because the code is relatively simple. short, and fast, and has no "bad cases". For _generating_ primes, the article uses other tricks to weed out sure-to-be-composite candidates before bothering with probable-prime testing. But the tricks they use in the paper are specific to a base 10 representation. In binary, it's common to skip candidates such that gcd(candidate, 2 * 3 * 5 * ...) > 1 where the product-of-small-primes second argument is picked to be the largest such that bigint division by it is exceptionally fast (usually one "digit" in whatever internal power-of-2 base is used by the bigint implementation). I did find an article from 2016 that mentioned M77232917 (a 23,249,425
digit Mersenne prime) and M74207281 (22,338,618 digit Mersenne prime) with the latter taking a month to check with the Lucas-Lehmer test.
More, it's a special version of Lucas-Lehmer (LL) that only works for testing candidates of the form 2**p-1 where p is itself an odd prime. In that context, it's not probabilistic - it definitively answers whether 2**p-1 is prime. Fast as it is, Miller-Rabin is impractical for guessing about numbers of this size. And as enormously faster as _it_ is, the special version of LL used for Mersenne testing is largely hand-coded in assembler to squeeze out every possible machine cycle. I don't want any of that in Python's standard library either ;-) Here! You should enjoy this account of the even larger Mersenne prime discovered this year: https://www.mersenne.org/primes/press/M77232917.html Note that the hyper-optimized testing code they use is so historically error-prone (& so good at provoking HW errors on overheated machines!) that they now verify each new Mersenne prime by 4 entirely different LL implementations, on 4 different machines.
[Tim]
Steven's numbers are pretty baffling to me, since these are all composite and so iterating Miller-Rabin "should get out" pretty fast:
[Steven]
That's because you haven't seen my code :-)
I'd be baffled anyway ;-) about 4.7 seconds to test 2**800 + 1;
I got 1.71 msec, over a thousand times faster.
less than a tenth of a millisecond to test 2**900 + 1;
I got 2.32 msec, over 20 times _slower_(!). and about 8.6 seconds to test 2**1000 + 1. And I got 3.14 msec, again over a thousand times faster. It's over-engineered, class-based, and written as a learning exercise.
While my code had no fluff at all. It's hard to believe you could add enough cruft to slow it down by a factor of 1000, but pretty much impossible to believe adding cruft could make it way faster in the middle case above. Or do you first try division by small primes? 2**900+1 could get out cheap in that case, because it happens to be divisible by 17.
No argument from me. There's a reason I wrote my `pp()` to begin with too -) State-of-the-art code for these things requires planet-sized brains and commitment, but even a little expert knowledge can yield simple code that's an enormous improvement over any "dead obvious" approach.
On Sat, Jul 14, 2018 at 03:39:25PM -0500, Tim Peters wrote:
While my code had no fluff at all. It's hard to believe you could add enough cruft to slow it down by a factor of 1000,
There's a factor of 20 because my computer is older and slower than yours. (I ran your version, and it was ~20 times slower on my computer than the numbers you give.) So the cruft factor is only 50 times :)
Yes indeed. My approach is: 1. trial division by small primes; 2. Miller-Rabin by a minimal set of provably deterministic witnesses (if such a set exists, which it does for n < 2**64); 3. or by the first twelve primes if no such known set exists; 4. followed by 30 random witnesses. 12+30 = 42 Miller-Rabin tests may be overkill :-) -- Steve
participants (7)
-
Danilo J. S. Bellini -
George Leslie-Waksman -
Greg Ewing -
Jeroen Demeyer -
Steve Barnes -
Steven D'Aprano -
Tim Peters