Pseudoprimes and Strong Pseudoprimes

Pseudoprimes

The tests for primality using pseudoprimes and strong  pseudoprimes can be easily implemented
    in Mathematica. Here are some functions which can be used to find pseudoprimes, strong pseodoprimes and implement a probabilisitic primality test.   First, to check if a number is a  pseudoprime to some base, we can define

pspQ[a_,n_]:= PowerMod[a,n-1,n] ==1 && !PrimeQ[n]

Note that we are using the built-in primeQ to  test if the number is composite.  To find the smallest pseudoprime to base 5 which is larger than 5, we can type

n=5;
 While[ !pspQ[5,n], n++];
 Print[n]

We verify  if  124 is a pseudoprime to base 5.

PowerMod[5,123,124]

The code can be modified to find pseudoprimes with different bases and conditions.  Regarding Carmichael numbers, it is not practical to verify the definition for a large number. Instead it is better to use the characterization of a Carmichael number in terms of its prime factorization. Naturally, this makes the discovery of large Carmichael numbers more difficult.

Strong Pseudoprimes

The  strong pseudoprime test for n  requires determining the even part of n-1. 
 We want to use this test to find large probable primes, hence we  should not
 incorporate the built-in primality test into the strong pseudoprime test. To 
eliminate the primes, one should use a program for trial division, which is suitable
 only up to some limit like 10^{12}.  Here is a program which tests if a
 number n  passes the strong pseudoprime test to base a. The program is designed to print the intermediate  computations. You may delete the print statements if you like. 
 
  First, we have to factor  n-1= 2^r s. This is the prupose of the first  While loop.
  The program then prints s and r.  Then we compute the terms b_0= a^s \pmod{n}
 and its squares until we reach the last value,  a^{n-1}\pmod{n} or we get 1. 
 We first check if  a^s\equiv 1\pmod{n}. If this is the case, then no further  information can be obtained from the computation and we declare the number to a porbable strong pseudoprime. If this is not the case, the While  loop computes succesive squares terminates    if one of values is 
1  or -1  or we have computed all the values. In the last case, if i=r
 and we haven't encountered 1 or -1 before then the number must be composite. 
If  i<r  and  b  is 1, then the number must be composite, because the previous 
value was not +1 or -1, Otherwise, the algorithm would have stopped at that 
value.

probableSpsp[a_,n_]:= Module[ {s=n-1,r=0,i=0,b,
                     spsp=True},
        
             While[Mod[s,2]==0, s=s/2;r++];    
             Print["s=",s,"\t", " and r=", r];
         
             b=PowerMod[a,s,n];
             Print[i, "\t", b]; 
             
             If[b!=1,
               While[!(b ==1 || b==n-1)&& i<r ,
                  b=PowerMod[b,2,n];i++;
                   Print[i,"\t",b];
                   ];
             If[i==r , spsp=False,
                If[b!=n-1, spsp=False]; 
               ];
             ];
                 
                     
           Return[spsp]]

We test it on a few numbers. Naturally, the test  will return true for all primes.

probableSpsp[2,1037]

probableSpsp[2,15790321]

probableSpsp[3,15790321]

probableSpsp[5,15790321]

The above computation shows that $15790321$ is actually prime as it is not one of the thirteen strong pseudoprimes to base $2$,  $3$ and $5$ which are   less than $25 \times 10^9$.

 We test in a  few more cases to verify the correctness of the procedure.

probableSpsp[2,911 1037]

probableSpsp[2, 2047]

probableSpsp[3,2047]

We can combine the test for several bases by making the function {\bf Listable}

SetAttributes[probableSpsp, Listable]

probableSpsp[{2,3,5,7}, 960946321]

The number reveals its compositeness  in the test to base $7$.