Pseudoprimes - Page 3

J.F. · 2009-06-02, 19:00

Time for an update.

I have teamed up with Will Galway, we are still working on completing the 64-bit search. In response to a request of Ernst: the 10^17 search is now complete!
I've placed a database + some stats here:
http://www.janfeitsma.nl/math/psp2/database

Readers are invited to double-check ranges - please inform me which range you have double-checked.

CRGreathouse · 2009-06-02, 19:27

Are you submitting

Code:

%H A055550 Jan Feitsma, <a href="http://www.janfeitsma.nl/math/psp2/database">Pseudoprime database</a>
%S A055550 0, 0, 3, 22, 78, 245, 750, 2057, 5597, 14884, 38975, 101629, 264239, 687007, 1801533, 4744920, 12603988

to Dr. Sloane?

J.F. · 2009-06-02, 21:47

That's a bit premature, I think. Better to do so while referring to a website with details about the actual computation, or a paper. So I'll keep this in mind.

Btw, that sequence currently does not list Galway's counts (terms 14 and 15), do you happen to know why?

CRGreathouse · 2009-06-03, 00:11

Quote:

Originally Posted by J.F.

Btw, that sequence currently does not list Galway's counts (terms 14 and 15), do you happen to know why?

I don't. I'm trying to remember if it ever did -- somehow I thought so, and the entry intimates as much. But I'm not sure.

flouran · 2009-06-04, 04:37

Quote:

Originally Posted by J.F.

Time for an update.

I have teamed up with Will Galway, we are still working on completing the 64-bit search. In response to a request of Ernst: the 10^17 search is now complete!
I've placed a database + some stats here:
http://www.janfeitsma.nl/math/psp2/database

Readers are invited to double-check ranges - please inform me which range you have double-checked.

I am a bit inexperienced in the topic, J.F., but do you have a copy of the source-code of the software program you used to count the pseudoprimes? You can PM me if you want.

J.F. · 2009-06-04, 14:23

Well, we are not ready to release the used source code in public. If we made it available now, it would be difficult to use as some sub-programs require quite large datafiles (8GB and some); plus there are compiling issues (optimized assembly code; MPI).

I'd also rather see doublechecking with independent programs. Can you program yourself? If so, and if you have considerable computing resources, then running the following would be useful:

Code:

for some r < 10^5, for which the Mersenne number 2^r-1 is not completely factored
   set plim to 2^64/(r+1) if r is even
   set plim to 2^64/(2r+1) if r is odd
   for all primes p with p%r == 1 and p <= plim
      if powermod(2, r-1, p) == 1
         if multiplicativeorder(2, p) == r
            report (r, p)
         end
      end
   end
end

WraithX · 2009-06-13, 13:58

Ok, so after programming the above pseudocode myself, and talking with J.F. about this, we've found there needs to be one small change made to the pseudocode to make it valid. The r-1 in the powermod needs to be r, like so:

Code:

for some r < 10^5, for which the Mersenne number 2^r-1 is not completely factored
   set plim to 2^64/(r+1) if r is even
   set plim to 2^64/(2r+1) if r is odd
   for all primes p with p%r == 1 and p <= plim
      if powermod(2, r, p) == 1
         if multiplicativeorder(2, p) == r
            report (r, p)
         end
      end
   end
end

Now, with this, here are some results found recently.

Those found by J.F. with 90100 <= r <= 90110 (with his own program)
90100 9020473494301
90100 1066333501
90100 2081940701
90100 177744504701
90102 270307
90102 1292783497
90105 8614038001
90105 903933361
90106 2252651
90106 90107
90109 429882465647

A few that I have found:
90111 none
90113 none
90115 none
90117 720936001
90119 180239
90121 none

The reason I am running on odd values of r is because it takes about 4 times longer to check even values of r with my program. Here are my timings so far:
r=90100 runtime= 0d 3h 37m 28.9470s
r=90111 runtime= 0d 1h 8m 33.5780s
r=90113 runtime= 0d 0h 40m 14.5620s
r=90115 runtime= 0d 0h 50m 7.1710s
r=90117 runtime= 0d 1h 7m 36.7120s
r=90119 runtime= 0d 0h 40m 5.0100s
r=90121 runtime= 0d 0h 39m 58.2910s

