mersenneforum.org  

Go Back   mersenneforum.org > Factoring Projects > Factoring

Reply
 
Thread Tools
Old 2016-09-16, 07:11   #1
TheoryQuest1
 
Sep 2016

1 Posts
Default 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 ?
TheoryQuest1 is offline   Reply With Quote
Old 2016-09-16, 10:04   #2
LaurV
Romulan Interpreter
 
LaurV's Avatar
 
Jun 2011
Thailand

22·7·11·29 Posts
Default

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
LaurV is offline   Reply With Quote
Old 2016-09-16, 13:21   #3
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

2·3·151 Posts
Default

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.
danaj is offline   Reply With Quote
Old 2016-09-17, 18:05   #4
xilman
Bamboozled!
 
xilman's Avatar
 
"π’‰Ίπ’ŒŒπ’‡·π’†·π’€­"
May 2003
Down not across

5×2,039 Posts
Default

Quote:
Originally Posted by danaj View Post
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.
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.
xilman is offline   Reply With Quote
Old 2016-09-17, 19:04   #5
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

90610 Posts
Default

Quote:
Originally Posted by xilman View Post
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!
danaj is offline   Reply With Quote
Old 2016-09-17, 19:13   #6
xilman
Bamboozled!
 
xilman's Avatar
 
"π’‰Ίπ’ŒŒπ’‡·π’†·π’€­"
May 2003
Down not across

5·2,039 Posts
Default

Quote:
Originally Posted by danaj View Post
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.
xilman is offline   Reply With Quote
Old 2016-09-17, 23:48   #7
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

90610 Posts
Default

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.
danaj is offline   Reply With Quote
Old 2016-09-18, 17:27   #8
xilman
Bamboozled!
 
xilman's Avatar
 
"π’‰Ίπ’ŒŒπ’‡·π’†·π’€­"
May 2003
Down not across

27D316 Posts
Default

Quote:
Originally Posted by danaj View Post
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
xilman is offline   Reply With Quote
Old 2016-09-18, 17:50   #9
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

90610 Posts
Default

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).
danaj is offline   Reply With Quote
Old 2016-09-18, 19:35   #10
henryzz
Just call me Henry
 
henryzz's Avatar
 
"David"
Sep 2007
Cambridge (GMT/BST)

34·71 Posts
Default

@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.
henryzz is offline   Reply With Quote
Old 2016-09-19, 16:08   #11
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

2·3·151 Posts
Default

Quote:
Originally Posted by henryzz View Post
@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.
danaj is offline   Reply With Quote
Reply

Thread Tools


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

All times are UTC. The time now is 22:35.

Fri Nov 27 22:35:20 UTC 2020 up 78 days, 19:46, 3 users, load averages: 1.45, 1.44, 1.41

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.