Carmichael 3 strong pseudoprimes
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.
You are encouraged to solve this task according to the task description, using any language you may know.
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 (where ) for all up to 61 (see page 7 of Notes by G.J.O Jameson March 2010 for solutions).
Pseudocode:
For a given
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
Ada
Uses the Miller_Rabin package from Miller-Rabin primality test#ordinary integers. <lang Ada>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;</lang>
- 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
D
This uses the third D entry of the Sieve of Eratosthenes Task. <lang d>import std.stdio, sieve_of_eratosthenes3;
int mod(in int n, in int m) pure nothrow {
return ((n % m) + m) % m;
}
void main() {
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); } } }
}</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
<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>
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.) <lang perl6>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}"; } } }
}</lang>
- 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
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.
(It was necessary to use modulus instead of remainder 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
The Carmichael numbers were grouped by the first Carmichael number.
──────── 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| next if (g*(p-1)) % d != 0 or (-1*p*p) % h3 != d % h3 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}" } } 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
Tcl
Using the primality tester from the Miller-Rabin task... <lang tcl>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
}</lang> Demonstrating: <lang tcl>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"
}</lang>
- 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