I will post more later.

Jeff Gilchrist · 2009-06-16, 17:14

Thanks to J.F.'s database of Fermat base 2 pseudoprimes up to 10^17 I have been able to confirm that does not exist any Baille-PSW pseudoprimes below 10^17: http://gilchrist.ca/jeff/factoring/pseudoprimes.html

Now to start working on helping him extend the database...

J.F. · 2013-03-01, 15:00

It has been far too long since Galway and I 'published' the 2^64 (S)PSP results [1,2] without proof.
I feel embarrassed when I see people using the list, while running the risk that it is not correct.

With this post I hope to resolve the 'preliminary' status and facilitate a proper acceptance of the result.
I put down my algorithm notes at [1] in the hope that someone finds the time to verify them to be correct, or spot a mistake. In case a mistake is found, I imagine this should only have minor impact on the published pseudoprime list.

History, background.

In 2008 I stumbled upon some work from Galway, which caught my interest hobby wise. He presented an extension of the base-2 PSP list up to 10^15. I modified his algorithms and tried to prove correctness of the result. In the mean time, I used the RUG computing facilities to create the 2^64 list.

Unfortunately, Galway and myself never found the time to finalize the project in the form of a paper.
I am largely to blame for this - I considered the result satisfying, even though it lacked closure in the form of a rigourous proof.
Other hobbies and my new career took up most of my attention, I just did not want to spend any significant effort on a paper.

Please accept my apologies and try to make the best of what I have put on my website.

Sincerely,
Jan Feitsma

Links (I was surprised by the number of references to my work, hence this letter).

[1] The pseudoprimes up to 2^64.
http://www.janfeitsma.nl/math/psp2/index
This is my playground website, which has finally been updated with the algorithms used to form the resulting list.

[2] Tables of pseudoprimes and related data
http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html
Galway graciously hosts the complete list since it wouldn't fit in my corner of the web.

[3] Computing the pseudoprimes up to 10^13
http://www.nt.ntnu.no/users/skoge/pr...f/31_Mikus.pdf
Counting base-2 PSP and base-3 PSP.

[4] The best known SPRP bases sets
http://miller-rabin.appspot.com/
A rather new website listing the current records, aiming to improve on them.

[5] Pseudoprime Enumeration with Probabilistic Primality Tests
http://gilchrist.ca/jeff/factoring/pseudoprimes.html
Note: Jeff Gilchrist and David Cleaver have spent quite some effort to reproduce my results, assuming the algorithms would be correct.
Many thanks again!

[6] A055550: Number of Poulet numbers (or pseudoprimes to base 2, A001567) less than 10^n.
http://oeis.org/A055550

[7] Sneaky Composites
http://wstein.org/edu/2010/414/projects/junge.pdf
Hey I am not a woman! Jan is a male name in The Netherlands. :)

[8] Pseudoprimes
http://www.mersenneforum.org/showthread.php?t=11443
I am "J.F." by the way.

danaj · 2013-06-23, 19:41

On the Pseudoprime Enumeration with Probabilistic Primality Tests page, the "Extra Strong Lucas" numbers are from T.R. Nicely's implementation. I don't have the Mo/Jones preprint, but I have looked at Grantham (1998 and 2000) and it isn't immediately obvious to me that Nicely's implementation follows. Nicely indicates that he chooses D when Jacobi(D,n) != 0, and then explains for 4 paragraphs of comments why this gives unsatisfactory results. However one can use Jacobi(D,n) = -1 instead: see OEIS A217719 -- start with P=3 Q=1 and increase P until (D|n) = -1, also checking the gcd. Wikipedia's current Lucas pseudoprime article uses the same process.

No PSP-2 (or SPSP-2) up to 2^64 is an extra strong Lucas pseudoprime with this method. As Nicely found with his code, the ES test runs about 20-25% faster because we can skip a mulmod per iteration.

Code:

         PSP-2  SPSP-2     LSP    SLSP    ESLP
1e 9      5597    1282    5485    1415     943
1e10     14884    3291   15352    3622    2346
1e11     38975    8607   42505    9714    6235
1e12    101629   22407  116928   25542   16231

I also found the numbers starting to differ for the "# PSP Miller-Rabin base 2" starting at 4e11. I've tried my C and C+GMP implementation as well as two other people's and all of them give me 15325 SPRP-2's below 4.0e11 rather than the 15758 in that page's table. After more experimentation all giving me the same results, I looked at the big downloadable annotated table and it agrees with my numbers. I think something is wrong with the second "Pseudoprime Table" column 2 (#PSP Miller-Rabin base 2) -- starting at 4.0e11 the numbers are all 433 too high.

joe53 · 2013-07-06, 12:07

My formula for the number of strong PSP-2 up to N:
z0(N) = c * sqrt(N)/y/y/v where c=56 approximately.
y means log N, v means log log N.
z means the actual value the number of strong PSP-2 up to N
with c=55.4 :
N z0(N) z(N) millions
---------------------------
10^12 0.022 0.022
10^13 0.058 0.059
10^14 0.153 0.156
10^15 0.415 0.419
10^16 1.132 1.136
10^17 3.118 3.115
10^18 8.659 8.647
10^19 24.225 24.220
on the other hand:
N z c
------------------------------
2^35 5511 55.799
2^40 23270 56.678
2^45 100342 56.619
2^50 441270 56.005
2^55 1988905 55.445
2^60 9212942 55.323
2^61 12552513 55.336
2^62 17114780 55.352
2^63 23355139 55.383
2^64 31894014 55.421
The number of primes up to N: Pi(N) = N/log N approx.
Does it means that if I take a m-digit (m>=8) random odd number, and passed the Sprp-2 test, then the probability(P) of composite less then
P < 1/10^(m/2) ? ( or P(N) < 1/sqrt(N) )
If passed then prime or spsp:
P= #spsp(N) / (#spsp(N)+Pi(N))
Strong pseudoprimes-2 counts
limit #spsp Pi(N) P
-------------------------------------------
1e8 488 5761455 8,469e-5 < 0,0001
1e10 3291 455052511 7,232e-6 < 0,00001
1e12 22407 37607912018 5,958e-7 < 0,000001
1e14 156251 3204941750802 4,875e-8 < 0,0000001
1e16 1135860 279238341033925 4,068e-9 < 0,00000001
1e18 8646507 24739954287740860 3,495e-10 < 0,000000001
If it's true(?) for big numbers then
if I take a 1000 digit random odd number(n), and passed the Sprp-2 test,
then the P(n composite) < 1/10^500 ?
If I do other 10 Sprp test, then P< 1/10^500 * 1/4^10 ~ 1/10^506 ?

2009-06-04, 14:23	#28
J.F. Jun 2008 2³×3² Posts	Well, we are not ready to release the used source code in public. If we made it available now, it would be difficult to use as some sub-programs require quite large datafiles (8GB and some); plus there are compiling issues (optimized assembly code; MPI). I'd also rather see doublechecking with independent programs. Can you program yourself? If so, and if you have considerable computing resources, then running the following would be useful: Code: for some r < 10^5, for which the Mersenne number 2^r-1 is not completely factored set plim to 2^64/(r+1) if r is even set plim to 2^64/(2r+1) if r is odd for all primes p with p%r == 1 and p <= plim if powermod(2, r-1, p) == 1 if multiplicativeorder(2, p) == r report (r, p) end end end end

2009-06-13, 13:58	#29
WraithX Mar 2006 111011111₂ Posts	Ok, so after programming the above pseudocode myself, and talking with J.F. about this, we've found there needs to be one small change made to the pseudocode to make it valid. The r-1 in the powermod needs to be r, like so: Code: for some r < 10^5, for which the Mersenne number 2^r-1 is not completely factored set plim to 2^64/(r+1) if r is even set plim to 2^64/(2r+1) if r is odd for all primes p with p%r == 1 and p <= plim if powermod(2, r, p) == 1 if multiplicativeorder(2, p) == r report (r, p) end end end end Now, with this, here are some results found recently. Those found by J.F. with 90100 <= r <= 90110 (with his own program) 90100 9020473494301 90100 1066333501 90100 2081940701 90100 177744504701 90102 270307 90102 1292783497 90105 8614038001 90105 903933361 90106 2252651 90106 90107 90109 429882465647 A few that I have found: 90111 none 90113 none 90115 none 90117 720936001 90119 180239 90121 none The reason I am running on odd values of r is because it takes about 4 times longer to check even values of r with my program. Here are my timings so far: r=90100 runtime= 0d 3h 37m 28.9470s r=90111 runtime= 0d 1h 8m 33.5780s r=90113 runtime= 0d 0h 40m 14.5620s r=90115 runtime= 0d 0h 50m 7.1710s r=90117 runtime= 0d 1h 7m 36.7120s r=90119 runtime= 0d 0h 40m 5.0100s r=90121 runtime= 0d 0h 39m 58.2910s I will post more later.

2013-06-23, 19:41	#32
danaj "Dana Jacobsen" Feb 2011 Bangkok, TH 1614₈ Posts	On the Pseudoprime Enumeration with Probabilistic Primality Tests page, the "Extra Strong Lucas" numbers are from T.R. Nicely's implementation. I don't have the Mo/Jones preprint, but I have looked at Grantham (1998 and 2000) and it isn't immediately obvious to me that Nicely's implementation follows. Nicely indicates that he chooses D when Jacobi(D,n) != 0, and then explains for 4 paragraphs of comments why this gives unsatisfactory results. However one can use Jacobi(D,n) = -1 instead: see OEIS A217719 -- start with P=3 Q=1 and increase P until (D\|n) = -1, also checking the gcd. Wikipedia's current Lucas pseudoprime article uses the same process. No PSP-2 (or SPSP-2) up to 2^64 is an extra strong Lucas pseudoprime with this method. As Nicely found with his code, the ES test runs about 20-25% faster because we can skip a mulmod per iteration. Code: PSP-2 SPSP-2 LSP SLSP ESLP 1e 9 5597 1282 5485 1415 943 1e10 14884 3291 15352 3622 2346 1e11 38975 8607 42505 9714 6235 1e12 101629 22407 116928 25542 16231 I also found the numbers starting to differ for the "# PSP Miller-Rabin base 2" starting at 4e11. I've tried my C and C+GMP implementation as well as two other people's and all of them give me 15325 SPRP-2's below 4.0e11 rather than the 15758 in that page's table. After more experimentation all giving me the same results, I looked at the big downloadable annotated table and it agrees with my numbers. I think something is wrong with the second "Pseudoprime Table" column 2 (#PSP Miller-Rabin base 2) -- starting at 4.0e11 the numbers are all 433 too high.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Pseudoprimes in special fields	devarajkandadai	Number Theory Discussion Group	7	2017-12-06 01:46
On generating Strong PseudoPrimes DataBase?	TheoryQuest1	Factoring	10	2016-09-19 16:08
Odd Fibonacci pseudoprimes	efeuvete	Math	7	2013-05-26 11:24

2009-06-02, 19:00	#23
J.F. Jun 2008 72₁₀ Posts	Time for an update. I have teamed up with Will Galway, we are still working on completing the 64-bit search. In response to a request of Ernst: the 10^17 search is now complete! I've placed a database + some stats here: http://www.janfeitsma.nl/math/psp2/database Readers are invited to double-check ranges - please inform me which range you have double-checked.

2009-06-02, 21:47	#25
J.F. Jun 2008 2³·3² Posts	That's a bit premature, I think. Better to do so while referring to a website with details about the actual computation, or a paper. So I'll keep this in mind. Btw, that sequence currently does not list Galway's counts (terms 14 and 15), do you happen to know why?

2009-06-16, 17:14	#30
Jeff Gilchrist Jun 2003 Ottawa, Canada 3·17·23 Posts	Thanks to J.F.'s database of Fermat base 2 pseudoprimes up to 10^17 I have been able to confirm that does not exist any Baille-PSW pseudoprimes below 10^17: http://gilchrist.ca/jeff/factoring/pseudoprimes.html Now to start working on helping him extend the database...

2013-03-01, 15:00	#31
J.F. Jun 2008 1001000₂ Posts	It has been far too long since Galway and I 'published' the 2^64 (S)PSP results [1,2] without proof. I feel embarrassed when I see people using the list, while running the risk that it is not correct. With this post I hope to resolve the 'preliminary' status and facilitate a proper acceptance of the result. I put down my algorithm notes at [1] in the hope that someone finds the time to verify them to be correct, or spot a mistake. In case a mistake is found, I imagine this should only have minor impact on the published pseudoprime list. History, background. In 2008 I stumbled upon some work from Galway, which caught my interest hobby wise. He presented an extension of the base-2 PSP list up to 10^15. I modified his algorithms and tried to prove correctness of the result. In the mean time, I used the RUG computing facilities to create the 2^64 list. Unfortunately, Galway and myself never found the time to finalize the project in the form of a paper. I am largely to blame for this - I considered the result satisfying, even though it lacked closure in the form of a rigourous proof. Other hobbies and my new career took up most of my attention, I just did not want to spend any significant effort on a paper. Please accept my apologies and try to make the best of what I have put on my website. Sincerely, Jan Feitsma Links (I was surprised by the number of references to my work, hence this letter). [1] The pseudoprimes up to 2^64. http://www.janfeitsma.nl/math/psp2/index This is my playground website, which has finally been updated with the algorithms used to form the resulting list. [2] Tables of pseudoprimes and related data http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html Galway graciously hosts the complete list since it wouldn't fit in my corner of the web. [3] Computing the pseudoprimes up to 10^13 http://www.nt.ntnu.no/users/skoge/pr...f/31_Mikus.pdf Counting base-2 PSP and base-3 PSP. [4] The best known SPRP bases sets http://miller-rabin.appspot.com/ A rather new website listing the current records, aiming to improve on them. [5] Pseudoprime Enumeration with Probabilistic Primality Tests http://gilchrist.ca/jeff/factoring/pseudoprimes.html Note: Jeff Gilchrist and David Cleaver have spent quite some effort to reproduce my results, assuming the algorithms would be correct. Many thanks again! [6] A055550: Number of Poulet numbers (or pseudoprimes to base 2, A001567) less than 10^n. http://oeis.org/A055550 [7] Sneaky Composites http://wstein.org/edu/2010/414/projects/junge.pdf Hey I am not a woman! Jan is a male name in The Netherlands. :) [8] Pseudoprimes http://www.mersenneforum.org/showthread.php?t=11443 I am "J.F." by the way.

2013-07-06, 12:07	#33
joe53 Jun 2013 2²·3 Posts	My formula for the number of strong PSP-2 up to N: z0(N) = c * sqrt(N)/y/y/v where c=56 approximately. y means log N, v means log log N. z means the actual value the number of strong PSP-2 up to N with c=55.4 : N z0(N) z(N) millions --------------------------- 10^12 0.022 0.022 10^13 0.058 0.059 10^14 0.153 0.156 10^15 0.415 0.419 10^16 1.132 1.136 10^17 3.118 3.115 10^18 8.659 8.647 10^19 24.225 24.220 on the other hand: N z c ------------------------------ 2^35 5511 55.799 2^40 23270 56.678 2^45 100342 56.619 2^50 441270 56.005 2^55 1988905 55.445 2^60 9212942 55.323 2^61 12552513 55.336 2^62 17114780 55.352 2^63 23355139 55.383 2^64 31894014 55.421 The number of primes up to N: Pi(N) = N/log N approx. Does it means that if I take a m-digit (m>=8) random odd number, and passed the Sprp-2 test, then the probability(P) of composite less then P < 1/10^(m/2) ? ( or P(N) < 1/sqrt(N) ) If passed then prime or spsp: P= #spsp(N) / (#spsp(N)+Pi(N)) Strong pseudoprimes-2 counts limit #spsp Pi(N) P ------------------------------------------- 1e8 488 5761455 8,469e-5 < 0,0001 1e10 3291 455052511 7,232e-6 < 0,00001 1e12 22407 37607912018 5,958e-7 < 0,000001 1e14 156251 3204941750802 4,875e-8 < 0,0000001 1e16 1135860 279238341033925 4,068e-9 < 0,00000001 1e18 8646507 24739954287740860 3,495e-10 < 0,000000001 If it's true(?) for big numbers then if I take a 1000 digit random odd number(n), and passed the Sprp-2 test, then the P(n composite) < 1/10^500 ? If I do other 10 Sprp test, then P< 1/10^500 * 1/4^10 ~ 1/10^506 ?