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

import Foundation
 
extension BinaryInteger {
  @inlinable
  public var isPrime: Bool {
    if self == 0 || self == 1 {
      return false
    } else if self == 2 {
      return true
    }
 
    let max = Self(ceil((Double(self).squareRoot())))
 
    for i in stride(from: 2, through: max, by: 1) {
      if self % i == 0 {
        return false
      }
    }
 
    return true
  }
}
 
@inlinable
public func carmichael<T: BinaryInteger & SignedNumeric>(p1: T) -> [(T, T, T)] {
  func mod(_ n: T, _ m: T) -> T { (n % m + m) % m }
 
  var res = [(T, T, T)]()
 
  guard p1.isPrime else {
    return res
  }
 
  for h3 in stride(from: 2, to: p1, by: 1) {
    for d in stride(from: 1, to: h3 + p1, by: 1) {
      if (h3 + p1) * (p1 - 1) % d != 0 || mod(-p1 * p1, h3) != d % h3 {
        continue
      }
 
      let p2 = 1 + (p1 - 1) * (h3 + p1) / d
 
      guard p2.isPrime else {
        continue
      }
 
      let p3 = 1 + p1 * p2 / h3
 
      guard p3.isPrime && (p2 * p3) % (p1 - 1) == 1 else {
        continue
      }
 
      res.append((p1, p2, p3))
    }
  }
 
  return res
}
 
 
let res =
  (1..<62)
    .lazy
    .filter({ $0.isPrime })
    .map(carmichael)
    .filter({ !$0.isEmpty })
    .flatMap({ $0 })
 
for c in res {
  print(c)
}