Carmichael 3 strong pseudoprimes

A lot of composite numbers can be detected by the Miller Rabin Test, but there are some that evade it. 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 $Prime_{1}<Prime_{2}<Prime_{3}$ ) for all $Prime_{1}$ up to 61 (see page 7 of Notes by G.J.O Jameson March 2010 for solutions).

Pseudocode:
For a given $Prime_{1}$

for 1 < h3 < Prime₁

g=h3+Prime₁

for 0 < d < h3+Prime₁

if (h3+Prime₁)*(Prime₁-1) mod d == 0 and -Prime₁ squared mod h3 == d mod h3

then

Prime₂ = 1 + ((Prime₁-1) * g/d)

next d if Prime₂ is not prime

Prime₃ = 1 + (Prime₁*Prime₂/h3)

next d if Prime₃ is not prime

next d if (Prime₂*Prime₃) mod (Prime₁-1) not equal 1

Prime₁ * Prime₂ * Prime₃ is a Carmichael Number

D

From the Python entry, with several changes. Imports the third (extensible) D entry of the Sieve of Eratosthenes Task. <lang d>import std.stdio, std.typecons, sieve_of_eratosthenes3;

struct Carmichael {

 immutable int p1;

 static int mod(in int n, in int m) pure nothrow {
   return ((n % m) + m) % m;
 }

 int opApply(immutable int delegate(in Tuple!(int,int,int)) dg) {
   int result;

   if (IsPrime(p1))
     foreach (immutable h3; 2 .. p1) {
       immutable g = h3 + p1;
       foreach (immutable d; 1 .. g)
         if ((g * (p1 - 1)) % d == 0 && mod(-p1 * p1, h3) == d % h3) {
           immutable p2 = 1 + ((p1 - 1) * g / d);
           if (IsPrime(p2)) {
             immutable p3 = 1 + (p1 * p2 / h3);
             if (IsPrime(p3) && (p2 * p3) % (p1 - 1) == 1) {
               // A int[3] literal heap-allocates.
               const tri = Tuple!(int,int,int)(p1, p2, p3);
               result = dg(tri);
               if (result) break;
             }
           }
         }
     }

   return result;
 }

}

void main() {

 foreach (immutable n; 0 .. 62)
   if (IsPrime(n))
     foreach (const c; Carmichael(n))
       writefln("%(%d x %)", [c.tupleof]);

}</lang>

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

Haskell

Translation of: Ruby

Library: primes

Works with: GHC version 7.4.1

Works with: primes version 0.2.1.0

<lang haskell>#!/usr/bin/runhaskell

import Data.Numbers.Primes import 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</lang>

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)

Mathematica

<lang mathematica>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]]</lang>

Python

<lang python>class Isprime():

   
   Extensible sieve of Erastosthenes
   
   >>> 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)))</lang>

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)

REXX

Note that REXX's version of modulus (//) is really a remainder function, so a version of
the modulus function was hard-coded below (when using a negative value).

numbers in order of calculation

