On generating Strong PseudoPrimes DataBase?

[TheoryQuest1]

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 ?

[LaurV]

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:

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.

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).

/*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:

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:

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 ...

[danaj]

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.

[xilman]

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.

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.

[danaj]

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!

[xilman]

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.

[danaj]

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.

[xilman]

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

[danaj]

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).

[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.

[danaj]

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.