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 "seed7_05.s7i";
 
const func boolean: isPrime (in integer: number) is func
  result
    var boolean: prime is FALSE;
  local
    var integer: upTo is 0;
    var integer: testNum is 3;
  begin
    if number = 2 then
      prime := TRUE;
    elsif odd(number) and number > 2 then
      upTo := sqrt(number);
      while number rem testNum <> 0 and testNum <= upTo do
        testNum +:= 2;
      end while;
      prime := testNum > upTo;
    end if;
  end func;
 
const proc: main is func
  local
    var integer: p1 is 0;
    var integer: h3 is 0;
    var integer: g is 0;
    var integer: d is 0;
    var integer: p2 is 0;
    var integer: p3 is 0;
  begin
    for p1 range 2 to 61 do
      if isPrime(p1) then
        for h3 range 2 to p1 do
          g := h3 + p1;
          for d range 1 to pred(g) do
            if (g * pred(p1)) mod d = 0 and -p1 ** 2 mod h3 = d mod h3 then
              p2 := 1 + pred(p1) * g div d;
              if isPrime(p2) then
                p3 := 1 + p1 * p2 div h3;
                if isPrime(p3) and (p2 * p3) mod pred(p1) = 1 then
                  writeln(p1 <& " * " <& p2 <& " * " <& p3 <& " = " <& p1*p2*p3);
                end if;
              end if;
            end if;
          end for;
        end for;
      end if;
    end for;
  end func;