<lang rexx>/*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 1=='f1'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 /*a method of prime memoization. */

   do p=3  to N  by 2;  if \isPrime(p) then iterate  /*Not prime? Skip.*/
   pm=p-1; nps=-p*p; @.=0; min=1e9; max=0 /*some handy-dandy variables.*/
            do h3=2  to  pm;  g=h3+p  /*find Carmichael #s for this P. */
              do d=1  to g-1
              if g*pm//d\==0                 then iterate
              if ((nps//h3)+h3)//h3\==d//h3  then iterate
              q=1+pm*g%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           /*bump the Carmichael # counter. */
              min=min(min,q);  max=max(max,q);  @.q=r   /*build a list.*/
              end   /*d*/
            end     /*h3*/
                                      /*display a list of some Carm #s.*/
        do j=min to max by 2; if @.j==0 then iterate   /*one of the #s?*/
        say '──────── a Carmichael number: '   p   times  j   times   @.j
        end   /*j*/
   say                                /*show bueatification 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 1 if wordpos(x,'2 3 5 7 11 13')\==0 then do; !.x=1; return 1; end if x<17 then return 0; if x//2==0 then return 0; if x//3==0 then return 0 if right(x,1)==5 then return 0; if x//7==0 then return 0

                do i=11 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</lang> output when using the default input

──────── a Carmichael number:  3 ∙ 11 ∙ 17

──────── a Carmichael number:  5 ∙ 29 ∙ 73
──────── a Carmichael number:  5 ∙ 17 ∙ 29
──────── a Carmichael number:  5 ∙ 13 ∙ 17

──────── a Carmichael number:  7 ∙ 19 ∙ 67
──────── a Carmichael number:  7 ∙ 31 ∙ 73
──────── a Carmichael number:  7 ∙ 13 ∙ 31
──────── a Carmichael number:  7 ∙ 23 ∙ 41
──────── a Carmichael number:  7 ∙ 73 ∙ 103
──────── a Carmichael number:  7 ∙ 13 ∙ 19


──────── a Carmichael number:  13 ∙ 61 ∙ 397
──────── a Carmichael number:  13 ∙ 37 ∙ 241
──────── a Carmichael number:  13 ∙ 97 ∙ 421
──────── a Carmichael number:  13 ∙ 37 ∙ 97
──────── a Carmichael number:  13 ∙ 37 ∙ 61

──────── 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 ∙ 991 ∙ 15361
──────── a Carmichael number:  31 ∙ 61 ∙ 631
──────── a Carmichael number:  31 ∙ 151 ∙ 1171
──────── a Carmichael number:  31 ∙ 61 ∙ 271
──────── a Carmichael number:  31 ∙ 61 ∙ 211
──────── a Carmichael number:  31 ∙ 271 ∙ 601
──────── a Carmichael number:  31 ∙ 181 ∙ 331

──────── a Carmichael number:  37 ∙ 109 ∙ 2017
──────── a Carmichael number:  37 ∙ 73 ∙ 541
──────── a Carmichael number:  37 ∙ 613 ∙ 1621
──────── a Carmichael number:  37 ∙ 73 ∙ 181
──────── a Carmichael number:  37 ∙ 73 ∙ 109

──────── a Carmichael number:  41 ∙ 1721 ∙ 35281
──────── a Carmichael number:  41 ∙ 881 ∙ 12041
──────── a Carmichael number:  41 ∙ 101 ∙ 461
──────── a Carmichael number:  41 ∙ 241 ∙ 761
──────── a Carmichael number:  41 ∙ 241 ∙ 521
──────── a Carmichael number:  41 ∙ 73 ∙ 137
──────── a Carmichael number:  41 ∙ 61 ∙ 101

──────── a Carmichael number:  43 ∙ 631 ∙ 13567
──────── a Carmichael number:  43 ∙ 271 ∙ 5827
──────── a Carmichael number:  43 ∙ 127 ∙ 2731
──────── a Carmichael number:  43 ∙ 127 ∙ 1093
──────── a Carmichael number:  43 ∙ 211 ∙ 757
──────── a Carmichael number:  43 ∙ 631 ∙ 1597
──────── a Carmichael number:  43 ∙ 127 ∙ 211
──────── a Carmichael number:  43 ∙ 211 ∙ 337
──────── a Carmichael number:  43 ∙ 433 ∙ 643
──────── a Carmichael number:  43 ∙ 547 ∙ 673
──────── a Carmichael number:  43 ∙ 3361 ∙ 3907

──────── a Carmichael number:  47 ∙ 3359 ∙ 6073
──────── a Carmichael number:  47 ∙ 1151 ∙ 1933
──────── a Carmichael number:  47 ∙ 3727 ∙ 5153

──────── a Carmichael number:  53 ∙ 157 ∙ 2081
──────── a Carmichael number:  53 ∙ 79 ∙ 599
──────── a Carmichael number:  53 ∙ 157 ∙ 521

──────── a Carmichael number:  59 ∙ 1451 ∙ 2089

──────── a Carmichael number:  61 ∙ 421 ∙ 12841
──────── a Carmichael number:  61 ∙ 181 ∙ 5521
──────── a Carmichael number:  61 ∙ 1301 ∙ 19841
──────── a Carmichael number:  61 ∙ 277 ∙ 2113
──────── a Carmichael number:  61 ∙ 181 ∙ 1381
──────── a Carmichael number:  61 ∙ 541 ∙ 3001
──────── a Carmichael number:  61 ∙ 661 ∙ 2521
──────── a Carmichael number:  61 ∙ 271 ∙ 571
──────── a Carmichael number:  61 ∙ 241 ∙ 421
──────── a Carmichael number:  61 ∙ 3361 ∙ 4021


69  Carmichael numbers found.

numbers in sorted order

With a few lines of code (and using sparse arrays), the Carmichael numbers can be shown in order. <lang rexx>/*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 1=='f1'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 /*a method of prime memoization. */

                                      /*Carmichael numbers aren't even.*/
   do p=3  to N  by 2;  if \isPrime(p) then iterate  /*Not prime? Skip.*/
   pm=p-1; nps=-p*p; @.=0; min=1e9; max=0 /*some handy-dandy variables.*/

            do h3=2  to  pm;  g=h3+p  /*find Carmichael #s for this P. */
              do d=1  to g-1
              if g*pm//d\==0                 then iterate
              if ((nps//h3)+h3)//h3\==d//h3  then iterate
              q=1+pm*g%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           /*bump the Carmichael # counter. */
              min=min(min,q);  max=max(max,q);  @.q=r   /*build a list.*/
              end   /*d*/
            end     /*h3*/
                                      /*display a list of some Carm #s.*/
        do j=min to max by 2; if @.j==0 then iterate   /*one of the #s?*/
        say '──────── a Carmichael number: '   p   times  j   times   @.j
        end   /*j*/
   say                                /*show bueatification 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 1 if wordpos(x,'2 3 5 7 11 13')\==0 then do; !.x=1; return 1; end if x<17 then return 0; if x//2==0 then return 0; if x//3==0 then return 0 if right(x,1)==5 then return 0; if x//7==0 then return 0

                do i=11 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</lang> 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

<lang ruby># Generate Charmichael Numbers

Nigel_Galloway
November 30th., 2012.

require 'prime' Integer.each_prime(61) {|p|

 (2...p).each {|h3|
   g = h3+p
   (1...g).each {|d|
     if (g*(p-1)) % d == 0 and (-1*p*p) % h3 == d % h3 then
       q = 1+((p-1)*g/d)
       next if not q.prime?
       r = 1+(p*q/h3)
       next if not r.prime? or not (q*r) % (p-1) == 1
       puts "#{p} X #{q} X #{r}" 
     end
   }
 }
 puts ""

}</lang>

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