mersenneforum.org On generating Strong PseudoPrimes DataBase?
 Register FAQ Search Today's Posts Mark Forums Read

 2016-09-16, 07:11 #1 TheoryQuest1   Sep 2016 1 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 ?
 2016-09-16, 10:04 #2 LaurV Romulan Interpreter     "name field" Jun 2011 Thailand 2×5×1,013 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) If you run this 10 times will give you different numbers each time, because the bases are random. You can increase the number of bases, or use other limits. You also can "strength" it by adding a test for Sarrus numbers (weak 2-PRP), depend what do you want to use the numbers for. Code: n=3; cnt=0; while(n<1000000, if(Mod(2,n)^(n-1)==1 && ispseudoprime(n,1) && !isprime(n), print(n)); n+=2) If you need indeed "strong" PRPs, to a defined set of bases, you can make your own function very easy. Here are some toys I put together long time ago (without much optimizations). Code: /*check if a number is b-PRP, if b missing then do 2-PRP test*/ isPRP(n=2,b=2)= { return(Mod(b,n)^(n-1)==1); } vecSPRP(n=2,b=2)= { my(v,s,d); d=n-1; 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]==n-1,v[5]=-1); for(r=1,s, v[r+5]=lift(d=d^2); if(v[r+5]==n-1,v[r+5]=-1); ); return(v); } isSPRP(n=2,b=2)= { my(s,d); d=n-1; s=0; while(!(d%2),d>>=1;s++); d=Mod(b,n)^d; if(d==1 || d==-1, return(1) ); for(r=1,s-1, 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((n-1)%(fct[i,1]-1), return(0))); return(1); } /*return next odd prp (slow!)*/ nextPRP(n=2,b=2)= { n+=1-n%2; while(!isPRP(n,b),n+=2); return(n); } /*return next odd sprp (slow!)*/ nextSPRP(n=2,b=2)= { n+=1-n%2; while(!isSPRP(n,b),n+=2); return(n); } /*return next "pure" pseudoprime (slow!)*/ nextPSP(n=2,b=2)= { n+=1-n%2; while(!isPRP(n,b)||isprime(n),n+=2); return(n); } /*return next "pure" strong pseudoprime (slow!)*/ nextSPSP(n=2,b=2)= { n+=1-n%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))==13||z==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==13||n>limt, until((n=nextPSP(n+1,13))%4==3||n>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 ); } /*multi-base 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) ); } /*multi-base 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))==13||z==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 } You can use it like: 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. or: 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 2016-09-16 at 10:20
 2016-09-16, 13:21 #3 danaj   "Dana Jacobsen" Feb 2011 Bangkok, TH 32·101 Posts Code: perl -Mntheory=:all -Mbigint -E 'foroddcomposites { say if is_strong_pseudoprime(\$_, 17) } 2**35, 2**35 + 1e7' Change the base and endpoints as desired. You'll find 14 different pseudoprimality tests in that module, along with 3 general-form proof methods and a few special-form ones. Also see: Pseudoprime Statistics, Tables, and Data For some source code: - https://github.com/danaj/Math-Prime-...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 base-2 pseudoprime database which has been used to create deterministic M-R 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 64-bit 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 random-base M-R tests at large sizes, and one can construct a very efficient and convincing compositeness test.
2016-09-17, 18:05   #4
xilman
Bamboozled!

"𒉺𒌌𒇷𒆷𒀭"
May 2003
Down not across

2CF016 Posts

Quote:
 Originally Posted by danaj 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.

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.

2016-09-17, 19:04   #5
danaj

"Dana Jacobsen"
Feb 2011
Bangkok, TH

32×101 Posts

Quote:
 Originally Posted by xilman Are you sure about that? 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.
I have looked at the internals for GMP's mpz_probab_prime, and it just uses the public interfaces and very similar code (it tosses in a Fermat test at the start, which didn't measure any faster than a M-R test on my machines). There is no public interface for Montgomery math that I'm aware of. It's come up on various forums for years. There are some internal functions that could be used, albeit they could change in any release.

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!

2016-09-17, 19:13   #6
xilman
Bamboozled!

"𒉺𒌌𒇷𒆷𒀭"
May 2003
Down not across

24×719 Posts

Quote:
 Originally Posted by danaj I have looked at the internals for GMP's mpz_probab_prime, and it just uses the public interfaces and very similar code (it tosses in a Fermat test at the start, which didn't measure any faster than a M-R test on my machines). There is no public interface for Montgomery math that I'm aware of. It's come up on various forums for years. There are some internal functions that could be used, albeit they could change in any release. 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!
I was thinking more along the lines of using the public API to implement Montgomery arithmetic and thence Miller-Rabin.

GMP-ECM has a Montgomery arithmetic implementation but it doesn't seem very easy to extract it for general use.

 2016-09-17, 23:48 #7 danaj   "Dana Jacobsen" Feb 2011 Bangkok, TH 32×101 Posts I see. GMP-ECM 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 64-bit 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 64-bit.
2016-09-18, 17:27   #8
xilman
Bamboozled!

"𒉺𒌌𒇷𒆷𒀭"
May 2003
Down not across

24·719 Posts

Quote:
 Originally Posted by danaj My 64-bit 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 64-bit.
Is your code publicly available? I'd be interested in taking a look.

Thanks,

Paul

 2016-09-18, 17:50 #9 danaj   "Dana Jacobsen" Feb 2011 Bangkok, TH 32×101 Posts https://github.com/danaj/Math-Prime-Util 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).
 2016-09-18, 19:35 #10 henryzz Just call me Henry     "David" Sep 2007 Liverpool (GMT/BST) 7×857 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.
2016-09-19, 16:08   #11
danaj

"Dana Jacobsen"
Feb 2011
Bangkok, TH

32·101 Posts

Quote:
 Originally Posted by henryzz @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.
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.

 Similar Threads Thread Thread Starter Forum Replies Last Post devarajkandadai Number Theory Discussion Group 7 2017-12-06 01:46 CRGreathouse Computer Science & Computational Number Theory 36 2013-11-21 07:47 efeuvete Math 7 2013-05-26 11:24 flouran Math 3 2009-06-01 05:04 amcfarlane Math 35 2006-09-02 23:02

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

Thu Oct 6 15:03:03 UTC 2022 up 49 days, 12:31, 1 user, load averages: 1.19, 1.45, 1.49