20160916, 07:11  #1 
Sep 2016
1_{8} Posts 
On generating Strong PseudoPrimes DataBase?
Is there a program/library avalible to generate a DataBase of Random Strong PseudoPrimes with respect to a Input Base (the BitSize of Numbers is an Input).
Of course one can code it, but I was wondering if a duplication of effort can be avoided ? 
20160916, 10:04  #2 
Romulan Interpreter
Jun 2011
Thailand
5^{2}·7·53 Posts 
pari/gp.
You can use the ispseudoprime(x,k) function to generate strong probable primes to random k bases and print those which are not prime (the name of the function is misleading, it will pass primes too). Something like: Code:
n=3; cnt=0; while(n<100000, if(ispseudoprime(n,1) && !isprime(n), print(n)); n+=2) Code:
n=3; cnt=0; while(n<1000000, if(Mod(2,n)^(n1)==1 && ispseudoprime(n,1) && !isprime(n), print(n)); n+=2) Code:
/*check if a number is bPRP, if b missing then do 2PRP test*/ isPRP(n=2,b=2)= { return(Mod(b,n)^(n1)==1); } vecSPRP(n=2,b=2)= { my(v,s,d); d=n1; s=0; while(!(d%2),d>>=1;s++); v=vector(s+5); v[1]=n; v[2]=b; v[3]=s; v[4]=d; v[5]=lift(d=Mod(b,n)^d); if(v[5]==n1,v[5]=1); for(r=1,s, v[r+5]=lift(d=d^2); if(v[r+5]==n1,v[r+5]=1); ); return(v); } isSPRP(n=2,b=2)= { my(s,d); d=n1; s=0; while(!(d%2),d>>=1;s++); d=Mod(b,n)^d; if(d==1  d==1, return(1) ); for(r=1,s1, if((d=d^2)==1, return(1), if(d==1, return(0) ) ) ); return(0); } isPSP(n=2,b=2)= { return(isPRP(n,b) & !isprime(n)); } isSPSP(n=2,b=2)= { return(isSPRP(n,b) & !isprime(n)); } /*very slow, factors the number and uses Korselt's criterion to check if Carmichael*/ isCarm(n=2)= { if(!issquarefree(n), return(0)); fct=factorint(n); for(i=1,#fct~,if((n1)%(fct[i,1]1), return(0))); return(1); } /*return next odd prp (slow!)*/ nextPRP(n=2,b=2)= { n+=1n%2; while(!isPRP(n,b),n+=2); return(n); } /*return next odd sprp (slow!)*/ nextSPRP(n=2,b=2)= { n+=1n%2; while(!isSPRP(n,b),n+=2); return(n); } /*return next "pure" pseudoprime (slow!)*/ nextPSP(n=2,b=2)= { n+=1n%2; while(!isPRP(n,b)isprime(n),n+=2); return(n); } /*return next "pure" strong pseudoprime (slow!)*/ nextSPSP(n=2,b=2)= { n+=1n%2; while(!isSPRP(n,b)isprime(n),n+=2); return(n); } /*complex mod function, must have second parameter as a positive integer*/ cmod(a=3+2*I,m)= { return(real(a)%m+imag(a)%m*I); } /*raise a complex number at an integer power modulo that integer, avoiding large numbers*/ /* n must be positive! and aa must be complex! it does not check its parameters, for speed reasons caution: n is the power, m is the mod! */ cpowmod(aa=3+2*I,n=3,m)= { if(m==0, m=n); if(n<10, return(cmod(aa^n,m)), /*else*/ return(cmod(cpowmod(aa^10,n\10,m)*cmod(aa^(n%10),m),m)) ); } /*complex base PRP test, b=base, n=number to test, positive integer*/ /*this would return true for all primes and few composite pseudoprimes*/ cbIsPrime(b=3+2*I,n)= { if(n%2==0, return(0)); if(n%4==3, return(cpowmod(b,n,n)==cmod(real(b)imag(b)*I,n)), /*else*/ return(cpowmod(b,n,n)==cmod(real(b)+imag(b)*I,n)) ); } /* nr=0; cnt=0; forstep(n=11,5*10^8,2,if(((z=(Mod(3,n)+Mod(2,n)*I)^(n+1))==13z==5+12*I)&&!isprime(n),nr++; print(n" : "n%4" ")); if(n>cnt,cnt+=10^6; printf("...%d...%c",n,13))); nr */ go_complex(n=1*10^10,limt=2*10^10)= { until(q==13n>limt, until((n=nextPSP(n+1,13))%4==3n>limt,); print(n", "q=lift((Mod(3,n)+Mod(2,n)*I)^(n+1))) ); } /*generates the pseudoprime list starting from n, write it in the file if asked */ genPSPlist(n=341, b=2, file="psp_2_list.txt")= { if(n%2==0,n++); while(1, while(!isPRP(n,b),n+=2); if(!isprime(n), /*a pseudeprime base b is found*/ if(isCarm(n),txtline="C",txtline=" "); if(isSPSP(n,b),txtline=Str(txtline,"S"),txtline=Str(txtline," ")); txtline=Str(txtline," "n" "); fct=factorint(n); for(i=1,#fct~,for(j=1,fct[i,2],txtline=Str(txtline,fct[i,1]" "))); print(txtline); if(file, write(file,txtline)) ); n+=2 ); } /*multibase PRP test, bv is a vector of bases*/ genMultiPSPlist(n=341, bvec=[2,3], file="psp_2_3_list.txt")= { cnt=0; if(n%2==0,n++); while(1, b=1; while(b<=#bvec,if(!isPRP(n,bvec[b]),break,b++)); if(b>#bvec && !isprime(n), /*a pseudeprime base all bvec[b] is found*/ fct=factorint(n); txtline=Str(n" "); for(i=1,#fct~,for(j=1,fct[i,2],txtline=Str(txtline,fct[i,1]" "))); print(txtline" "); if(file, write(file,txtline)) ); n+=2; if(cnt++>1000000,printf("... %d ...%c",n,13); cnt=0) ); } /*multibase SPRP test, bv is a vector of bases*/ /*ex: genMultiSPSPlist(,[31,73],0) would print 9080191, 15560651, etc */ genMultiSPSPlist(n=561, bvec=[2,3], file="spsp_2_3_list.txt")= { cnt=0; if(n%2==0,n++); while(1, b=1; while(b<=#bvec,if(!isSPRP(n,bvec[b]),break,b++)); if(b>#bvec && !isprime(n), /*a pseudeprime base all bvec[b] is found*/ fct=factorint(n); txtline=Str(n" "); for(i=1,#fct~,for(j=1,fct[i,2],txtline=Str(txtline,fct[i,1]" "))); print(txtline" "); if(file, write(file,txtline)) ); n+=2; if(cnt++>1000000,printf("... %d ...%c",n,13); cnt=0) ); } /* playing around gp > k=0; kk=0; n=1 %1 = 1 gp > until(k>kk, until(cnt>0, until(!isprime(n),n+=2); print("N="n"\t"factorint(n)); cnt=0; for(i=2,n/2,if(isSPRP(n,i),print(" : "i);cnt++ ))); k=[cnt,2*cnt/n,2.0*cnt/n][3]); kk=k N=9 Mat([3, 2]) N=15 [3, 1; 5, 1] N=21 [3, 1; 7, 1] N=25 Mat([5, 2]) : 7 %2 = 0.08000000000000000000000000000 /*some work example*/ /* gp > { b=3+2*I; cnt=0; nr=0; n=2; while((n=nextPSP(n++,13))<=10^7, if(((z=cpowmod(b,n+1,n))==13z==5+12*I)&&!isprime(n), nr++; print(n" : "n%4" ") ); if(n>cnt, cnt+=1000; printf("...%d...%c",n,13) ) ); print(" "); nr } */ /* start must be 3 (mod 4), step is 2 or 4, ex: go(11,10^6,2) */ go(start=11,stop=10^6,step=2)= { nr=0; cnt=0; forstep(n=start,stop,step, if(((z=(Mod(3,n)+Mod(2,n)*I)^(n+1))==13  z==5+12*I) && !isprime(n), nr++; print(n" : "n%4" ") ); if(n>cnt, cnt+=10^6; printf("...%d...%c",n,13) ) ); nr } Code:
gp > n=2; while(n<10^4, print(n=nextPSP(n,2)); n++) 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321 8481 8911 10261 time = 13 ms. gp > n=2; while(n<10^4, print(n=nextSPSP(n,2)); n++) 2047 3277 4033 4681 8321 15841 time = 27 ms. Code:
gp > genMultiSPSPlist(,[2,3],0) 1373653 829 1657 1530787 619 2473 1987021 997 1993 2284453 1069 2137 3116107 883 3529 5173601 929 5569 6787327 1303 5209 11541307 1699 6793 13694761 2617 5233 15978007 1999 7993 16070429 1637 9817 16879501 1453 11617 ... 18000579 ... Last fiddled with by LaurV on 20160916 at 10:20 
20160916, 13:21  #3 
"Dana Jacobsen"
Feb 2011
Bangkok, TH
2^{2}·227 Posts 
Code:
perl Mntheory=:all Mbigint E 'foroddcomposites { say if is_strong_pseudoprime($_, 17) } 2**35, 2**35 + 1e7' Also see: Pseudoprime Statistics, Tables, and Data For some source code:  https://github.com/danaj/MathPrime...er/primality.c  https://sourceforge.net/projects/mpzaprcl/files/  http://rosettacode.org/mw/index.php?...primality_test The speed of the strong pseudoprime test in the GMP internals, my code, and Cleaver's is essentially the same, as there's little room for changing anything if using GMP and its public interface. This is not the case with the Lucas test code however. In other library choices, gwnum (used in PFGW) is slower than GMP until ~2000 digits, but is faster above. Note that below 2^64, we already have the wonderful base2 pseudoprime database which has been used to create deterministic MR sets, hashed sets (e.g. one test for under 2^32, two tests below 2^49, three below 2^64), and of course verify BPSW as correct for all 64bit inputs. Also see the Sorenson/Webster 2015 arXiv paper on deterministic bases. For larger inputs, there are so many numbers that I'm not sure what the random SPRP data will really tell you. There are some papers that discuss probabilities at different bit sizes. BPSW is simple, in regular use, and still has no counterexamples known. Combine that with the probabilities for randombase MR tests at large sizes, and one can construct a very efficient and convincing compositeness test. 
20160917, 18:05  #4  
Bamboozled!
"πΊππ·π·π"
May 2003
Down not across
2^{5}×331 Posts 
Quote:
Coincidentally I've been using GMP through its Perl interface for testing composites. AFAIK, GMP does not use Montgomery arithmetic for its pseudoprimality test. Either I'm wrong or its speed could be (approximately) doubled. I should find out either way. 

20160917, 19:04  #5  
"Dana Jacobsen"
Feb 2011
Bangkok, TH
38C_{16} Posts 
Quote:
It could be sped up in a few ways, but I don't think there there are faster ways using the standard API. Please prove me wrong with code examples! 

20160917, 19:13  #6  
Bamboozled!
"πΊππ·π·π"
May 2003
Down not across
2^{5}·331 Posts 
Quote:
GMPECM has a Montgomery arithmetic implementation but it doesn't seem very easy to extract it for general use. 

20160917, 23:48  #7 
"Dana Jacobsen"
Feb 2011
Bangkok, TH
908_{10} Posts 
I see. GMPECM uses a lot of __gmpn_* internal functions inside its implementation, making it a maintenance issue for portability (which may or may not be important for ones application).
My 64bit code uses Montgomery math with a tiny bit of assembly for most of its primality tests, and it makes a big difference. But that's only good for 64bit. 
20160918, 17:27  #8 
Bamboozled!
"πΊππ·π·π"
May 2003
Down not across
10100101100000_{2} Posts 

20160918, 17:50  #9 
"Dana Jacobsen"
Feb 2011
Bangkok, TH
908_{10} Posts 
https://github.com/danaj/MathPrimeUtil
montmath.h has the core, written by Wojciech Izykowski. We worked together on the the inverse code. mulmod.h is the standard basic modular math. primality.[hc] has the various tests. I recently updated the Pollard/Brent routine in factor.c to use the mont code, which was a nice speedup (e.g. my code for OEIS A275598 ran about 2x faster in the 10^11+ range). 
20160918, 19:35  #10 
Just call me Henry
"David"
Sep 2007
Cambridge (GMT/BST)
11×23^{2} Posts 
@Danaj, would you be willing to consider adding GWNUM math to your library for large enough numbers. This should be possible using it as a library. This could help the larger searches for prime gaps.

20160919, 16:08  #11 
"Dana Jacobsen"
Feb 2011
Bangkok, TH
1110001100_{2} Posts 
To keep portability, I could do it via an added define, so one would edit the Makefile pointing to the library and setting the define. Probably just for PSP, SPSP, and Lucas (so pseudoprime, strong pseudoprime, and BPSW). It probably would still be a little slower than PFGW  how much is an interesting question.

Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Pseudoprimes in special fields  devarajkandadai  Number Theory Discussion Group  7  20171206 01:46 
Pseudoprimes  CRGreathouse  Computer Science & Computational Number Theory  36  20131121 07:47 
Odd Fibonacci pseudoprimes  efeuvete  Math  7  20130526 11:24 
Strong Pseudoprime Distribution  flouran  Math  3  20090601 05:04 
Strong Pseudoprime Question  amcfarlane  Math  35  20060902 23:02 