Carmichael 3 strong pseudoprimes

From Rosetta Code
Jump to: navigation, search
Task
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

[edit] Ada

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

[edit] 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

[edit] Haskell

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.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
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)

[edit] 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)%d
end
 
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
->

[edit] 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]]

[edit] 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,

[edit] 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

[edit] 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

[edit] 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)

[edit] 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)
 

[edit] 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.)

[edit] numbers in order of calculation

/*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

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.

[edit] numbers in sorted order

With a few lines of code (and using sparse arrays), the Carmichael numbers can be shown in order.

/*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

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.

[edit] Ruby

Works with: Ruby version 1.9
# 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 ""
}
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

[edit] 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

[edit] 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

[edit] zkl

Using the Miller-Rabin primality test in lib GMP.

var BN=Import("zklBigNum"), bi=BN(0); // gonna recycle bi
primes:=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
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox