- //Code by Warren D. Smith, March 2013. Compiles with gcc.

//Speed, on 2GHz intel gcc 4.2.1, is about 25 minutes to test 2^30 consecutive integers...

//i.e. about 666666 prime-tests per second... which I must say I consider disappointing.

uint hash63778tab[128] = {

//these 50 numbers are the only nonprime odd spsp{2,63778} below 2^32:

979363153, 0, 0, 643552909, 0, 8725753, 0, 0, 0, 0, 329153653, 48191653, 0, 8036033, 0,

3116107, 0, 0, 839280691, 0, 450866021, 0, 0, 0, 0, 0, 0, 2269093, 0, 0, 0, 0, 0, 2693739751,

0, 2717428033, 1464568381, 0, 0, 0, 0, 1714322377, 1140573601, 0, 0, 63318169, 0, 0,

4157008813, 0, 0, 220729, 0, 0, 0, 0, 0, 3207297773, 0, 0, 0, 385319089, 1930534453, 0, 0,

0, 1054999441, 787085857, 0, 0, 2017509601, 0, 0, 967287751, 1272866167, 0, 3548378341,

3274264033, 0, 0, 0, 2206095589, 2766006253, 453366029, 176030977, 126132553, 0, 0, 0, 0,

3242533897, 0, 1121176981, 0, 2181961, 0, 0, 17208601, 0, 0, 0, 0, 48316969, 741751,

724969087, 0, 0, 0, 43661257, 0, 0, 0, 0, 0, 0, 1153164097, 1104194521, 0, 3323829169, 0,

765378241, 22591301, 1894909141, 3900327241, 60547831, 0, 1328256247, 0 };

uint isprime32(uint32 p){ //TRUE iff p is a prime and <=32 bits long.

uint hash;

if(p<4){ if(p>=2){ return(TRUE); } return(FALSE); }

if((p&1)==0){ return(FALSE); }

if( !strongpsp32(p, 2) ){ return(FALSE); }

if(p>=2047){

if( !strongpsp32(p, 63778) ){ return(FALSE); }

hash = ((p*2869299540U)>>25)&127;

if(hash63778tab[hash]==p){ return(FALSE); }

}

return(TRUE);

}

//Code for the strongpsp32 and underlying routines are:

//This can be sped up (no loop needed) on many processors using

//"count trailing zeros" instruction and barrel-shifter.

#define evenize(e,s) for(s=0; (e&1)==0; e/=2){ s++; }

uint32 mulmod32(uint32 a, uint32 b, uint32 c){ //returns d = a*b mod c

uint32 d;

uint64 q;

q = a; q *= b; d = q%c; return(d);

}

//computes g^e mod m.

uint32 modularexp32(uint32 g, uint32 e, uint32 m){

uint32 U, L, q, A, gpow[16];

int j;

uint es;

if(e<=1){ if(e==0){ return(1); } return(g); }

gpow[0] = 1;

gpow[1] = g;

for(j=2; j<8; j++){

gpow[j] = mulmod32(g, gpow[j-1], m);

}

for(j=30; j>=0; j-=3){ //speed up with assembly, count leading zeros??

es = e>>j;

if(es!=0) break;

}

A = gpow[es&7];

while(j>0){

A = mulmod32(A, A, m);

A = mulmod32(A, A, m);

A = mulmod32(A, A, m);

j -= 3;

es = (e>>j)&7;

A = mulmod32(A, gpow[es], m);

}

return(A);

}

uint strongpsp32(uint32 p, uint B){

int s;

uint32 A, pm1, e, U, L, q;

assert(B>1 and (p&1)); //assumes p is odd and p>B>=2...

//p<B can also be ok but you better know what you are doing

pm1 = p-1;

e = pm1;

evenize(e,s);

A = modularexp32((uint32)B, e, p);

if(A==1){ return(TRUE); }

for(;;){

if(A == pm1){ return(TRUE); }

s--;

if(s==0){ return(FALSE); }

A = mulmod32(A,A,p);

}

} - On 3/13/2013 5:47 AM, Phil Carmody wrote:
> No I couldn't. That's so overly verbose and redundant it makes me twitch, I can

Apologies, my Pari/GP script-fu is definitely not on par with DJB's. I just

> barely bring myself to repeat it!

> "lift(Mod(p*s, lift(znorder(Mod(2,s))))) == 1" is just

> "p*s % znorder(Mod(2,s)) == 1"

> Having 3 exit conditions to the loop is overkill too.

wanted to provide a script so that people could plug and chug an r,p pair to see

what psp's would be generated. Also, I didn't know that the % (mod) operator

still worked in Pari. I thought everything had to be done with Mod(). Thanks

for that insight.

> The latter worries me a bit, as it might imply wasted effort. I'm trying to

And apologies here too, I mis-remembered a statement from his Category S page

> picture how these duplicates arise. Given a n, the maximal prime factor p|s is

> uniquely defined, and r as order_2(p) is uniquely defined. Therefore n can only

> appear with pair (r,p)?

and mis-spoke by applying it to the Category E psp's.

On 3/14/2013 11:02 AM, WarrenS wrote:

> 2. Consulting the Cunningham project pages,

> http://homes.cerias.purdue.edu/~ssw/cun/index.html

> every Mersenne-form number 2^r - 1 now is fully factored if

> r<929. Apparently the first two open cases are r=929 and 947

> yielding 214 and 217 digit numbers to factor.

An update here: M929 has been factored, and a group of people have already

started factoring M947. You can find the factor for M929 here:

http://homes.cerias.purdue.edu/~ssw/cun/page125

And, you can see who is factoring which Cunningham number here:

http://homes.cerias.purdue.edu/~ssw/cun/who

And more importantly, you can find all known*1 factors for all important*2

numbers in the online factor database here:

factordb.com

Once there, you can type in 2^929-1, and it will show you all the factors of

that number and that it is Fully Factored (FF). Currently, you can type in

2^947-1 and see that it is a CF, composite number, with factors known, but not

yet fully factored.

*1 = All known factors that have been stored into the factordb.

*2 = All numbers that people are interested in and store in the factordb.

Also, the factordb stores prime numbers too. Below 300 digits it will just

prove the number prime, and above that it will accept Primo certificates and

verify them locally.

-David C.