A lot of composite numbers can be separated from primes by Fermat's Little Theorem, but there are some that completely confound it.

The   Miller Rabin Test   uses a combination of Fermat's Little Theorem and Chinese Division Theorem to overcome this.

The purpose of this task is to investigate such numbers using a method based on   Carmichael numbers,   as suggested in   Notes by G.J.O Jameson March 2010.



Task

Find Carmichael numbers of the form:

  Prime1 × Prime2 × Prime3 
where    (Prime1 < Prime2 < Prime3)    for all    Prime1    up to   61.

(See page 7 of   Notes by G.J.O Jameson March 2010   for solutions.)



Pseudocode

For a given   



P
r
i
m

e

1




{\displaystyle Prime_{1}}



for 1 < h3 < Prime1
    for 0 < d < h3+Prime1
         if (h3+Prime1)*(Prime1-1) mod d == 0 and -Prime1 squared mod h3 == d mod h3
         then
               Prime2 = 1 + ((Prime1-1) * (h3+Prime1)/d)
               next d if Prime2 is not prime
               Prime3 = 1 + (Prime1*Prime2/h3)
               next d if Prime3 is not prime
               next d if (Prime2*Prime3) mod (Prime1-1) not equal 1
               Prime1 * Prime2 * Prime3 is a Carmichael Number



related task
Chernick's Carmichael numbers

# sieve of Eratosthene: sets s[i] to TRUE if i is prime, FALSE otherwise #
PROC sieve = ( REF[]BOOL s )VOID:
     BEGIN
        # start with everything flagged as prime                             # 
        FOR i TO UPB s DO s[ i ] := TRUE OD;
        # sieve out the non-primes                                           #
        s[ 1 ] := FALSE;
        FOR i FROM 2 TO ENTIER sqrt( UPB s ) DO
            IF s[ i ] THEN FOR p FROM i * i BY i TO UPB s DO s[ p ] := FALSE OD FI
        OD
     END # sieve # ;
 
# construct a sieve of primes up to the maximum number required for the task #
# For Prime1, we need to check numbers up to around 120 000                  #
INT max number = 200 000;
[ 1 : max number ]BOOL is prime;
sieve( is prime );
 
# Find the Carmichael 3 Stromg Pseudoprimes for Prime1 up to 61              #
 
FOR prime1 FROM 2 TO 61 DO
    IF is prime[ prime 1 ] THEN
        FOR h3 TO prime1 - 1 DO
            FOR d TO ( h3 + prime1 ) - 1 DO
                IF   ( h3 + prime1 ) * ( prime1 - 1 ) MOD d = 0
                AND ( - ( prime1 * prime1 ) ) MOD h3 = d MOD h3
                THEN
                    INT prime2 = 1 + ( ( prime1 - 1 ) * ( h3 + prime1 ) OVER d );
                    IF is prime[ prime2 ] THEN
                        INT prime3 = 1 + ( prime1 * prime2 OVER h3 );
                        IF is prime[ prime3 ] THEN 
                            IF ( prime2 * prime3 ) MOD ( prime1 - 1 ) = 1 THEN
                                print( ( whole( prime1, 0 ), " ", whole( prime2, 0 ), " ", whole( prime3, 0 ), newline ) )
                            FI
                        FI
                    FI
                FI
            OD
        OD
    FI
OD