Carmichael 3 strong pseudoprimes: Difference between revisions
(→{{header|Fortran}}: Ah, division.) |
|||
Line 386: | Line 386: | ||
MOD(3 + MOD(D,3),3) 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 |
MOD(3 + MOD(D,3),3) 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 |
||
</pre> |
</pre> |
||
Thus, <code>MOD(-1, |
Thus, <code>MOD(-1,3) + MOD(3,3) = -1 + 0</code> but <code>MOD(3,3) = 0</code> so that ''(a + b) mod n'' does not come out as ''[(a mod n) + (b mod n)] mod n'' as might be expected. However, this is consistent with integer division that truncates. If ''q·f = n/d'' where q is the integer part and f the fractional part of the division, then the remainder is ''fd'' and has a sign too. |
||
Integer arithmetic can be a delicate matter, not just because of overflows. With division, sometimes the formula will never produce a remainder, as in n(n + 1)/2, and sometimes the divisor is always a factor for more complex reasons. But if a term produces a fractional part, is it to be discarded or not? Some languages (such as Pascal) insist on "div" instead of "/" so that one thinks about truncating division, others (such as pl/i) just generate the fractional part and proceed: perhaps the result will round? truncate? to the desired value, or perhaps other terms will cancel out the fractional parts and all will be well, or, perhaps they will combine and the result will be out by one, ''etc''. And some formulae require the fractional parts, as in the expression for the n'th Fibonacci number: |
Integer arithmetic can be a delicate matter, not just because of overflows. With division, sometimes the formula will never produce a remainder, as in n(n + 1)/2, and sometimes the divisor is always a factor for more complex reasons. But if a term produces a fractional part, is it to be discarded or not? Some languages (such as Pascal) insist on "div" instead of "/" so that one thinks about truncating division, others (such as pl/i) just generate the fractional part and proceed: perhaps the result will round? truncate? to the desired value, or perhaps other terms will cancel out the fractional parts and all will be well, or, perhaps they will combine and the result will be out by one, ''etc''. And some formulae require the fractional parts, as in the expression for the n'th Fibonacci number: |
Revision as of 02:38, 17 July 2016
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.
- Task
Find Carmichael numbers of the form:
- Prime1 × Prime2 × Prime3
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
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
ALGOL 68
Uses the Sieve of Eratosthenes code from the Smith Numbers task with an increased upper-bound (included here for convenience). <lang algol68># sieve of Eratosthene: sets s[i] to TRUE if i is prime, FALSE otherwise # PROC sieve = ( REF[]BOOL s )VOID:
BEGIN # start with everything flagged as prime # FOR i TO UPB s DO s[ i ] := TRUE OD; # sieve out the non-primes # s[ 1 ] := FALSE; FOR i FROM 2 TO ENTIER sqrt( UPB s ) DO IF s[ i ] THEN FOR p FROM i * i BY i TO UPB s DO s[ p ] := FALSE OD FI OD END # sieve # ;
- construct a sieve of primes up to the maximum number required for the task #
- For Prime1 up to 61 , we need to check numbers up to around 120 000 #
INT max number = 200 000; [ 1 : max number ]BOOL is prime; sieve( is prime );
- Find the Carmichael 3 Strong Pseudoprimes for Prime1 up to 61 #
FOR prime1 FROM 2 TO 61 DO
FOR h3 TO prime1 - 1 DO FOR d TO ( h3 + prime1 ) - 1 DO IF ( h3 + prime1 ) * ( prime1 - 1 ) MOD d = 0 AND ( - ( prime1 * prime1 ) ) MOD h3 = d MOD h3 THEN INT prime2 = 1 + ( ( prime1 - 1 ) * ( h3 + prime1 ) OVER d ); IF is prime[ prime2 ] THEN INT prime3 = 1 + ( prime1 * prime2 OVER h3 ); IF is prime[ prime3 ] THEN IF ( prime2 * prime3 ) MOD ( prime1 - 1 ) = 1 THEN print( ( whole( prime1, 0 ), " ", whole( prime2, 0 ), " ", whole( prime3, 0 ), newline ) ) FI FI FI FI OD OD
OD</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 9 89 401 9 29 53 ... 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
C
<lang 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;
} </lang>
- 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
<lang 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); } } }
}</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
EchoLisp
<lang scheme>
- 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)))
</lang>
- Output:
<lang scheme> (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 </lang>
Fortran
This is F77 style, and directly translates the given calculation. It turns out that the normal integers suffice for the demonstration, except for just one of the products of the three primes: 41x1721x35281 = 2489462641, which is bigger than 2147483647, the 32-bit limit. Fortunately, INTEGER*8 variables are also available, so the extension is easy. Otherwise, one would have to mess about with using two integers in a bignum style, one holding say the millions, and the second the number up to a million.
Directly translating the given formulae involves a trap. There is an expression -Prime1 squared mod h3, and translating this to MOD(-P1**2,H3)
may not work. The behaviour of the MOD function for negative numbers varies between language, compiler, and cpu. For positive numbers there is no confusion. However, in some cases one wants an affine shift as in day-of-the-week calculations from a daynumber, D: given MOD(D,7), one wants the same result with MOD(D - 7,7) or any other such shift and indeed, this happens with some versions of MOD. But not all, as Prof. Knuth remarks in his description of the calculation of the date for Easter. The MOD function can also work as follows:
D -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 MOD(D,3) -1 0 -2 -1 0 -2 -1 0 1 2 0 1 2 0 1 MOD(3 + MOD(D,3),3) 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1
Thus, MOD(-1,3) + MOD(3,3) = -1 + 0
but MOD(3,3) = 0
so that (a + b) mod n does not come out as [(a mod n) + (b mod n)] mod n as might be expected. However, this is consistent with integer division that truncates. If q·f = n/d where q is the integer part and f the fractional part of the division, then the remainder is fd and has a sign too.
Integer arithmetic can be a delicate matter, not just because of overflows. With division, sometimes the formula will never produce a remainder, as in n(n + 1)/2, and sometimes the divisor is always a factor for more complex reasons. But if a term produces a fractional part, is it to be discarded or not? Some languages (such as Pascal) insist on "div" instead of "/" so that one thinks about truncating division, others (such as pl/i) just generate the fractional part and proceed: perhaps the result will round? truncate? to the desired value, or perhaps other terms will cancel out the fractional parts and all will be well, or, perhaps they will combine and the result will be out by one, etc. And some formulae require the fractional parts, as in the expression for the n'th Fibonacci number: where each term generates an infinite non-repeating fractional part, and, they cancel exactly, even for large values of n...
Here I am supposing that (h3+Prime1)*(Prime1-1) mod d == 0 and -Prime1 squared mod h3 == d mod h3 is evaluated as {[(h3 + Prime1)*(Prime1 - 1) mod d] == 0} and {[-(Prime1 squared) mod h3] == d mod h3} which is to say it is [-(Prime1 squared)] mod h3 not -[(Prime1 squared) mod h3] where the first case involves the MOD(negative number) while the second is -MOD(positive number) and the only way that result could match d mod h3 is when both are zero.
Testing with the peculiar sign ignored via .AND. (MOD(P1**2,H3) .EQ. MOD(D,H3))) THEN
gives limited results:
Carmichael numbers that are the product of three primes: P1 x P2 x P3 = C 3 11 17 561 5 29 73 10585 7 19 67 8911 13 61 397 314821 13 37 241 115921 19 43 409 334153 31 991 15361 471905281 37 109 2017 8134561 41 1721 35281 2489462641 43 631 13567 368113411 43 271 5827 67902031 43 127 2731 14913991 61 421 12841 329769721 61 181 5521 60957361
Which turns out to be when the two terms have remainders of one. So, using the double MOD approach - which gives the same result for either style of MOD... <lang Fortran> LOGICAL FUNCTION ISPRIME(N) !Ad-hoc, since N is not going to be big...
INTEGER N !Despite this intimidating allowance of 32 bits... INTEGER F !A possible factor. ISPRIME = .FALSE. !Most numbers aren't prime. DO F = 2,SQRT(DFLOAT(N)) !Wince... IF (MOD(N,F).EQ.0) RETURN !Not even avoiding even numbers beyond two. END DO !Nice and brief, though. ISPRIME = .TRUE. !No factor found. END FUNCTION ISPRIME !So, done. Hopefully, not often.
PROGRAM CHASE INTEGER P1,P2,P3 !The three primes to be tested. INTEGER H3,D !Assistants. INTEGER MSG !File unit number. MSG = 6 !Standard output. WRITE (MSG,1) !A heading would be good. 1 FORMAT ("Carmichael numbers that are the product of three primes:" & /" P1 x P2 x P3 =",9X,"C") DO P1 = 2,61 !Step through the specified range. IF (ISPRIME(P1)) THEN !Selecting only the primes. DO H3 = 2,P1 - 1 !For 1 < H3 < P1. DO D = 1,H3 + P1 - 1 !For 0 < D < H3 + P1. IF (MOD((H3 + P1)*(P1 - 1),D).EQ.0 !Filter. & .AND. (MOD(H3 + MOD(-P1**2,H3),H3) .EQ. MOD(D,H3))) THEN !Beware MOD for negative numbers! MOD(-P1**2, may surprise... P2 = 1 + (P1 - 1)*(H3 + P1)/D !Candidate for the second prime. IF (ISPRIME(P2)) THEN !Is it prime? P3 = 1 + P1*P2/H3 !Yes. Candidate for the third prime. IF (ISPRIME(P3)) THEN !Is it prime? IF (MOD(P2*P3,P1 - 1).EQ.1) THEN !Yes! Final test. WRITE (MSG,2) P1,P2,P3, INT8(P1)*P2*P3 !Result! 2 FORMAT (3I6,I12) END IF END IF END IF END IF END DO END DO END IF END DO END</lang>
Output:
Carmichael numbers that are the product of three primes: P1 x P2 x P3 = C 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
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)
Icon and Unicon
The following works in both languages. <lang unicon>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</lang>
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 ->
Java
<lang java>public class Test {
static int mod(int n, int m) { return ((n % m) + m) % m; }
static boolean isPrime(int n) { if (n == 2 || n == 3) return true; else if (n < 2 || n % 2 == 0 || n % 3 == 0) return false; for (int div = 5, inc = 2; Math.pow(div, 2) <= n; div += inc, inc = 6 - inc) if (n % div == 0) return false; return true; }
public static void main(String[] args) { for (int p = 2; p < 62; p++) { if (!isPrime(p)) continue; for (int h3 = 2; h3 < p; h3++) { int g = h3 + p; for (int d = 1; d < g; d++) { if ((g * (p - 1)) % d != 0 || mod(-p * p, h3) != d % h3) continue; int q = 1 + (p - 1) * g / d; if (!isPrime(q)) continue; int r = 1 + (p * q / h3); if (!isPrime(r) || (q * r) % (p - 1) != 1) continue; System.out.printf("%d x %d x %d%n", p, q, r); } } } }
}</lang>
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
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 <lang Julia> 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 </lang>
Main <lang Julia> hi = 61 car = carmichael(hi)
curp = 0 tcnt = 0 print("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))
end println("\n\n", size(car)[2], " results in total.") </lang>
- 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
<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>
PARI/GP
<lang parigp>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]", ")))</lang>
- 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
<lang perl>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;</lang>
- 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.) <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
PicoLisp
<lang 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)</lang>
PL/I
<lang 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;</lang> 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
<lang 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)))</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)
Racket
<lang 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))))))
</lang> Output: <lang racket> (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) </lang>
REXX
slightly optimized
Note that REXX's version of modulus (//
) is really a remainder function.
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. <lang rexx>/*REXX program calculates Carmichael 3─strong pseudoprimes (up to and including N). */ numeric digits 18 /*handle big dig #s (9 is the default).*/ parse arg N .; if N== then N=61 /*allow user to specify for the search.*/ tell= N>0; N=abs(N) /*N>0? Then display Carmichael numbers*/
- =0 /*number of Carmichael numbers so far. */
@.=0; @.2=1; @.3=1; @.5=1; @.7=1; @.11=1; @.13=1; @.17=1; @.19=1; @.23=1; @.29=1; @.31=1
/*[↑] prime number memoization array. */ do p=3 to N by 2; pm=p-1; bot=0; top=0 /*step through some (odd) prime numbers*/ if \isPrime(p) then iterate; nps=-p*p /*is P a prime? No, then skip it.*/ !.=0 /*the list of Carmichael #'s (so far).*/ do h3=2 to pm; g=h3+p /*find Carmichael #s for this prime. */ gPM=g*pm; npsH3=((nps//h3)+h3)//h3 /*define a couple of shortcuts for pgm.*/ /* [↓] perform some weeding of D values*/ 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 #=#+1; !.q=r /*bump Carmichael counter; add to array*/ if bot==0 then bot=q; bot=min(bot,q); top=max(top,q) end /*d*/ end /*h3*/ $= /*display a list of some Carmichael #s.*/ do j=bot to top by 2 while tell; if !.j\==0 then $=$ p"∙"j'∙'!.j end /*j*/
if $\== then say 'Carmichael number: ' strip($) end /*p*/
say say '──────── ' # " Carmichael numbers found." exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ isPrime: parse arg x; if @.x then return 1 /*X a known prime?*/
if x<37 then return 0; if x//2==0 then return 0; if x// 3==0 then return 0 parse var x -1 _; if _==5 then return 0; if x// 7==0 then return 0 do k=11 by 6 until k*k>x; if x// k ==0 then return 0 if x//(k+2)==0 then return 0 end /*i*/ @.x=1; return 1</lang>
output when using the default input:
Carmichael number: 3∙11∙17 Carmichael number: 5∙13∙17 5∙17∙29 5∙29∙73 Carmichael number: 7∙13∙19 7∙19∙67 7∙23∙41 7∙31∙73 7∙73∙103 Carmichael number: 13∙37∙61 13∙61∙397 13∙97∙421 Carmichael number: 17∙41∙233 17∙353∙1201 Carmichael number: 19∙43∙409 19∙199∙271 Carmichael number: 23∙199∙353 Carmichael number: 29∙113∙1093 29∙197∙953 Carmichael number: 31∙61∙211 31∙151∙1171 31∙181∙331 31∙271∙601 31∙991∙15361 Carmichael number: 37∙73∙109 37∙109∙2017 37∙613∙1621 Carmichael number: 41∙61∙101 41∙73∙137 41∙101∙461 41∙241∙521 41∙881∙12041 41∙1721∙35281 Carmichael number: 43∙127∙211 43∙211∙337 43∙271∙5827 43∙433∙643 43∙547∙673 43∙631∙1597 43∙3361∙3907 Carmichael number: 47∙1151∙1933 47∙3359∙6073 47∙3727∙5153 Carmichael number: 53∙79∙599 53∙157∙521 Carmichael number: 59∙1451∙2089 Carmichael number: 61∙181∙1381 61∙241∙421 61∙271∙571 61∙277∙2113 61∙421∙12841 61∙541∙3001 61∙661∙2521 61∙1301∙19841 61∙3361∙4021 ──────── 69 Carmichael numbers found.
output when using the input of: -1000
──────── 1038 Carmichael numbers found.
output when using the input of: -10000
──────── 8716 Carmichael numbers found.
more optimized
This REXX version (pre-)generates a number of primes to assist the isPrime function. <lang rexx>/*REXX program calculates Carmichael 3─strong pseudoprimes (up to and including N). */ numeric digits 18 /*handle big dig #s (9 is the default).*/ parse arg N .; if N== then N=61 /*allow user to specify for the search.*/ tell= N>0; N=abs(N) /*N>0? Then display Carmichael numbers*/
- =0; @.=0 /*number of Carmichael numbers so far. */
@.2=1; @.3=1; @.5=1; @.7=1; @.11=1; @.13=1; @.17=1; @.19=1; @.23=1; @.29=1; @.31=1; @.37=1 HP=37; do i=HP+2 by 2 for N*20; if isPrime(i) then do; @.i=1; HP=i; end; end /*i*/ HP=HP+2
/*[↑] prime number memoization array. */ do p=3 to N by 2; pm=p-1; bot=0; top=0 /*step through some (odd) prime numbers*/ if \isPrime(p) then iterate; nps=-p*p /*is P a prime? No, then skip it.*/ !.=0 /*the list of Carmichael #'s (so far).*/ do h3=2 to pm; g=h3+p /*find Carmichael #s for this prime. */ gPM=g*pm; npsH3=((nps//h3)+h3)//h3 /*define a couple of shortcuts for pgm.*/ /* [↓] perform some weeding of D values*/ 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 #=#+1; !.q=r /*bump Carmichael counter; add to array*/ if bot==0 then bot=q; bot=min(bot,q); top=max(top,q) end /*d*/ end /*h3*/ $= /*display a list of some Carmichael #s.*/ do j=bot to top by 2 while tell; if !.j\==0 then $=$ p"∙"j'∙'!.j end /*j*/
if $\== then say 'Carmichael number: ' strip($) end /*p*/
say say '──────── ' # " Carmichael numbers found." exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ isPrime: parse arg x; if @.x then return 1 /*X a known prime?*/
if x<HP then return 0; if x//2==0 then return 0; if x// 3==0 then return 0 parse var x -1 _; if _==5 then return 0; if x// 7==0 then return 0 if x//11==0 then return 0; if x//13==0 then return 0 if x//17==0 then return 0; if x//19==0 then return 0 if x//23==0 then return 0; if x//29==0 then return 0 if x//31==0 then return 0; if x//37==0 then return 0
c=x /* ___*/ b=1; do while b<=x; b=b*4; end /*these two lines compute integer √ X */ s=0; do while b>1; b=b%4; _=c-s-b; s=s%2; if _>=0 then do; c=_; s=s+b; end; end
do k=41 by 6 to s; parse var k -1 _ if _\==5 then if x// k ==0 then return 0 if _\==3 then if x//(k+2)==0 then return 0 end /*k*/ /*K will never be divisible by three.*/ @.x=1; return 1 /*Define a new prime (X). Indicate so.*/</lang>
output when the following us used for input: 1000
Carmichael number: 3∙11∙17 Carmichael number: 5∙13∙17 5∙17∙29 5∙29∙73 Carmichael number: 7∙13∙19 7∙19∙67 7∙23∙41 7∙31∙73 7∙73∙103 Carmichael number: 13∙37∙61 13∙61∙397 13∙97∙421 Carmichael number: 17∙41∙233 17∙353∙1201 Carmichael number: 19∙43∙409 19∙199∙271 Carmichael number: 23∙199∙353 Carmichael number: 29∙113∙1093 29∙197∙953 Carmichael number: 31∙61∙211 31∙151∙1171 31∙181∙331 31∙271∙601 31∙991∙15361 Carmichael number: 37∙73∙109 37∙109∙2017 37∙613∙1621 Carmichael number: 41∙61∙101 41∙73∙137 41∙101∙461 41∙241∙521 41∙881∙12041 41∙1721∙35281 Carmichael number: 43∙127∙211 43∙211∙337 43∙271∙5827 43∙433∙643 43∙547∙673 43∙631∙1597 43∙3361∙3907 Carmichael number: 47∙1151∙1933 47∙3359∙6073 47∙3727∙5153 Carmichael number: 53∙79∙599 53∙157∙521 Carmichael number: 59∙1451∙2089 Carmichael number: 61∙181∙1381 61∙241∙421 61∙271∙571 61∙277∙2113 61∙421∙12841 61∙541∙3001 61∙661∙2521 61∙1301∙19841 61∙3361∙4021 Carmichael number: 67∙331∙463 67∙2311∙51613 Carmichael number: 71∙271∙521 71∙421∙491 71∙631∙701 71∙701∙5531 71∙911∙9241 Carmichael number: 73∙157∙2293 73∙379∙523 73∙601∙21937 73∙937∙13681 Carmichael number: 79∙2237∙25247 Carmichael number: 83∙2953∙4019 83∙6971∙289297 Carmichael number: 89∙353∙617 89∙617∙3433 89∙881∙7129 89∙4049∙120121 Carmichael number: 97∙193∙1249 97∙673∙769 97∙769∙10657 97∙1201∙38833 97∙2113∙5857 Carmichael number: 101∙151∙251 101∙1301∙43801 Carmichael number: 103∙647∙1361 103∙1123∙6427 103∙3571∙183907 103∙3877∙4591 103∙5407∙185641 Carmichael number: 107∙743∙1061 107∙3181∙26183 Carmichael number: 109∙163∙379 109∙229∙4993 109∙241∙2389 109∙379∙919 109∙433∙2053 109∙541∙2269 109∙1297∙12853 109∙1657∙6229 109∙1801∙4789 Carmichael number: 113∙337∙449 113∙827∙18691 113∙6833∙85793 113∙8737∙22961 Carmichael number: 127∙631∙26713 127∙659∙1373 127∙991∙3313 127∙2143∙30241 127∙16633∙422479 Carmichael number: 131∙491∙4021 131∙521∙2731 131∙571∙1871 131∙911∙2341 131∙1171∙38351 131∙1301∙8971 131∙4421∙115831 131∙5851∙191621 131∙17291∙1132561 Carmichael number: 137∙409∙14009 137∙953∙5441 137∙3673∙20129 137∙5441∙11833 Carmichael number: 139∙691∙829 139∙3359∙66701 139∙4003∙92737 139∙4969∙8971 139∙17389∙21391 Carmichael number: 149∙593∙29453 149∙1481∙3109 Carmichael number: 151∙211∙541 151∙331∙3571 151∙601∙751 151∙701∙1451 151∙751∙8101 151∙2551∙192601 151∙4951∙53401 Carmichael number: 157∙313∙1093 157∙937∙6397 157∙1093∙15601 157∙2017∙28789 157∙4993∙11701 157∙26053∙409033 Carmichael number: 163∙379∙2377 163∙487∙811 163∙811∙1297 163∙1297∙42283 163∙1621∙37747 Carmichael number: 167∙4483∙34031 167∙29383∙490697 Carmichael number: 173∙1291∙31907 173∙10321∙255077 173∙36809∙155317 Carmichael number: 179∙1069∙4451 179∙3739∙7121 179∙9257∙57139 179∙10859∙485941 Carmichael number: 181∙271∙9811 181∙541∙3061 181∙631∙811 181∙733∙66337 181∙1693∙43777 181∙2953∙3637 181∙3061∙9721 181∙5881∙9421 181∙11161∙15661 181∙16561∙999181 Carmichael number: 191∙421∙431 191∙571∙15581 191∙1901∙7411 191∙3877∙56963 191∙12541∙342191 Carmichael number: 193∙257∙1601 193∙401∙11057 193∙577∙5569 193∙1249∙2593 193∙2689∙30529 193∙7681∙211777 193∙61057∙94273 Carmichael number: 199∙397∙4159 199∙859∙2311 199∙937∙20719 199∙991∙32869 199∙4159∙8713 199∙8713∙82567 Carmichael number: 211∙281∙491 211∙421∙631 211∙631∙66571 211∙1051∙9241 211∙1741∙4651 211∙1831∙4111 211∙2311∙54181 211∙2521∙3571 211∙3221∙35771 211∙3571∙9661 211∙4019∙11159 211∙4201∙98491 211∙32341∙70351 211∙68041∙127051 Carmichael number: 223∙1777∙23311 223∙2221∙5107 223∙5107∙37963 223∙14653∙79699 223∙25087∙1864801 Carmichael number: 227∙1583∙6781 227∙2713∙5651 227∙7459∙423299 227∙21019∙91757 Carmichael number: 229∙457∙2053 229∙571∙4219 229∙1597∙182857 229∙2243∙73379 229∙7069∙32377 Carmichael number: 233∙5569∙185369 Carmichael number: 239∙409∙3911 239∙1429∙1667 239∙1667∙6427 239∙32369∙234431 Carmichael number: 241∙337∙4513 241∙433∙52177 241∙1201∙2161 241∙1361∙4001 241∙2161∙3361 241∙3361∙32401 241∙5521∙110881 241∙6481∙780961 241∙7321∙588121 241∙84961∙181201 Carmichael number: 251∙751∙3251 251∙2251∙56501 251∙3251∙4001 251∙4751∙22501 251∙16001∙803251 251∙22501∙297251 251∙31751∙2656501 Carmichael number: 257∙641∙1153 257∙769∙49409 257∙67073∙3447553 257∙78593∙403969 Carmichael number: 263∙787∙2621 263∙1049∙3407 263∙6551∙12577 263∙71527∙1881161 Carmichael number: 269∙1877∙126229 269∙4289∙384581 269∙10453∙65393 269∙15277∙21977 Carmichael number: 271∙541∙811 271∙811∙2971 271∙1171∙7741 271∙1801∙16831 271∙2161∙65071 271∙4591∙4861 271∙4861∙77491 271∙8191∙1109881 271∙8641∙47791 271∙9631∙52201 271∙10531∙1426951 271∙14797∙1336663 Carmichael number: 277∙829∙7177 277∙1381∙1933 277∙3313∙9661 277∙6073∙31741 277∙18493∙88321 277∙19597∙36433 Carmichael number: 281∙421∙701 281∙617∙673 281∙1009∙4649 281∙1321∙23201 281∙4201∙9521 281∙7121∙26681 281∙9521∙13721 281∙25621∙84701 Carmichael number: 283∙4231∙598687 283∙17203∙58657 Carmichael number: 293∙877∙4673 293∙1607∙31391 293∙3943∙5987 Carmichael number: 307∙613∙919 307∙919∙141067 307∙1531∙3673 307∙2143∙13159 307∙3673∙225523 307∙6427∙246637 307∙17443∙153001 307∙18973∙1941571 Carmichael number: 311∙1117∙26723 311∙1303∙2357 311∙2791∙21701 311∙3659∙7069 311∙23251∙33791 311∙26041∙323951 311∙28211∙165541 311∙44641∙52391 Carmichael number: 313∙521∙23297 313∙937∙58657 313∙1093∙6709 313∙1249∙55849 313∙3433∙38377 313∙3793∙395737 313∙5449∙12097 313∙6577∙8761 313∙7177∙70201 313∙9049∙472057 313∙12637∙359581 313∙49297∙5143321 313∙51481∙947857 313∙66457∙184081 313∙129169∙400297 Carmichael number: 317∙18013∙41081 317∙104281∙2542853 Carmichael number: 331∙661∙991 331∙947∙24113 331∙991∙4621 331∙1321∙17491 331∙2311∙12541 331∙2971∙49171 331∙3331∙551281 331∙4051∙18121 331∙4621∙14851 331∙37181∙1758131 331∙37951∙897271 331∙41141∙316691 Carmichael number: 337∙421∙47293 337∙449∙21617 337∙673∙1009 337∙953∙8681 337∙1009∙3697 337∙1597∙12517 337∙2473∙11113 337∙3697∙12097 337∙16369∙1379089 337∙19489∙597073 337∙35393∙40433 337∙58129∙2176609 Carmichael number: 347∙3461∙92383 347∙4153∙29411 Carmichael number: 349∙929∙7541 349∙1741∙2089 349∙3191∙12239 349∙4177∙20533 349∙122149∙21315001 Carmichael number: 353∙617∙19801 353∙1153∙5153 353∙1321∙66617 353∙13217∙77761 353∙130241∙2704417 Carmichael number: 359∙43319∙3887881 359∙46183∙592133 Carmichael number: 367∙733∙1831 367∙1831∙9883 367∙5003∙42701 367∙9151∙419803 367∙24889∙51607 367∙28183∙574621 Carmichael number: 373∙1117∙1861 373∙1613∙150413 373∙5581∙1040857 373∙16741∙81097 373∙139501∙26016937 Carmichael number: 379∙631∙9199 379∙757∙4159 379∙2269∙24571 379∙2539∙21871 379∙6427∙202987 379∙9829∙17011 379∙10639∙268813 Carmichael number: 383∙33617∙40111 383∙38201∙860647 383∙74873∙3186263 Carmichael number: 389∙3299∙6791 Carmichael number: 397∙1783∙4951 397∙2971∙51283 397∙4357∙8317 397∙30097∙56629 397∙55837∙852589 397∙79201∙10480933 397∙99793∙370261 Carmichael number: 401∙641∙2161 401∙1201∙1601 401∙2161∙216641 401∙2801∙9601 401∙9521∙19681 401∙9601∙70001 401∙15601∙18401 401∙18401∙567601 401∙161201∙32320801 Carmichael number: 409∙2857∙6529 409∙6121∙96289 409∙6529∙22441 409∙7039∙575791 409∙35089∙683401 409∙36721∙114649 Carmichael number: 419∙15467∙47653 419∙22573∙78167 419∙47653∙539639 Carmichael number: 421∙631∙11551 421∙701∙2381 421∙3851∙85331 421∙7561∙289381 421∙9661∙15121 421∙13441∙209581 421∙18481∙39901 421∙20231∙54251 421∙35533∙7479697 421∙42589∙208489 421∙89041∙12495421 Carmichael number: 431∙1721∙29671 431∙1979∙142159 431∙8171∙55901 431∙13331∙168991 Carmichael number: 433∙937∙11593 433∙1297∙2161 433∙2161∙16417 433∙2593∙48817 433∙2953∙21673 433∙3457∙6481 433∙3697∙55201 433∙6481∙87697 Carmichael number: 439∙3067∙673207 439∙3943∙45553 439∙9199∙2019181 439∙10513∙17959 439∙64679∙7098521 439∙96799∙14164921 Carmichael number: 443∙1327∙4421 443∙2029∙4967 443∙74257∙143651 443∙102103∙2380613 443∙167077∙236471 443∙251057∙889747 Carmichael number: 449∙2689∙3137 449∙50849∙4566241 449∙145601∙325249 449∙202049∙45360001 Carmichael number: 457∙3877∙93253 457∙5701∙8893 457∙7297∙32377 457∙15733∙19381 457∙21433∙163249 457∙28729∙55633 457∙71593∙2337001 457∙73721∙1203233 457∙114001∙1211593 Carmichael number: 461∙691∙1151 461∙1013∙38917 461∙1381∙159161 461∙3541∙23321 461∙5981∙24841 461∙26681∙4099981 Carmichael number: 463∙2927∙15401 463∙6007∙39733 463∙214831∙49733377 463∙218527∙10117801 Carmichael number: 467∙141199∙474389 Carmichael number: 479∙57839∙219881 Carmichael number: 487∙1459∙8263 487∙1531∙2683 487∙1621∙1783 487∙1783∙108541 487∙8263∙9721 487∙12637∙32563 487∙17011∙2761453 487∙26731∙110323 487∙51517∙69499 Carmichael number: 491∙1471∙10781 491∙6959∙569479 491∙16661∙154351 491∙41651∙46061 491∙122501∙6683111 491∙386611∙637001 Carmichael number: 499∙997∙4483 499∙10459∙39841 Carmichael number: 503∙5021∙21587 Carmichael number: 509∙3557∙41149 509∙7621∙23369 509∙11939∙110491 509∙86869∙11054081 Carmichael number: 521∙1301∙8581 521∙21841∙41081 Carmichael number: 523∙1567∙163909 523∙6091∙1592797 523∙9397∙140419 523∙15661∙481807 523∙38629∙69427 523∙155557∙1114471 523∙193663∙462493 Carmichael number: 541∙811∙3511 541∙1621∙7561 541∙6661∙257401 541∙7561∙54541 541∙12421∙197641 541∙16561∙814501 Carmichael number: 547∙1093∙2731 547∙2731∙6553 547∙6553∙35491 547∙7333∙235951 547∙26209∙186187 547∙52963∙827737 547∙158341∙2624623 Carmichael number: 557∙1669∙42257 557∙38921∙7226333 Carmichael number: 563∙28663∙329333 Carmichael number: 569∙2273∙117577 569∙13633∙1108169 569∙17609∙25561 569∙21017∙37489 569∙22153∙787817 Carmichael number: 571∙661∙16411 571∙2281∙2851 571∙2851∙13681 571∙6841∙43891 571∙13681∙1562371 571∙65323∙18649717 Carmichael number: 577∙757∙39709 577∙1153∙6337 577∙5569∙100417 577∙6337∙26497 577∙20161∙646273 577∙32833∙37441 577∙53857∙181729 577∙79777∙86689 577∙339841∙15083713 577∙559297∙819073 Carmichael number: 587∙5861∙33403 587∙9377∙54499 587∙12893∙36919 587∙49811∙3654883 Carmichael number: 593∙21017∙31081 593∙35521∙3009137 593∙176417∙34871761 Carmichael number: 599∙2393∙84319 599∙120199∙17999801 599∙179999∙35939801 599∙266111∙547769 599∙368369∙12979591 Carmichael number: 601∙1201∙1801 601∙1801∙541201 601∙3001∙200401 601∙3121∙38281 601∙3301∙5101 601∙4201∙4801 601∙4801∙412201 601∙5101∙278701 601∙6151∙7951 601∙9001∙386401 601∙19801∙28201 601∙52201∙3921601 601∙99901∙923701 Carmichael number: 607∙1213∙9091 607∙4243∙1287751 607∙21817∙322999 607∙24847∙1885267 607∙61813∙7504099 607∙186649∙12588439 607∙370873∙45023983 607∙373903∙22695913 Carmichael number: 613∙919∙2143 613∙1021∙312937 613∙1327∙73951 613∙1429∙23053 613∙2857∙17341 613∙7549∙87313 613∙9181∙2813977 613∙12241∙111997 613∙51817∙213181 613∙246637∙783361 613∙364753∙386173 Carmichael number: 617∙661∙1013 617∙8009∙705937 617∙16633∙120737 617∙29569∙2606297 617∙59753∙81929 617∙73613∙133981 617∙129361∙6139673 617∙137369∙1629937 617∙383153∙47281081 Carmichael number: 619∙1237∙4327 619∙2267∙26987 619∙5563∙1721749 619∙28429∙703903 619∙53149∙56239 619∙92083∙452377 619∙398611∙9490009 Carmichael number: 631∙1471∙46411 631∙5881∙90511 631∙26209∙82279 631∙32831∙67481 Carmichael number: 641∙4481∙7681 641∙12161∙26881 641∙17921∙370561 641∙19841∙176641 Carmichael number: 643∙107857∙2391451 Carmichael number: 647∙4523∙19381 647∙64601∙75583 647∙188633∙532951 647∙444449∙7013623 Carmichael number: 653∙13367∙2909551 653∙176041∙732197 Carmichael number: 659∙2633∙5923 659∙23689∙624443 659∙27919∙34781 659∙30269∙92779 659∙73039∙6876101 659∙92779∙1329161 Carmichael number: 661∙991∙131011 661∙1321∙4621 661∙2131∙4231 661∙3191∙6491 661∙3301∙12541 661∙4621∙763621 661∙5281∙81181 661∙22111∙1623931 661∙22441∙95701 661∙138821∙152681 Carmichael number: 673∙1009∙14449 673∙2017∙3361 673∙3361∙12097 673∙13441∙1292257 673∙40801∙155137 673∙231841∙9178177 Carmichael number: 677∙2029∙85853 677∙4733∙1602121 677∙6761∙25013 677∙45293∙511057 Carmichael number: 683∙8867∙16369 683∙11161∙206027 683∙15749∙32303 683∙42967∙2934647 683∙94117∙9183131 Carmichael number: 691∙7591∙2622691 691∙16561∙2288731 691∙31051∙71761 691∙34501∙2648911 691∙69691∙3009781 691∙743131∙1330321 Carmichael number: 701∙2801∙10501 701∙3701∙1297201 701∙3851∙899851 701∙6301∙7001 701∙18401∙58901 701∙41651∙2245951 701∙44101∙170801 701∙46901∙319201 701∙52501∙296801 701∙53201∙632101 Carmichael number: 709∙4957∙12037 709∙7789∙16993 709∙9677∙21713 709∙36109∙5120257 709∙210277∙819157 Carmichael number: 719∙97649∙190271 Carmichael number: 727∙1453∙2179 727∙2179∙792067 727∙2663∙193601 727∙3631∙8713 727∙4423∙321553 727∙176903∙32152121 727∙308551∙1823713 727∙651223∙2784937 Carmichael number: 733∙5857∙84181 733∙13177∙47581 733∙18301∙789097 733∙22571∙2363507 733∙25621∙9390097 733∙150427∙1238911 733∙271573∙22118113 733∙631717∙3561913 Carmichael number: 739∙821∙4019 739∙3691∙454609 739∙10333∙2545363 739∙62731∙1783009 739∙152029∙1321759 Carmichael number: 743∙6679∙225569 743∙6997∙9011 743∙596569∙7266407 Carmichael number: 751∙2251∙10501 751∙2851∙237901 751∙21751∙181501 751∙109751∙649001 751∙123001∙1338751 751∙153001∙1767751 751∙191251∙10259251 751∙318751∙2418001 Carmichael number: 757∙2017∙18397 757∙2269∙858817 757∙15121∙3815533 757∙27541∙79273 757∙32257∙2219869 757∙33013∙59221 757∙184843∙633151 757∙627481∙6506893 Carmichael number: 761∙2129∙31769 761∙2281∙3041 761∙3041∙771401 761∙6841∙19001 761∙8969∙1137569 761∙13681∙101081 761∙19001∙1032841 761∙41801∙497041 761∙230281∙1184081 761∙251941∙339341 761∙314641∙497801 Carmichael number: 769∙6529∙9601 769∙41729∙697601 Carmichael number: 773∙22003∙122363 773∙44777∙47093 Carmichael number: 787∙3931∙9433 787∙5503∙45589 787∙106373∙3348623 Carmichael number: 797∙2389∙476009 797∙3583∙16319 797∙5573∙11941 797∙21493∙428249 797∙58109∙7718813 797∙148853∙859681 Carmichael number: 809∙5657∙9697 809∙78781∙176549 809∙82013∙22116173 809∙176549∙2197357 809∙453289∙1171601 Carmichael number: 811∙1621∙438211 811∙4051∙19441 811∙4591∙744661 811∙6481∙17011 811∙19441∙3153331 811∙77761∙1189891 811∙86131∙478441 Carmichael number: 821∙1231∙6971 821∙15581∙42641 821∙137597∙6275953 Carmichael number: 823∙2467∙4111 823∙4111∙23017 823∙4933∙9043 823∙27127∙637873 823∙341953∙31269703 Carmichael number: 827∙2243∙2833 827∙4957∙5783 827∙24781∙476603 827∙101009∙2880499 827∙691363∙57175721 Carmichael number: 829∙1657∙17389 829∙9109∙15733 829∙10949∙2269181 829∙24841∙1872109 829∙140761∙5556709 Carmichael number: 839∙5867∙223747 Carmichael number: 853∙2557∙4261 853∙7669∙594697 853∙12781∙5451097 853∙17041∙309277 853∙19597∙185737 Carmichael number: 857∙6421∙127973 857∙10273∙160073 857∙95873∙115561 857∙796937∙9229393 Carmichael number: 859∙2861∙3719 859∙8581∙9439 859∙9439∙150151 859∙27457∙66067 859∙321751∙1039039 Carmichael number: 863∙24137∙38791 863∙28447∙153437 863∙38791∙62927 863∙56893∙68099 Carmichael number: 877∙1753∙56941 877∙3067∙30223 877∙6133∙8761 877∙24091∙7042603 877∙36793∙6453493 877∙263677∙8894029 Carmichael number: 881∙2861∙840181 881∙22441∙57641 881∙130241∙16391761 Carmichael number: 883∙2647∙44101 883∙8191∙267877 883∙11467∙35281 883∙15877∙824671 883∙16633∙358219 883∙21757∙3842287 883∙30871∙134947 883∙42337∙216091 883∙126127∙161407 883∙260191∙114874327 883∙403957∙10808911 883∙507151∙531847 Carmichael number: 887∙14177∙50503 Carmichael number: 907∙7853∙16007 907∙137713∙24981139 Carmichael number: 911∙2003∙912367 911∙9283∙1208117 911∙9311∙55441 911∙11831∙898171 911∙16381∙28211 911∙30941∙4026751 911∙55511∙12642631 911∙167441∙204751 911∙175631∙2962961 911∙185641∙1551551 911∙227501∙2328691 Carmichael number: 919∙8263∙949213 919∙15607∙170749 919∙60589∙11136259 919∙129439∙569161 919∙156979∙321301 919∙311203∙2918323 919∙877609∙21797911 Carmichael number: 929∙5569∙23201 929∙6961∙35729 929∙42689∙1071841 929∙139201∙307169 Carmichael number: 937∙1873∙70201 937∙6553∙7489 937∙7489∙1002457 937∙21529∙3362113 937∙38377∙5993209 937∙177841∙820873 Carmichael number: 941∙5171∙23971 941∙6581∙8461 941∙8461∙361901 941∙28201∙102461 941∙44651∙4668511 941∙209621∙1133641 941∙322891∙701711 941∙355321∙1732421 Carmichael number: 947∙29327∙1983763 947∙47129∙299539 947∙307451∙10398433 Carmichael number: 953∙2857∙9521 953∙5881∙18257 953∙17137∙69497 953∙52361∙159937 953∙159937∙2771273 Carmichael number: 967∙1289∙25439 967∙1933∙4831 967∙4831∙11593 967∙26083∙5044453 967∙62791∙7589863 967∙88873∙1909783 967∙156493∙30265747 Carmichael number: 971∙3881∙753691 971∙8731∙44621 971∙12611∙3061321 971∙110581∙635351 971∙142591∙2387171 971∙169751∙648931 971∙1324051∙3263081 Carmichael number: 977∙2441∙794953 977∙5857∙12689 977∙6833∙39041 977∙17569∙41969 977∙478241∙155747153 Carmichael number: 983∙3929∙8839 983∙8839∙1241249 983∙970217∙190744663 Carmichael number: 991∙4951∙58411 991∙10111∙501001 991∙16831∙26731 991∙56431∙607861 991∙99991∙5215321 991∙118801∙206911 991∙138403∙336997 991∙167311∙312841 991∙338581∙890011 991∙658351∙1924561 Carmichael number: 997∙1993∙56773 997∙8467∙367027 997∙12451∙4137883 997∙17929∙130477 997∙29383∙450691 997∙167329∙15166093 997∙1002973∙99996409 ──────── 1038 Carmichael numbers found.
Ruby
<lang ruby># 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 puts
end</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
Rust
<lang rust> fn is_prime(n: i64) -> bool {
if n > 1 { (2..((n / 2) + 1)).all(|x| n % x != 0) } else { false }
}
// The modulo operator actually calculates the remainder. fn modulo(n: i64, m: i64) -> i64 {
((n % m) + m) % m
}
fn carmichael(p1: i64) -> Vec<(i64, i64, i64)> {
let mut results = Vec::new(); if !is_prime(p1) { return results; }
for h3 in 2..p1 { for d in 1..(h3 + p1) { if (h3 + p1) * (p1 - 1) % d != 0 || modulo(-p1 * p1, h3) != d % h3 { continue; }
let p2 = 1 + ((p1 - 1) * (h3 + p1) / d); if !is_prime(p2) { continue; }
let p3 = 1 + (p1 * p2 / h3); if !is_prime(p3) || ((p2 * p3) % (p1 - 1) != 1) { continue; }
results.push((p1, p2, p3)); } }
results
}
fn main() {
(1..62) .filter(|&x| is_prime(x)) .map(carmichael) .filter(|x| !x.is_empty()) .flat_map(|x| x) .inspect(|x| println!("{:?}", x)) .count(); // Evaluate entire iterator
} </lang>
- Output:
(3, 11, 17) (5, 29, 73) (5, 17, 29) (5, 13, 17) . . . (61, 661, 2521) (61, 271, 571) (61, 241, 421) (61, 3361, 4021)
Seed7
The function isPrime below is borrowed from the Seed7 algorithm collection.
<lang seed7>$ 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;</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
Sidef
<lang ruby>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);</lang>
- 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... <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
zkl
Using the Miller-Rabin primality test in lib GMP. <lang zkl>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 }</lang> <lang>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) });</lang>
- 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
ZX Spectrum Basic
<lang zxbasic>10 FOR p=2 TO 61 20 LET n=p: GO SUB 1000 30 IF NOT n THEN GO TO 200 40 FOR h=1 TO p-1 50 FOR d=1 TO h-1+p 60 IF NOT (FN m((h+p)*(p-1),d)=0 AND FN w(-p*p,h)=FN m(d,h)) THEN GO TO 180 70 LET q=INT (1+((p-1)*(h+p)/d)) 80 LET n=q: GO SUB 1000 90 IF NOT n THEN GO TO 180 100 LET r=INT (1+(p*q/h)) 110 LET n=r: GO SUB 1000 120 IF (NOT n) OR ((FN m((q*r),(p-1))<>1)) THEN GO TO 180 130 PRINT p;" ";q;" ";r 180 NEXT d 190 NEXT h 200 NEXT p 210 STOP 1000 IF n<4 THEN LET n=(n>1): RETURN 1010 IF (NOT FN m(n,2)) OR (NOT FN m(n,3)) THEN LET n=0: RETURN 1020 LET i=5 1030 IF NOT ((i*i)<=n) THEN LET n=1: RETURN 1040 IF (NOT FN m(n,i)) OR NOT FN m(n,(i+2)) THEN LET n=0: RETURN 1050 LET i=i+6 1060 GO TO 1030 2000 DEF FN m(a,b)=a-(INT (a/b)*b): REM Mod function 2010 DEF FN w(a,b)=FN m(FN m(a,b)+b,b): REM Mod function modified </lang>