mersenneforum.org Fast isPrime() for n < 2^32
 Register FAQ Search Today's Posts Mark Forums Read

2011-04-21, 11:36   #23
R.D. Silverman

Nov 2003

22×5×373 Posts

Quote:
 Originally Posted by CRGreathouse Of course it matters for this thread, since that's the whole topic. But I would say that it matters to me personally: I frequently find myself in the position of needing to test the primality of many (trillions) of small numbers, usually below 64 bits. At the moment I use PARI which takes about 3 microseconds. (Well, if you abuse it enough you can get to that speed anyway; the current version doesn't yet know about Jan's result and so if you ask it to prove primality it takes more like a millisecond.) At those speeds even just a billion numbers take an hour. Reducing that to 15 or 20 minutes would be great.
Agreed.

 2011-04-22, 11:01 #24 Xitami   Apr 2010 2×7 Posts :) PARI/GP can write programs in C Code: ifelse(n,p,k)={ my(c); if(p0,print); for(i=0,n,print1(" ")); c=(p+k)>>1; print1("if(a1262["c"]
2011-04-22, 14:22   #25
Xitami

Apr 2010

2×7 Posts

Quote:
 Originally Posted by ATH [...]But I just took the list of the 31,894,014 SPRP base2 < 2^64 [...]
1. where i can find this list?
2. how many pseudoprimes survive test to two bases ?

 2011-04-22, 18:27 #26 ATH Einyen     Dec 2003 Denmark 1011010001002 Posts The full list is here: http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html but this include all 118,968,378 Fermat base2 PRP + Miller Rabin SPRP base2. The best constant second base I found so far is 9375 where 1,147,253 (0.9643%) survives. More exotic choices does a bit better. For number n if you use base: floor(n/5)+49 then 1,092,753 survives. Base: floor(Sqrt(n))+1: 1,023,767 survives.
2011-04-26, 13:14   #27
Xitami

Apr 2010

2·7 Posts

Quote:
 Originally Posted by zerothbase [...]This is definitely the weakest part of the entire code.[...]
Code:

Code:
typedef unsigned long long int u64;

inline u64 mulmod(u64 i, u64 j, u64 k){
u64 p, r = 0;
while(j > 0){
if(j&1) //  r = (i + r) % k
if(0==r) r=i;
else {  r=k-r;
if(i>=r) r=i-r;
else r=k-r+i;}
// i = (i + i) % k
if(i>k-i) i=(i-k)+i;
else i+=i;
j>>=1;}
return r;}
no MOD, no overflow

2011-04-27, 18:33   #28
Jeff Gilchrist

Jun 2003

100100100012 Posts

Quote:
 Originally Posted by ATH This has already been done: http://www.trnicely.net/misc/bpsw.html "Furthermore, preliminary analysis by Gilchrist of Feitsma's database, further extended, indicates that no Baillie-PSW pseudoprime (standard or strong) exists below 2^64 (24 October 2009; double checking of this result is in progress), approximately 1.8446744e19." But I just took the list of the 31,894,014 SPRP base2 < 2^64 and did a Lucas-Selfridge test on them, it took around 15min, and none of them are lucas PRP. So assuming doublechecking proves that the SPRP list is complete there is BPSW pseudoprime below 2^64.
Jan Feitsma has posted the results as "preliminary" because nobody has verified the math of his algorithm to ensure that all the base 2 psp up to 2^64 have been enumerated. David Cleaver (WraithX on these boards) and I verified the pseudoprime database with Jan's algorithm using independently written code so we have double checked that his algorithm generated the same set of psp using different code. But Jan has not proven that his algorithm in fact enumerates 100% of the psp up to 2^64. Jan was originally in contact with Wil Galway (the one hosting the 2^64 database now) to review his math but for whatever reason his algorithm was never evaluated.

Back in 2009 when David and I finished double checking the psp database, I was able to confirm there were no BPSW psp up to 2^64 with his database (assuming we did in fact enumerate all of them) and posted that here: http://gilchrist.ca/jeff/factoring/pseudoprimes.html

Charles Greathouse has also confirmed no BPSW psp exist using the BPSW implementation in PARI with Jan Feitsma's database so two different implementations of BPSW have also been tried.

