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

#include <iomanip>
#include <iostream>
 
int mod(int n, int d) {
    return (d + n % d) % d;
}
 
bool is_prime(int n) {
    if (n < 2)
        return false;
    if (n % 2 == 0)
        return n == 2;
    if (n % 3 == 0)
        return n == 3;
    for (int p = 5; p * p <= n; p += 4) {
        if (n % p == 0)
            return false;
        p += 2;
        if (n % p == 0)
            return false;
    }
    return true;
}
 
void print_carmichael_numbers(int prime1) {
    for (int h3 = 1; h3 < prime1; ++h3) {
        for (int d = 1; d < h3 + prime1; ++d) {
            if (mod((h3 + prime1) * (prime1 - 1), d) != 0
                || mod(-prime1 * prime1, h3) != mod(d, h3))
                continue;
            int prime2 = 1 + (prime1 - 1) * (h3 + prime1)/d;
            if (!is_prime(prime2))
                continue;
            int prime3 = 1 + prime1 * prime2/h3;
            if (!is_prime(prime3))
                continue;
            if (mod(prime2 * prime3, prime1 - 1) != 1)
                continue;
            unsigned int c = prime1 * prime2 * prime3;
            std::cout << std::setw(2) << prime1 << " x "
                << std::setw(4) << prime2 << " x "
                << std::setw(5) << prime3 << " = "
                << std::setw(10) << c << '\n';
        }
    }
}
 
int main() {
    for (int p = 2; p <= 61; ++p) {
        if (is_prime(p))
            print_carmichael_numbers(p);
    }
    return 0;
}