# Carmichael 3 strong pseudoprimes

Carmichael 3 strong pseudoprimes
You are encouraged to solve this task according to the task description, using any language you may know.

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.

The objective is to find Carmichael numbers of the form $Prime_1 \times Prime_2 \times Prime_3$ (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 Prime1

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

## Contents

Uses the Miller_Rabin package from Miller-Rabin primality test#ordinary integers.

with Ada.Text_IO, Miller_Rabin; procedure Nemesis is    type Number is range 0 .. 2**40-1; -- sufficiently large for the task    function Is_Prime(N: Number) return Boolean is      package MR is new Miller_Rabin(Number); use MR;   begin      return MR.Is_Prime(N) = Probably_Prime;   end Is_Prime; begin   for P1 in Number(2) .. 61 loop      if Is_Prime(P1) then         for H3 in Number(1) .. P1 loop            declare               G: Number := H3 + P1;               P2, P3: Number;            begin               Inner:               for D in 1 .. G-1 loop                  if ((H3+P1) * (P1-1)) mod D = 0 and then                    (-(P1 * P1)) mod H3 = D mod H3                  then                     P2 := 1 + ((P1-1) * G / D);                     P3 := 1 +(P1*P2/H3);                     if Is_Prime(P2) and then Is_Prime(P3)                       and then (P2*P3) mod (P1-1) = 1                     then                       Ada.Text_IO.Put_Line                        ( Number'Image(P1) & " *"   & Number'Image(P2) & " *" &                          Number'Image(P3) & "  = " & Number'Image(P1*P2*P3) );                     end if;                  end if;               end loop Inner;            end;         end loop;      end if;   end loop;end Nemesis;
Output:
 3 * 11 * 17  =  561
5 * 29 * 73  =  10585
5 * 17 * 29  =  2465
5 * 13 * 17  =  1105
7 * 19 * 67  =  8911

... (the full output is 69 lines long) ...

61 * 271 * 571  =  9439201
61 * 241 * 421  =  6189121
61 * 3361 * 4021  =  824389441

## C

 #include <stdio.h> /* C's % operator actually calculates the remainder of a / b so we need a * small adjustment so it works as expected for negative values */#define mod(n,m) ((((n) % (m)) + (m)) % (m)) int is_prime(unsigned int n){    if (n <= 3) {        return n > 1;    }    else if (!(n % 2) || !(n % 3)) {        return 0;    }    else {        unsigned int i;        for (i = 5; i*i <= n; i += 6)            if (!(n % i) || !(n % (i + 2)))                return 0;        return 1;    }} void carmichael3(int p1){    if (!is_prime(p1)) return;     int h3, d, p2, p3;    for (h3 = 1; h3 < p1; ++h3) {        for (d = 1; d < h3 + p1; ++d) {            if ((h3 + p1)*(p1 - 1) % d == 0 && mod(-p1 * p1, h3) == d % h3) {                p2 = 1 + ((p1 - 1) * (h3 + p1)/d);                if (!is_prime(p2)) continue;                p3 = 1 + (p1 * p2 / h3);                if (!is_prime(p3) || (p2 * p3) % (p1 - 1) != 1) continue;                printf("%d %d %d\n", p1, p2, p3);            }        }    }} int main(void){    int p1;    for (p1 = 2; p1 < 62; ++p1)        carmichael3(p1);    return 0;}
Output:
3 11 17
5 29 73
5 17 29
5 13 17
7 19 67
7 31 73
.
.
.
61 181 1381
61 541 3001
61 661 2521
61 271 571
61 241 421
61 3361 4021


## D

enum mod = (in int n, in int m) pure nothrow @nogc=> ((n % m) + m) % m; bool isPrime(in uint n) pure nothrow @nogc {  if (n == 2 || n == 3)    return true;  else if (n < 2 || n % 2 == 0 || n % 3 == 0)    return false;  for (uint div = 5, inc = 2; div ^^ 2 <= n;     div += inc, inc = 6 - inc)    if (n % div == 0)      return false;  return true;} void main() {  import std.stdio;   foreach (immutable p; 2 .. 62) {    if (!p.isPrime) continue;    foreach (immutable h3; 2 .. p) {      immutable g = h3 + p;      foreach (immutable d; 1 .. g) {        if ((g * (p - 1)) % d != 0 || mod(-p * p, h3) != d % h3)          continue;        immutable q = 1 + (p - 1) * g / d;        if (!q.isPrime) continue;        immutable r = 1 + (p * q / h3);        if (!r.isPrime || (q * r) % (p - 1) != 1) continue;        writeln(p, " x ", q, " x ", r);      }    }  }}
Output:
3 x 11 x 17
5 x 29 x 73
5 x 17 x 29
5 x 13 x 17
7 x 19 x 67
7 x 31 x 73
7 x 13 x 31
7 x 23 x 41
7 x 73 x 103
7 x 13 x 19
13 x 61 x 397
13 x 37 x 241
13 x 97 x 421
13 x 37 x 97
13 x 37 x 61
17 x 41 x 233
17 x 353 x 1201
19 x 43 x 409
19 x 199 x 271
23 x 199 x 353
29 x 113 x 1093
29 x 197 x 953
31 x 991 x 15361
31 x 61 x 631
31 x 151 x 1171
31 x 61 x 271
31 x 61 x 211
31 x 271 x 601
31 x 181 x 331
37 x 109 x 2017
37 x 73 x 541
37 x 613 x 1621
37 x 73 x 181
37 x 73 x 109
41 x 1721 x 35281
41 x 881 x 12041
41 x 101 x 461
41 x 241 x 761
41 x 241 x 521
41 x 73 x 137
41 x 61 x 101
43 x 631 x 13567
43 x 271 x 5827
43 x 127 x 2731
43 x 127 x 1093
43 x 211 x 757
43 x 631 x 1597
43 x 127 x 211
43 x 211 x 337
43 x 433 x 643
43 x 547 x 673
43 x 3361 x 3907
47 x 3359 x 6073
47 x 1151 x 1933
47 x 3727 x 5153
53 x 157 x 2081
53 x 79 x 599
53 x 157 x 521
59 x 1451 x 2089
61 x 421 x 12841
61 x 181 x 5521
61 x 1301 x 19841
61 x 277 x 2113
61 x 181 x 1381
61 x 541 x 3001
61 x 661 x 2521
61 x 271 x 571
61 x 241 x 421
61 x 3361 x 4021

## EchoLisp

 ;; charmichaΓ«l numbers up to N-th prime ; 61 is 18-th prime(define (charms (N 18) local: (h31 0) (Prime2 0) (Prime3 0))(for* ((Prime1 (primes N))       (h3 (in-range 1 Prime1))       (d  (+ h3 Prime1)))      (set! h31 (+ h3 Prime1))      #:continue (!zero? (modulo (* h31 (1- Prime1)) d))      #:continue (!= (modulo d h3) (modulo (- (* Prime1 Prime1)) h3))      (set! Prime2 (1+ ( * (1- Prime1) (quotient h31 d))))      #:when (prime? Prime2)      (set! Prime3 (1+ (quotient (*  Prime1  Prime2)  h3)))      #:when (prime? Prime3)      #:when (= 1 (modulo (* Prime2 Prime3) (1- Prime1)))      (printf " π₯ %12d = %d x %d x %d"  (* Prime1 Prime2 Prime3) Prime1 Prime2 Prime3)))
Output:
 (charms 3)π₯          561 = 3 x 11 x 17π₯        10585 = 5 x 29 x 73π₯         2465 = 5 x 17 x 29π₯         1105 = 5 x 13 x 17 (charms 18);; skipped ....π₯    902645857 = 47 x 3727 x 5153π₯      2632033 = 53 x 53 x 937π₯     17316001 = 53 x 157 x 2081π₯      4335241 = 53 x 157 x 521π₯    178837201 = 59 x 1451 x 2089π₯    329769721 = 61 x 421 x 12841π₯     60957361 = 61 x 181 x 5521π₯      6924781 = 61 x 61 x 1861π₯      6924781 = 61 x 61 x 1861π₯     15247621 = 61 x 181 x 1381π₯     99036001 = 61 x 541 x 3001π₯    101649241 = 61 x 661 x 2521π₯      6189121 = 61 x 241 x 421π₯    824389441 = 61 x 3361 x 4021

Translation of: Ruby
Library: primes
Works with: GHC version 7.4.1
Works with: primes version 0.2.1.0
#!/usr/bin/runhaskell import Data.Numbers.Primesimport Control.Monad (guard) carmichaels = do  p <- takeWhile (<= 61) primes  h3 <- [2..(p-1)]  let g = h3 + p  d <- [1..(g-1)]  guard $(g * (p - 1)) mod d == 0 && (-1 * p * p) mod h3 == d mod h3 let q = 1 + (((p - 1) * g) div d) guard$ isPrime q  let r = 1 + ((p * q) div h3)  guard $isPrime r && (q * r) mod (p - 1) == 1 return (p, q, r) main = putStr$ unlines $map show carmichaels Output: (3,11,17) (5,29,73) (5,17,29) (5,13,17) (7,19,67) (7,31,73) (7,13,31) (7,23,41) (7,73,103) (7,13,19) (13,61,397) (13,37,241) (13,97,421) (13,37,97) (13,37,61) (17,41,233) (17,353,1201) (19,43,409) (19,199,271) (23,199,353) (29,113,1093) (29,197,953) (31,991,15361) (31,61,631) (31,151,1171) (31,61,271) (31,61,211) (31,271,601) (31,181,331) (37,109,2017) (37,73,541) (37,613,1621) (37,73,181) (37,73,109) (41,1721,35281) (41,881,12041) (41,101,461) (41,241,761) (41,241,521) (41,73,137) (41,61,101) (43,631,13567) (43,271,5827) (43,127,2731) (43,127,1093) (43,211,757) (43,631,1597) (43,127,211) (43,211,337) (43,433,643) (43,547,673) (43,3361,3907) (47,3359,6073) (47,1151,1933) (47,3727,5153) (53,157,2081) (53,79,599) (53,157,521) (59,1451,2089) (61,421,12841) (61,181,5521) (61,1301,19841) (61,277,2113) (61,181,1381) (61,541,3001) (61,661,2521) (61,271,571) (61,241,421) (61,3361,4021)  ## Icon and Unicon The following works in both languages. link "factors" procedure main(A) n := integer(!A) | 61 every write(carmichael3(!n))end procedure carmichael3(p1) every (isprime(p1), (h := 1+!(p1-1)), (d := !(h+p1-1))) do if (mod(((h+p1)*(p1-1)),d) = 0, mod((-p1*p1),h) = mod(d,h)) then { p2 := 1 + (p1-1)*(h+p1)/d p3 := 1 + p1*p2/h if (isprime(p2), isprime(p3), mod((p2*p3),(p1-1)) = 1) then suspend format(p1,p2,p3) }end procedure mod(n,d) return (d+n%d)%dend procedure format(p1,p2,p3) return left(p1||" * "||p2||" * "||p3,20)||" = "||(p1*p2*p3)end Output, with middle lines elided: ->c3sp 3 * 11 * 17 = 561 5 * 29 * 73 = 10585 5 * 17 * 29 = 2465 5 * 13 * 17 = 1105 7 * 19 * 67 = 8911 7 * 31 * 73 = 15841 7 * 13 * 31 = 2821 7 * 23 * 41 = 6601 7 * 73 * 103 = 52633 7 * 13 * 19 = 1729 13 * 61 * 397 = 314821 13 * 37 * 241 = 115921 ... 53 * 157 * 2081 = 17316001 53 * 79 * 599 = 2508013 53 * 157 * 521 = 4335241 59 * 1451 * 2089 = 178837201 61 * 421 * 12841 = 329769721 61 * 181 * 5521 = 60957361 61 * 1301 * 19841 = 1574601601 61 * 277 * 2113 = 35703361 61 * 181 * 1381 = 15247621 61 * 541 * 3001 = 99036001 61 * 661 * 2521 = 101649241 61 * 271 * 571 = 9439201 61 * 241 * 421 = 6189121 61 * 3361 * 4021 = 824389441 ->  ## Julia This solution is a straightforward implementation of the algorithm of the Jameson paper cited in the task description. Just for fun, I use Julia's capacity to accommodate Unicode identifiers to match some of the paper's symbols to the variables used in the carmichael function. Function  function carmichael{T<:Integer}(pmax::T) 0 < pmax || throw(DomainError()) car = T[] for p in primes(pmax) for hβ in 2:(p-1) m = (p - 1)*(hβ + p) pmh = mod(-p^2, hβ) for Ξ in 1:(hβ+p-1) m%Ξ==0 && Ξ%hβ==pmh || continue q = div(m, Ξ) + 1 isprime(q) || continue r = div((p*q - 1), hβ) + 1 isprime(r) && mod(q*r, (p-1))==1 || continue append!(car, [p, q, r]) end end end reshape(car, 3, div(length(car), 3))end  Main  hi = 61car = carmichael(hi) curp = 0tcnt = 0print("Carmichael 3 (p\u00d7q\u00d7r) Pseudoprimes, up to p = ", hi, ":")for j in sortperm(1:size(car)[2], by=x->(car[1,x], car[2,x], car[3,x])) p, q, r = car[:,j] c = prod(car[:,j]) if p != curp curp = p print(@sprintf("\n\np = %d\n ", p)) tcnt = 0 end if tcnt == 4 print("\n ") tcnt = 1 else tcnt += 1 end print(@sprintf("p\u00d7%d\u00d7%d = %d ", q, r, c))endprintln("\n\n", size(car)[2], " results in total.")  Output: Carmichael 3 (pΓqΓr) Pseudoprimes, up to p = 61: p = 3 pΓ11Γ17 = 561 p = 5 pΓ13Γ17 = 1105 pΓ17Γ29 = 2465 pΓ29Γ73 = 10585 p = 7 pΓ13Γ19 = 1729 pΓ13Γ31 = 2821 pΓ19Γ67 = 8911 pΓ23Γ41 = 6601 pΓ31Γ73 = 15841 pΓ73Γ103 = 52633 p = 13 pΓ37Γ61 = 29341 pΓ37Γ97 = 46657 pΓ37Γ241 = 115921 pΓ61Γ397 = 314821 pΓ97Γ421 = 530881 p = 17 pΓ41Γ233 = 162401 pΓ353Γ1201 = 7207201 p = 19 pΓ43Γ409 = 334153 pΓ199Γ271 = 1024651 p = 23 pΓ199Γ353 = 1615681 p = 29 pΓ113Γ1093 = 3581761 pΓ197Γ953 = 5444489 p = 31 pΓ61Γ211 = 399001 pΓ61Γ271 = 512461 pΓ61Γ631 = 1193221 pΓ151Γ1171 = 5481451 pΓ181Γ331 = 1857241 pΓ271Γ601 = 5049001 pΓ991Γ15361 = 471905281 p = 37 pΓ73Γ109 = 294409 pΓ73Γ181 = 488881 pΓ73Γ541 = 1461241 pΓ109Γ2017 = 8134561 pΓ613Γ1621 = 36765901 p = 41 pΓ61Γ101 = 252601 pΓ73Γ137 = 410041 pΓ101Γ461 = 1909001 pΓ241Γ521 = 5148001 pΓ241Γ761 = 7519441 pΓ881Γ12041 = 434932961 pΓ1721Γ35281 = 2489462641 p = 43 pΓ127Γ211 = 1152271 pΓ127Γ1093 = 5968873 pΓ127Γ2731 = 14913991 pΓ211Γ337 = 3057601 pΓ211Γ757 = 6868261 pΓ271Γ5827 = 67902031 pΓ433Γ643 = 11972017 pΓ547Γ673 = 15829633 pΓ631Γ1597 = 43331401 pΓ631Γ13567 = 368113411 pΓ3361Γ3907 = 564651361 p = 47 pΓ1151Γ1933 = 104569501 pΓ3359Γ6073 = 958762729 pΓ3727Γ5153 = 902645857 p = 53 pΓ79Γ599 = 2508013 pΓ157Γ521 = 4335241 pΓ157Γ2081 = 17316001 p = 59 pΓ1451Γ2089 = 178837201 p = 61 pΓ181Γ1381 = 15247621 pΓ181Γ5521 = 60957361 pΓ241Γ421 = 6189121 pΓ271Γ571 = 9439201 pΓ277Γ2113 = 35703361 pΓ421Γ12841 = 329769721 pΓ541Γ3001 = 99036001 pΓ661Γ2521 = 101649241 pΓ1301Γ19841 = 1574601601 pΓ3361Γ4021 = 824389441 69 results in total.  ## Mathematica / Wolfram Language Cases[Cases[ Cases[Table[{p1, h3, d}, {p1, Array[Prime, PrimePi@61]}, {h3, 2, p1 - 1}, {d, 1, h3 + p1 - 1}], {p1_Integer, h3_, d_} /; PrimeQ[1 + (p1 - 1) (h3 + p1)/d] && Divisible[p1^2 + d, h3] :> {p1, 1 + (p1 - 1) (h3 + p1)/d, h3}, Infinity], {p1_, p2_, h3_} /; PrimeQ[1 + Floor[p1 p2/h3]] :> {p1, p2, 1 + Floor[p1 p2/h3]}], {p1_, p2_, p3_} /; Mod[p2 p3, p1 - 1] == 1 :> Print[p1, "*", p2, "*", p3]] ## PARI/GP f(p)={ my(v=List(),q,r); for(h=2,p-1, for(d=1,h+p-1, if((h+p)*(p-1)%d==0 && Mod(p,h)^2==-d && isprime(q=(p-1)*(h+p)/d+1) && isprime(r=p*q\h+1)&&q*r%(p-1)==1, listput(v,p*q*r) ) ) ); Set(v)};forprime(p=3,67,v=f(p); for(i=1,#v,print1(v[i]", "))) Output: 561, 1105, 2465, 10585, 1729, 2821, 6601, 8911, 15841, 52633, 29341, 46657, 115921, 314821, 530881, 162401, 7207201, 334153, 1024651, 1615681, 3581761, 5444489, 399001, 512461, 1193221, 1857241, 5049001, 5481451, 471905281, 294409, 488881, 1461241, 8134561, 36765901, 252601, 410041, 1909001, 5148001, 7519441, 434932961, 2489462641, 1152271, 3057601, 5968873, 6868261, 11972017, 14913991, 15829633, 43331401, 67902031, 368113411, 564651361, 104569501, 902645857, 958762729, 2508013, 4335241, 17316001, 178837201, 6189121, 9439201, 15247621, 35703361, 60957361, 99036001, 101649241, 329769721, 824389441, 1574601601, 10267951, 163954561, 7991602081, ## Perl Library: ntheory use ntheory qw/forprimes is_prime vecprod/; forprimes { my$p = $_; for my$h3 (2 .. $p-1) { my$ph3 = $p +$h3;      for my $d (1 ..$ph3-1) {               # Jameseon procedure page 6         next if ((-$p*$p) % $h3) != ($d % $h3); next if (($p-1)*$ph3) %$d;         my $q = 1 + ($p-1)*$ph3 /$d;        # Jameson eq 7         next unless is_prime($q); my$r = 1 + ($p*$q-1) / $h3; # Jameson eq 6 next unless is_prime($r);         next unless ($q*$r) % ($p-1) == 1; printf "%2d x %5d x %5d = %s\n",$p,$q,$r,vecprod($p,$q,$r); } }} 3,61; Output:  3 x 11 x 17 = 561 5 x 29 x 73 = 10585 5 x 17 x 29 = 2465 5 x 13 x 17 = 1105 ... full output is 69 lines ... 61 x 661 x 2521 = 101649241 61 x 271 x 571 = 9439201 61 x 241 x 421 = 6189121 61 x 3361 x 4021 = 824389441  ## Perl 6 An almost direct translation of the pseudocode. We take the liberty of going up to 67 to show we aren't limited to 32-bit integers. (Perl 6 uses arbitrary precision in any case.) for (2..67).grep: *.is-prime -> \Prime1 { for 1 ^..^ Prime1 -> \h3 { my \g = h3 + Prime1; for 0 ^..^ h3 + Prime1 -> \d { if (h3 + Prime1) * (Prime1 - 1) %% d and -Prime1**2 % h3 == d % h3 { my \Prime2 = floor 1 + (Prime1 - 1) * g / d; next unless Prime2.is-prime; my \Prime3 = floor 1 + Prime1 * Prime2 / h3; next unless Prime3.is-prime; next unless (Prime2 * Prime3) % (Prime1 - 1) == 1; say "{Prime1} Γ {Prime2} Γ {Prime3} == {Prime1 * Prime2 * Prime3}"; } } }} Output: 3 Γ 11 Γ 17 == 561 5 Γ 29 Γ 73 == 10585 5 Γ 17 Γ 29 == 2465 5 Γ 13 Γ 17 == 1105 7 Γ 19 Γ 67 == 8911 7 Γ 31 Γ 73 == 15841 7 Γ 13 Γ 31 == 2821 7 Γ 23 Γ 41 == 6601 7 Γ 73 Γ 103 == 52633 7 Γ 13 Γ 19 == 1729 13 Γ 61 Γ 397 == 314821 13 Γ 37 Γ 241 == 115921 13 Γ 97 Γ 421 == 530881 13 Γ 37 Γ 97 == 46657 13 Γ 37 Γ 61 == 29341 17 Γ 41 Γ 233 == 162401 17 Γ 353 Γ 1201 == 7207201 19 Γ 43 Γ 409 == 334153 19 Γ 199 Γ 271 == 1024651 23 Γ 199 Γ 353 == 1615681 29 Γ 113 Γ 1093 == 3581761 29 Γ 197 Γ 953 == 5444489 31 Γ 991 Γ 15361 == 471905281 31 Γ 61 Γ 631 == 1193221 31 Γ 151 Γ 1171 == 5481451 31 Γ 61 Γ 271 == 512461 31 Γ 61 Γ 211 == 399001 31 Γ 271 Γ 601 == 5049001 31 Γ 181 Γ 331 == 1857241 37 Γ 109 Γ 2017 == 8134561 37 Γ 73 Γ 541 == 1461241 37 Γ 613 Γ 1621 == 36765901 37 Γ 73 Γ 181 == 488881 37 Γ 73 Γ 109 == 294409 41 Γ 1721 Γ 35281 == 2489462641 41 Γ 881 Γ 12041 == 434932961 41 Γ 101 Γ 461 == 1909001 41 Γ 241 Γ 761 == 7519441 41 Γ 241 Γ 521 == 5148001 41 Γ 73 Γ 137 == 410041 41 Γ 61 Γ 101 == 252601 43 Γ 631 Γ 13567 == 368113411 43 Γ 271 Γ 5827 == 67902031 43 Γ 127 Γ 2731 == 14913991 43 Γ 127 Γ 1093 == 5968873 43 Γ 211 Γ 757 == 6868261 43 Γ 631 Γ 1597 == 43331401 43 Γ 127 Γ 211 == 1152271 43 Γ 211 Γ 337 == 3057601 43 Γ 433 Γ 643 == 11972017 43 Γ 547 Γ 673 == 15829633 43 Γ 3361 Γ 3907 == 564651361 47 Γ 3359 Γ 6073 == 958762729 47 Γ 1151 Γ 1933 == 104569501 47 Γ 3727 Γ 5153 == 902645857 53 Γ 157 Γ 2081 == 17316001 53 Γ 79 Γ 599 == 2508013 53 Γ 157 Γ 521 == 4335241 59 Γ 1451 Γ 2089 == 178837201 61 Γ 421 Γ 12841 == 329769721 61 Γ 181 Γ 5521 == 60957361 61 Γ 1301 Γ 19841 == 1574601601 61 Γ 277 Γ 2113 == 35703361 61 Γ 181 Γ 1381 == 15247621 61 Γ 541 Γ 3001 == 99036001 61 Γ 661 Γ 2521 == 101649241 61 Γ 271 Γ 571 == 9439201 61 Γ 241 Γ 421 == 6189121 61 Γ 3361 Γ 4021 == 824389441 67 Γ 2311 Γ 51613 == 7991602081 67 Γ 331 Γ 7393 == 163954561 67 Γ 331 Γ 463 == 10267951 ## PicoLisp (de modulo (X Y) (% (+ Y (% X Y)) Y) ) (de prime? (N) (let D 0 (or (= N 2) (and (> N 1) (bit? 1 N) (for (D 3 T (+ D 2)) (T (> D (sqrt N)) T) (T (=0 (% N D)) NIL) ) ) ) ) ) (for P1 61 (when (prime? P1) (for (H3 2 (> P1 H3) (inc H3)) (let G (+ H3 P1) (for (D 1 (> G D) (inc D)) (when (and (=0 (% (* G (dec P1)) D) ) (= (modulo (* (- P1) P1) H3) (% D H3)) ) (let (P2 (inc (/ (* (dec P1) G) D) ) P3 (inc (/ (* P1 P2) H3)) ) (when (and (prime? P2) (prime? P3) (= 1 (modulo (* P2 P3) (dec P1))) ) (print (list P1 P2 P3)) ) ) ) ) ) ) ) )(prinl) (bye) ## PL/I Carmichael: procedure options (main, reorder); /* 24 January 2014 */ declare (Prime1, Prime2, Prime3, h3, d) fixed binary (31); put ('Carmichael numbers are:'); do Prime1 = 1 to 61; do h3 = 2 to Prime1; d_loop: do d = 1 to h3+Prime1-1; if (mod((h3+Prime1)*(Prime1-1), d) = 0) & (mod(-Prime1*Prime1, h3) = mod(d, h3)) then do; Prime2 = (Prime1-1) * (h3+Prime1)/d; Prime2 = Prime2 + 1; if ^is_prime(Prime2) then iterate d_loop; Prime3 = Prime1*Prime2/h3; Prime3 = Prime3 + 1; if ^is_prime(Prime3) then iterate d_loop; if mod(Prime2*Prime3, Prime1-1) ^= 1 then iterate d_loop; put skip edit (trim(Prime1), ' x ', trim(Prime2), ' x ', trim(Prime3)) (A); end; end; end; end; /* Uses is_prime from Rosetta Code PL/I. */ end Carmichael; Results: Carmichael numbers are: 3 x 11 x 17 5 x 29 x 73 5 x 17 x 29 5 x 13 x 17 7 x 19 x 67 7 x 31 x 73 7 x 13 x 31 7 x 23 x 41 7 x 73 x 103 7 x 13 x 19 9 x 89 x 401 9 x 29 x 53 13 x 61 x 397 13 x 37 x 241 13 x 97 x 421 13 x 37 x 97 13 x 37 x 61 17 x 41 x 233 17 x 353 x 1201 19 x 43 x 409 19 x 199 x 271 21 x 761 x 941 23 x 199 x 353 27 x 131 x 443 27 x 53 x 131 29 x 113 x 1093 29 x 197 x 953 31 x 991 x 15361 31 x 61 x 631 31 x 151 x 1171 31 x 61 x 271 31 x 61 x 211 31 x 271 x 601 31 x 181 x 331 35 x 647 x 7549 35 x 443 x 3877 37 x 109 x 2017 37 x 73 x 541 37 x 613 x 1621 37 x 73 x 181 37 x 73 x 109 41 x 1721 x 35281 41 x 881 x 12041 41 x 101 x 461 41 x 241 x 761 41 x 241 x 521 41 x 73 x 137 41 x 61 x 101 43 x 631 x 13567 43 x 271 x 5827 43 x 127 x 2731 43 x 127 x 1093 43 x 211 x 757 43 x 631 x 1597 43 x 127 x 211 43 x 211 x 337 43 x 433 x 643 43 x 547 x 673 43 x 3361 x 3907 47 x 3359 x 6073 47 x 1151 x 1933 47 x 3727 x 5153 49 x 313 x 5113 49 x 97 x 433 51 x 701 x 7151 53 x 157 x 2081 53 x 79 x 599 53 x 157 x 521 55 x 3079 x 84673 55 x 163 x 4483 55 x 1567 x 28729 55 x 109 x 1999 55 x 433 x 2647 55 x 919 x 3889 55 x 139 x 547 55 x 3889 x 12583 55 x 109 x 163 55 x 433 x 487 57 x 113 x 1289 57 x 113 x 281 57 x 4649 x 10193 59 x 1451 x 2089 61 x 421 x 12841 61 x 181 x 5521 61 x 1301 x 19841 61 x 277 x 2113 61 x 181 x 1381 61 x 541 x 3001 61 x 661 x 2521 61 x 271 x 571 61 x 241 x 421 61 x 3361 x 4021  ## Python class Isprime(): ''' Extensible sieve of Eratosthenes >>> isprime.check(11) True >>> isprime.multiples {2, 4, 6, 8, 9, 10} >>> isprime.primes [2, 3, 5, 7, 11] >>> isprime(13) True >>> isprime.multiples {2, 4, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 18, 20, 21, 22} >>> isprime.primes [2, 3, 5, 7, 11, 13, 17, 19] >>> isprime.nmax 22 >>> ''' multiples = {2} primes = [2] nmax = 2 def __init__(self, nmax): if nmax > self.nmax: self.check(nmax) def check(self, n): if type(n) == float: if not n.is_integer(): return False n = int(n) multiples = self.multiples if n <= self.nmax: return n not in multiples else: # Extend the sieve primes, nmax = self.primes, self.nmax newmax = max(nmax*2, n) for p in primes: multiples.update(range(p*((nmax + p + 1) // p), newmax+1, p)) for i in range(nmax+1, newmax+1): if i not in multiples: primes.append(i) multiples.update(range(i*2, newmax+1, i)) self.nmax = newmax return n not in multiples __call__ = check def carmichael(p1): ans = [] if isprime(p1): for h3 in range(2, p1): g = h3 + p1 for d in range(1, g): if (g * (p1 - 1)) % d == 0 and (-p1 * p1) % h3 == d % h3: p2 = 1 + ((p1 - 1)* g // d) if isprime(p2): p3 = 1 + (p1 * p2 // h3) if isprime(p3): if (p2 * p3) % (p1 - 1) == 1: #print('%i X %i X %i' % (p1, p2, p3)) ans += [tuple(sorted((p1, p2, p3)))] return ans isprime = Isprime(2) ans = sorted(sum((carmichael(n) for n in range(62) if isprime(n)), []))print(',\n'.join(repr(ans[i:i+5])[1:-1] for i in range(0, len(ans)+1, 5))) Output: (3, 11, 17), (5, 13, 17), (5, 17, 29), (5, 29, 73), (7, 13, 19), (7, 13, 31), (7, 19, 67), (7, 23, 41), (7, 31, 73), (7, 73, 103), (13, 37, 61), (13, 37, 97), (13, 37, 241), (13, 61, 397), (13, 97, 421), (17, 41, 233), (17, 353, 1201), (19, 43, 409), (19, 199, 271), (23, 199, 353), (29, 113, 1093), (29, 197, 953), (31, 61, 211), (31, 61, 271), (31, 61, 631), (31, 151, 1171), (31, 181, 331), (31, 271, 601), (31, 991, 15361), (37, 73, 109), (37, 73, 181), (37, 73, 541), (37, 109, 2017), (37, 613, 1621), (41, 61, 101), (41, 73, 137), (41, 101, 461), (41, 241, 521), (41, 241, 761), (41, 881, 12041), (41, 1721, 35281), (43, 127, 211), (43, 127, 1093), (43, 127, 2731), (43, 211, 337), (43, 211, 757), (43, 271, 5827), (43, 433, 643), (43, 547, 673), (43, 631, 1597), (43, 631, 13567), (43, 3361, 3907), (47, 1151, 1933), (47, 3359, 6073), (47, 3727, 5153), (53, 79, 599), (53, 157, 521), (53, 157, 2081), (59, 1451, 2089), (61, 181, 1381), (61, 181, 5521), (61, 241, 421), (61, 271, 571), (61, 277, 2113), (61, 421, 12841), (61, 541, 3001), (61, 661, 2521), (61, 1301, 19841), (61, 3361, 4021) ## Racket  #lang racket(require math) (for ([p1 (in-range 3 62)] #:when (prime? p1)) (for ([h3 (in-range 2 p1)]) (define g (+ p1 h3)) (let next ([d 1]) (when (< d g) (when (and (zero? (modulo (* g (- p1 1)) d)) (= (modulo (- (sqr p1)) h3) (modulo d h3))) (define p2 (+ 1 (quotient (* g (- p1 1)) d))) (when (prime? p2) (define p3 (+ 1 (quotient (* p1 p2) h3))) (when (and (prime? p3) (= 1 (modulo (* p2 p3) (- p1 1)))) (displayln (list p1 p2 p3 '=> (* p1 p2 p3)))))) (next (+ d 1))))))  Output:  (3 11 17 => 561)(5 29 73 => 10585)(5 17 29 => 2465)(5 13 17 => 1105)(7 19 67 => 8911)(7 31 73 => 15841)(7 23 41 => 6601)(7 73 103 => 52633)(13 61 397 => 314821)(13 97 421 => 530881)(13 37 97 => 46657)(13 37 61 => 29341)(17 41 233 => 162401)(17 353 1201 => 7207201)(19 43 409 => 334153)(19 199 271 => 1024651)(23 199 353 => 1615681)(29 113 1093 => 3581761)(29 197 953 => 5444489)(31 991 15361 => 471905281)(31 61 631 => 1193221)(31 151 1171 => 5481451)(31 61 271 => 512461)(31 61 211 => 399001)(31 271 601 => 5049001)(31 181 331 => 1857241)(37 109 2017 => 8134561)(37 73 541 => 1461241)(37 613 1621 => 36765901)(37 73 181 => 488881)(37 73 109 => 294409)(41 1721 35281 => 2489462641)(41 881 12041 => 434932961)(41 101 461 => 1909001)(41 241 761 => 7519441)(41 241 521 => 5148001)(41 73 137 => 410041)(41 61 101 => 252601)(43 631 13567 => 368113411)(43 127 1093 => 5968873)(43 211 757 => 6868261)(43 631 1597 => 43331401)(43 127 211 => 1152271)(43 211 337 => 3057601)(43 433 643 => 11972017)(43 547 673 => 15829633)(43 3361 3907 => 564651361)(47 3359 6073 => 958762729)(47 1151 1933 => 104569501)(47 3727 5153 => 902645857)(53 157 2081 => 17316001)(53 79 599 => 2508013)(53 157 521 => 4335241)(59 1451 2089 => 178837201)(61 421 12841 => 329769721)(61 1301 19841 => 1574601601)(61 277 2113 => 35703361)(61 541 3001 => 99036001)(61 661 2521 => 101649241)(61 271 571 => 9439201)(61 241 421 => 6189121)(61 3361 4021 => 824389441)  ## REXX Note that REXX's version of modulus (//) is really a remainder function, so a version of the modulus function was hard-coded below. (It was necessary to use modulus instead of remainder when using a negative value.) The Carmichael numbers are shown in numerical order. Some code optimization was done, while not necessary for the small default number (61), it was significant for larger numbers. /*REXX program calculates Carmichael 3-strong pseudoprimes (up to N).*/numeric digits 30 /*in case user wants bigger nums.*/parse arg N .; if N=='' then N=61 /*allow user to specify the limit*/if 7=='f7'x then times='af'x /*if EBCDIC machine, use a bullet*/ else times='f9'x /* " ASCII " " " " */carms=0 /*number of Carmichael #s so far.*/!.=0; !.2=1; !.3=1; !.5=1; !.7=1; !.11=1; !.13=1; !.17=1; !.19=1; !.23=1 /*[β] prime # memoization array.*/ do p=3 to N by 2; if \isPrime(p) then iterate /*Not prime? Skip*/ pm=p-1; nps=-p*p; bot=0; top=0 /*some handy-dandy REXX variables*/ @.=0 /*[β] Carmichael numbers are odd*/ do h3=2 to pm; g=h3+p /*find Carmichael #s for this P. */ gPM=g*pm; npsH3=((nps//h3)+h3)//h3 /*shortcuts.*/ do d=1 for g-1 if gPM//d \== 0 then iterate if npsH3 \== d//h3 then iterate q=1+gPM%d; if \isPrime(q) then iterate r=1+p*q%h3; if q*r//pm\==1 then iterate if \isPrime(r) then iterate carms=carms+1; @.q=r /*bump Carmichael #; add to array*/ if bot==0 then bot=q; bot=min(bot,q); top=max(top,q) end /*d*/ /* [β] find minimum and maximum.*/ end /*h3*/$=0                                /*display a list of some Carm #s.*/         do j=bot  to top  by 2;       if @.j==0  then iterate;        $=1 say 'ββββββββ a Carmichael number: ' p times j times @.j end /*j*/ if$  then say                     /*show beautification blank line.*/    end        /*p*/ say;     say carms ' Carmichael numbers found.'exit                                   /*stick a fork in it, we're done.*//*ββββββββββββββββββββββββββββββββββISPRIME subroutineββββββββββββββββββ*/isPrime: procedure expose !.;     parse arg x;   if !.x      then return 1if x<23 then return 0; if x//2==0 then return 0; if x// 3==0 then return 0if right(x,1)==5 then return 0;                  if x// 7==0 then return 0if x//11==0      then return 0;                  if x//13==0 then return 0if x//17==0      then return 0;                  if x//19==0 then return 0                   do i=23 by 6  until i*i>x; if x// i   ==0 then return 0                                              if x//(i+2)==0 then return 0                   end  /*i*/!.x=1;  return 1

output   when using the default input:

ββββββββ a Carmichael number:  3 β 11 β 17

ββββββββ a Carmichael number:  5 β 13 β 17
ββββββββ a Carmichael number:  5 β 17 β 29
ββββββββ a Carmichael number:  5 β 29 β 73

ββββββββ a Carmichael number:  7 β 13 β 19
ββββββββ a Carmichael number:  7 β 19 β 67
ββββββββ a Carmichael number:  7 β 23 β 41
ββββββββ a Carmichael number:  7 β 31 β 73
ββββββββ a Carmichael number:  7 β 73 β 103

ββββββββ a Carmichael number:  13 β 37 β 61
ββββββββ a Carmichael number:  13 β 61 β 397
ββββββββ a Carmichael number:  13 β 97 β 421

ββββββββ a Carmichael number:  17 β 41 β 233
ββββββββ a Carmichael number:  17 β 353 β 1201

ββββββββ a Carmichael number:  19 β 43 β 409
ββββββββ a Carmichael number:  19 β 199 β 271

ββββββββ a Carmichael number:  23 β 199 β 353

ββββββββ a Carmichael number:  29 β 113 β 1093
ββββββββ a Carmichael number:  29 β 197 β 953

ββββββββ a Carmichael number:  31 β 61 β 211
ββββββββ a Carmichael number:  31 β 151 β 1171
ββββββββ a Carmichael number:  31 β 181 β 331
ββββββββ a Carmichael number:  31 β 271 β 601
ββββββββ a Carmichael number:  31 β 991 β 15361

ββββββββ a Carmichael number:  37 β 73 β 109
ββββββββ a Carmichael number:  37 β 109 β 2017
ββββββββ a Carmichael number:  37 β 613 β 1621

ββββββββ a Carmichael number:  41 β 61 β 101
ββββββββ a Carmichael number:  41 β 73 β 137
ββββββββ a Carmichael number:  41 β 101 β 461
ββββββββ a Carmichael number:  41 β 241 β 521
ββββββββ a Carmichael number:  41 β 881 β 12041
ββββββββ a Carmichael number:  41 β 1721 β 35281

ββββββββ a Carmichael number:  43 β 127 β 211
ββββββββ a Carmichael number:  43 β 211 β 337
ββββββββ a Carmichael number:  43 β 271 β 5827
ββββββββ a Carmichael number:  43 β 433 β 643
ββββββββ a Carmichael number:  43 β 547 β 673
ββββββββ a Carmichael number:  43 β 631 β 1597
ββββββββ a Carmichael number:  43 β 3361 β 3907

ββββββββ a Carmichael number:  47 β 1151 β 1933
ββββββββ a Carmichael number:  47 β 3359 β 6073
ββββββββ a Carmichael number:  47 β 3727 β 5153

ββββββββ a Carmichael number:  53 β 79 β 599
ββββββββ a Carmichael number:  53 β 157 β 521

ββββββββ a Carmichael number:  59 β 1451 β 2089

ββββββββ a Carmichael number:  61 β 181 β 1381
ββββββββ a Carmichael number:  61 β 241 β 421
ββββββββ a Carmichael number:  61 β 271 β 571
ββββββββ a Carmichael number:  61 β 277 β 2113
ββββββββ a Carmichael number:  61 β 421 β 12841
ββββββββ a Carmichael number:  61 β 541 β 3001
ββββββββ a Carmichael number:  61 β 661 β 2521
ββββββββ a Carmichael number:  61 β 1301 β 19841
ββββββββ a Carmichael number:  61 β 3361 β 4021

69  Carmichael numbers found.


## Ruby

Works with: Ruby version 1.9
# Generate Charmichael Numbers require 'prime' Prime.each(61) do |p|  (2...p).each do |h3|    g = h3 + p    (1...g).each do |d|      next if (g*(p-1)) % d != 0 or (-p*p) % h3 != d % h3      q = 1 + ((p - 1) * g / d)      next unless q.prime?      r = 1 + (p * q / h3)      next unless r.prime? and (q * r) % (p - 1) == 1      puts "#{p} x #{q} x #{r}"     end  end  putsend
Output:
3 x 11 x 17

5 x 29 x 73
5 x 17 x 29
5 x 13 x 17

7 x 19 x 67
7 x 31 x 73
7 x 13 x 31
7 x 23 x 41
7 x 73 x 103
7 x 13 x 19

13 x 61 x 397
13 x 37 x 241
13 x 97 x 421
13 x 37 x 97
13 x 37 x 61

17 x 41 x 233
17 x 353 x 1201

19 x 43 x 409
19 x 199 x 271

23 x 199 x 353

29 x 113 x 1093
29 x 197 x 953

31 x 991 x 15361
31 x 61 x 631
31 x 151 x 1171
31 x 61 x 271
31 x 61 x 211
31 x 271 x 601
31 x 181 x 331

37 x 109 x 2017
37 x 73 x 541
37 x 613 x 1621
37 x 73 x 181
37 x 73 x 109

41 x 1721 x 35281
41 x 881 x 12041
41 x 101 x 461
41 x 241 x 761
41 x 241 x 521
41 x 73 x 137
41 x 61 x 101

43 x 631 x 13567
43 x 271 x 5827
43 x 127 x 2731
43 x 127 x 1093
43 x 211 x 757
43 x 631 x 1597
43 x 127 x 211
43 x 211 x 337
43 x 433 x 643
43 x 547 x 673
43 x 3361 x 3907

47 x 3359 x 6073
47 x 1151 x 1933
47 x 3727 x 5153

53 x 157 x 2081
53 x 79 x 599
53 x 157 x 521

59 x 1451 x 2089

61 x 421 x 12841
61 x 181 x 5521
61 x 1301 x 19841
61 x 277 x 2113
61 x 181 x 1381
61 x 541 x 3001
61 x 661 x 2521
61 x 271 x 571
61 x 241 x 421
61 x 3361 x 4021


## Seed7

The function isPrime below is borrowed from the Seed7 algorithm collection.

$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; Output: 3 * 11 * 17 = 561 5 * 29 * 73 = 10585 5 * 17 * 29 = 2465 5 * 13 * 17 = 1105 7 * 19 * 67 = 8911 7 * 31 * 73 = 15841 7 * 13 * 31 = 2821 7 * 23 * 41 = 6601 7 * 73 * 103 = 52633 7 * 13 * 19 = 1729 13 * 61 * 397 = 314821 13 * 37 * 241 = 115921 13 * 97 * 421 = 530881 13 * 37 * 97 = 46657 13 * 37 * 61 = 29341 17 * 41 * 233 = 162401 17 * 353 * 1201 = 7207201 19 * 43 * 409 = 334153 19 * 199 * 271 = 1024651 23 * 199 * 353 = 1615681 29 * 113 * 1093 = 3581761 29 * 197 * 953 = 5444489 31 * 991 * 15361 = 471905281 31 * 61 * 631 = 1193221 31 * 151 * 1171 = 5481451 31 * 61 * 271 = 512461 31 * 61 * 211 = 399001 31 * 271 * 601 = 5049001 31 * 181 * 331 = 1857241 37 * 109 * 2017 = 8134561 37 * 73 * 541 = 1461241 37 * 613 * 1621 = 36765901 37 * 73 * 181 = 488881 37 * 73 * 109 = 294409 41 * 1721 * 35281 = 2489462641 41 * 881 * 12041 = 434932961 41 * 101 * 461 = 1909001 41 * 241 * 761 = 7519441 41 * 241 * 521 = 5148001 41 * 73 * 137 = 410041 41 * 61 * 101 = 252601 43 * 631 * 13567 = 368113411 43 * 271 * 5827 = 67902031 43 * 127 * 2731 = 14913991 43 * 127 * 1093 = 5968873 43 * 211 * 757 = 6868261 43 * 631 * 1597 = 43331401 43 * 127 * 211 = 1152271 43 * 211 * 337 = 3057601 43 * 433 * 643 = 11972017 43 * 547 * 673 = 15829633 43 * 3361 * 3907 = 564651361 47 * 3359 * 6073 = 958762729 47 * 1151 * 1933 = 104569501 47 * 3727 * 5153 = 902645857 53 * 157 * 2081 = 17316001 53 * 79 * 599 = 2508013 53 * 157 * 521 = 4335241 59 * 1451 * 2089 = 178837201 61 * 421 * 12841 = 329769721 61 * 181 * 5521 = 60957361 61 * 1301 * 19841 = 1574601601 61 * 277 * 2113 = 35703361 61 * 181 * 1381 = 15247621 61 * 541 * 3001 = 99036001 61 * 661 * 2521 = 101649241 61 * 271 * 571 = 9439201 61 * 241 * 421 = 6189121 61 * 3361 * 4021 = 824389441  ## Sidef Translation of: Perl var ntheory = frequire('ntheory'); ntheory.forprimes({ |*p| p = Number.new(p[-1]); range(2, p-1).each { |h3| var ph3 = (p + h3); range(1, ph3-1).each { |d| ((-p * p) % h3) != (d % h3) && next; ((p-1)*ph3) % d && next; var q = 1+((p-1) * ph3 / d); ntheory.is_prime(q) || next; var r = 1+((p*q - 1)/h3); ntheory.is_prime(r) || next; (q*r) % (p-1) == 1 || next; printf("%2d x %5d x %5d = %s\n",p,q,r,ntheory.vecprod(p,q,r)); } }}, 3, 61); Output:  3 x 11 x 17 = 561 5 x 29 x 73 = 10585 5 x 17 x 29 = 2465 5 x 13 x 17 = 1105 ... full output is 69 lines ... 61 x 661 x 2521 = 101649241 61 x 271 x 571 = 9439201 61 x 241 x 421 = 6189121 61 x 3361 x 4021 = 824389441  ## Tcl Using the primality tester from the Miller-Rabin task... proc carmichael {limit {rounds 10}} { set carmichaels {} for {set p1 2} {$p1 <= $limit} {incr p1} { if {![miller_rabin$p1 $rounds]} continue for {set h3 2} {$h3 < $p1} {incr h3} { set g [expr {$h3 + $p1}] for {set d 1} {$d < $h3+$p1} {incr d} {		if {(($h3+$p1)*($p1-1))%$d != 0} continue		if {(-($p1**2))%$h3 != $d%$h3} continue 		set p2 [expr {1 + ($p1-1)*$g/$d}] if {![miller_rabin$p2 $rounds]} continue set p3 [expr {1 +$p1*$p2/$h3}]		if {![miller_rabin $p3$rounds]} continue 		if {($p2*$p3)%($p1-1) != 1} continue lappend carmichaels$p1 $p2$p3 [expr {$p1*$p2*$p3}] } } } return$carmichaels}

Demonstrating:

set results [carmichael 61 2]puts "[expr {[llength $results]/4}] Carmichael numbers found"foreach {p1 p2 p3 c}$results {    puts "$p1 x$p2 x $p3 =$c"}
Output:
69 Carmichael numbers found
3 x 11 x 17 = 561
5 x 29 x 73 = 10585
5 x 17 x 29 = 2465
5 x 13 x 17 = 1105
7 x 19 x 67 = 8911
7 x 31 x 73 = 15841
7 x 13 x 31 = 2821
7 x 23 x 41 = 6601
7 x 73 x 103 = 52633
7 x 13 x 19 = 1729
13 x 61 x 397 = 314821
13 x 37 x 241 = 115921
13 x 97 x 421 = 530881
13 x 37 x 97 = 46657
13 x 37 x 61 = 29341
17 x 41 x 233 = 162401
17 x 353 x 1201 = 7207201
19 x 43 x 409 = 334153
19 x 199 x 271 = 1024651
23 x 199 x 353 = 1615681
29 x 113 x 1093 = 3581761
29 x 197 x 953 = 5444489
31 x 991 x 15361 = 471905281
31 x 61 x 631 = 1193221
31 x 151 x 1171 = 5481451
31 x 61 x 271 = 512461
31 x 61 x 211 = 399001
31 x 271 x 601 = 5049001
31 x 181 x 331 = 1857241
37 x 109 x 2017 = 8134561
37 x 73 x 541 = 1461241
37 x 613 x 1621 = 36765901
37 x 73 x 181 = 488881
37 x 73 x 109 = 294409
41 x 1721 x 35281 = 2489462641
41 x 881 x 12041 = 434932961
41 x 101 x 461 = 1909001
41 x 241 x 761 = 7519441
41 x 241 x 521 = 5148001
41 x 73 x 137 = 410041
41 x 61 x 101 = 252601
43 x 631 x 13567 = 368113411
43 x 271 x 5827 = 67902031
43 x 127 x 2731 = 14913991
43 x 127 x 1093 = 5968873
43 x 211 x 757 = 6868261
43 x 631 x 1597 = 43331401
43 x 127 x 211 = 1152271
43 x 211 x 337 = 3057601
43 x 433 x 643 = 11972017
43 x 547 x 673 = 15829633
43 x 3361 x 3907 = 564651361
47 x 3359 x 6073 = 958762729
47 x 1151 x 1933 = 104569501
47 x 3727 x 5153 = 902645857
53 x 157 x 2081 = 17316001
53 x 79 x 599 = 2508013
53 x 157 x 521 = 4335241
59 x 1451 x 2089 = 178837201
61 x 421 x 12841 = 329769721
61 x 181 x 5521 = 60957361
61 x 1301 x 19841 = 1574601601
61 x 277 x 2113 = 35703361
61 x 181 x 1381 = 15247621
61 x 541 x 3001 = 99036001
61 x 661 x 2521 = 101649241
61 x 271 x 571 = 9439201
61 x 241 x 421 = 6189121
61 x 3361 x 4021 = 824389441


## zkl

Using the Miller-Rabin primality test in lib GMP.

var BN=Import("zklBigNum"), bi=BN(0); // gonna recycle biprimes:=T(2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61);var p2,p3;cs:=[[(p1,h3,d); primes; { [2..p1 - 1] }; // list comprehension      { [1..h3 + p1 - 1] },	{ ((h3 + p1)*(p1 - 1)%d == 0 and ((-p1*p1):mod(_,h3) == d%h3)) },//guard	{ (p2=1 + (p1 - 1)*(h3 + p1)/d):bi.set(_).probablyPrime() },//guard	{ (p3=1 + (p1*p2/h3)):bi.set(_).probablyPrime() },	 //guard	{ 1==(p2*p3)%(p1 - 1) };				 //guard   { T(p1,p2,p3) }  // return list of three primes in Carmichael number]];fcn mod(a,b) { m:=a%b; if(m<0) m+b else m }
cs.len().println(" Carmichael numbers found:");cs.pump(Console.println,fcn([(p1,p2,p3)]){   "%2d * %4d * %5d = %d".fmt(p1,p2,p3,p1*p2*p3) });
Output:
69 Carmichael numbers found:
3 *   11 *    17 = 561
5 *   29 *    73 = 10585
5 *   17 *    29 = 2465
5 *   13 *    17 = 1105
7 *   19 *    67 = 8911
...
61 *  181 *  1381 = 15247621
61 *  541 *  3001 = 99036001
61 *  661 *  2521 = 101649241
61 *  271 *   571 = 9439201
61 *  241 *   421 = 6189121
61 * 3361 *  4021 = 824389441