Last fiddled with by Jeff Gilchrist on 2011-04-27 at 18:39

2011-04-28, 04:03   #29
CRGreathouse

Aug 2006

11·13·41 Posts

Quote:
 Originally Posted by Jeff Gilchrist Charles Greathouse has also confirmed no BPSW psp exist using the BPSW implementation in PARI with Jan Feitsma's database so two different implementations of BPSW have also been tried.
I know that many people have checked this, but I don't know which implementations of BPSW were used (though as you point out, it's at least two).

 2011-05-24, 16:57 #30 JohnFullspeed   May 2011 France 7·23 Posts 2^32 primes Why are you not using an array with all the primes smaller tahnn2^32? You have only 200 000 000 of primes so using a simple RLE the file is 200 000 000 of byte Two hours to compute the llst (i can give the file) and you can factorized 20 digits number in less than one minute. For the classic m^k mod you can use a fast algo Code: function E_M(m, K, N: int64): int64; var b, a: int64; begin b := 1; a := m; if odd(k) then b := A; k := k shr 1; repeat A := (A * A) mod N; if odd(k) then b := (A * b) mod N; k := k shr 1; until k = 0; result := b; end; It's also good for RSA encryption John Last fiddled with by JohnFullspeed on 2011-05-24 at 17:00
2011-05-24, 18:38   #31
CRGreathouse

Aug 2006

11×13×41 Posts

Quote:
 Originally Posted by JohnFullspeed Two hours to compute the llst (i can give the file)
Two hours?!? It takes PARI 11 seconds on this machine. (Technically that's "just" to 4.26e9 since this computer is running 32-bit gp.)

Quote:
 Originally Posted by JohnFullspeed and you can factorized 20 digits number in less than one minute.
Sure, as long as it's not a semiprime where both factors are between 2^32 and 10^10. But you can do much better -- using rho and SQUFOF, PARI can factor worst-case numbers of that size in about 200 milliseconds, and ten times faster for typical numbers.

But primality testing is much faster than factoring! PARI only takes about 2 milliseconds to prove that a number is prime; typical composites are extremely easy to detect, taking only 20 *microseconds* by my calculations (I had to test 100,000 to get good timing!).

So trial division works at that size, but we're looking for something really fast here.

 2011-05-24, 20:51 #32 Xitami   Apr 2010 E16 Posts A := (A * A) mod N; John, how long is A*A?
 2011-05-25, 07:32 #33 JohnFullspeed   May 2011 France 7·23 Posts a:=a*a mod N 10000 bits if you want.... In fact you can use this mrethod for othezr problems as divisioins.. The method is: Let A and B be big positive integerss Let ? an sllower opoeraareationn A ? B = C you work on the binary representation each bit of A If te bit is 0 you make a specific code if the bit is 1 you make and other code you pass to, the Bit after unil A have bit Code: function E_M(m, K, N: int64): int64; var b, a: int64; begin b := 1; a := m; if odd(k) then b := A; k := k shr 1; <= first bit of tha operande repeat A := (A * A) mod N; <= in all casae Bit=0 or Bit=1 if odd(k) then <= the bit is 1 b := (A * b) mod N; k := k shr 1; <+ nexT bit until k = 0; <= no new bit (end) result := b; end;` if you remove the mod you have a speed exopanciation Is't easy to do a faster div.... A := (A * A); <= in all casae Bit=0 or Bit=1 if odd(k) then <= the bit is 1 b := (A * b); k := k shr 1; <+ nexT bit until k = 0; <= no new bit (end) You understood that A can have 1000 bios you make 'only' 1000+5000 divmod John

 Similar Threads Thread Thread Starter Forum Replies Last Post jasong jasong 35 2016-12-11 00:57 FloatingPoint Operation Billion Digits 39 2015-10-21 02:15 Andi47 Puzzles 20 2009-04-01 02:35 ixfd64 Hardware 1 2005-11-21 21:28 maheshexp Math 2 2004-05-29 01:54

All times are UTC. The time now is 15:52.

Sat Aug 8 15:52:07 UTC 2020 up 22 days, 11:38, 2 users, load averages: 2.79, 2.62, 2.46