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

 
(ns example
  (:gen-class))
 
(defn prime? [n]
  " Prime number test (using Java) "
  (.isProbablePrime (biginteger n) 16))
 
(defn carmichael [p1]
  " Triplets of Carmichael primes, with first element prime p1 "
  (if (prime? p1)
    (into [] (for [h3 (range 2 p1)
          :let [g (+ h3 p1)]
          d (range 1 g)
          :when (and (= (mod (* g (dec p1)) d) 0)
                     (= (mod (- (* p1 p1)) h3) (mod d h3)))
          :let [p2 (inc (quot (* (dec p1) g) d))]
          :when (prime? p2)
          :let [p3 (inc (quot (* p1 p2) h3))]
          :when (prime? p3)
          :when (= (mod (* p2 p3) (dec p1)) 1)]
         [p1 p2 p3]))))
 
; Generate Result
(def numbers (mapcat carmichael (range 2 62)))
(println (count numbers) "Carmichael numbers found:")
(doseq [t numbers]
  (println (format "%5d x %5d x %5d = %10d" (first t) (second t) (last t) (apply * t))))
 