(*
   written 2022/12/13
   ported  2025/12/17

   https://en.wikipedia.org/wiki/Carmichael_number

   According to Korselt's criterion, a number N is a Carmichael number
   if and only if it satisfies three properties.  First, it must have
   more than one prime factor.  Second, no prime factor can repeat.
   And third, for every prime p that divides N, p – 1 also divides N – 1.

   Consider again the number 561.  It's equal to 3 × 11 × 17, so it
   clearly satisfies the first two properties in Korselt's list.  To
   show the last property, subtract 1 from each prime factor to get
   2, 10 and 16.  In addition, subtract 1 from 561.  All three of the
   smaller numbers are divisors of 560.  The number 561 is therefore a
   Carmichael number.

   Carmichael numbers are relatively rare instances that violate
   Fermat's Little Theorem, which states that if b**n - b is divisible
   by n for all values of b, then n is prime.

   My initial algorithm:
      walk up the primes to half(n)
         - find a prime divisor p:
         - test if it is a repeat   (if yes: fail)
         - test if p-1 divides n-1  (if no : fail)
       return true if number of prime divisors > 1
*)

(* -------------------------------------- *
 * code from the Klein library            *
 * -------------------------------------- *)

function MOD(m: integer, n: integer): integer
    m - m/n * n

(* -------------------------------------- *
 * code to test if prime, from sieve      *
 * -------------------------------------- *)

function isPrime(n: integer): boolean
    not hasDivisorFrom(2, n)

function hasDivisorFrom(i: integer, n: integer): boolean
    if i < n then
       divides(i, n) or hasDivisorFrom(i+1, n)
    else
       false

function divides(a: integer, b: integer): boolean
    MOD(b, a) = 0

(* -------------------------------------- *
 * for finding next prime to consider     *
 * -------------------------------------- *)

function nextp(n: integer): integer
    if (isPrime(n))
       then n
       else nextp(n+1)

(* -------------------------------------- *
 * main loop: implement Korselt's def     *
 * -------------------------------------- *)

function examine_prime_divs(currn: integer, currp: integer,
                            lastpdiv: integer, numpdivs: integer,
                            n: integer)
    : boolean
    if (n/2 < currp) then                        (* done... *)
       1 < numpdivs                              (*   are there multiple p? *)
    else if (not divides(currp, currn)) then     (* not a divisor, so   *)
       examine_prime_divs(currn, nextp(currp+1), (*   go to next prime  *)
                          lastpdiv, numpdivs, n)
    else if (currp = lastpdiv) then              (* p is a repeat *)
       false
    else if (not divides(currp-1, n-1)) then     (* p-1 doesn't divide n-1 *)
       false
    else                                         (* p is a new prime divisor *)
       examine_prime_divs(currn/currp, currp,
                          currp, numpdivs+1, n)

(* -------------------------------------- *
 * main launches the loop                 *
 * -------------------------------------- *)

function main(n : integer): boolean
    examine_prime_divs(n, 2, -1, 0, n)

(* --------------------------------------*
 * this Python function could be adapted *
 * to generate Carmichael numbers <= n   *
 * --------------------------------------*

   # standard limit prevents testing up to 1000
   # this depth works up to at least n == 3000

   import sys
   sys.setrecursionlimit(5000)

   def test_up_to(max):
       for i in range(1, max+1):
           if main(i):
               print(i)
*)