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.
- 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, 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 Stromg Pseudoprimes for Prime1 up to 61 #
FOR prime1 FROM 2 TO 61 DO
IF is prime[ prime 1 ] THEN 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 FI
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 13 61 397 13 37 241 13 97 421 13 37 97 13 37 61 ... 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
Clojure
<lang lisp> (ns example
(:gen-class))
(defn prime? [n]
" Prime number test (using Java) " (.isProbablePrime (biginteger n) 16))
(defn carmichael [p1]
" Triplets of Carmichael primes, with first element prime p1 " (if (prime? p1) (into [] (for [h3 (range 2 p1) :let [g (+ h3 p1)] d (range 1 g) :when (and (= (mod (* g (dec p1)) d) 0) (= (mod (- (* p1 p1)) h3) (mod d h3))) :let [p2 (inc (quot (* (dec p1) g) d))] :when (prime? p2) :let [p3 (inc (quot (* p1 p2) h3))] :when (prime? p3) :when (= (mod (* p2 p3) (dec p1)) 1)] [p1 p2 p3]))))
- Generate Result
(def numbers (mapcat carmichael (range 2 62))) (println (count numbers) "Carmichael numbers found:") (doseq [t numbers]
(println (format "%5d x %5d x %5d = %10d" (first t) (second t) (last t) (apply * t))))
</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
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
Plan
This is F77 style, and directly translates the given calculation as per formula translation. 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.
Source
So, using the double MOD approach (see the Discussion) - 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
FreeBASIC
<lang freebasic>' version 17-10-2016 ' compile with: fbc -s console
' using a sieve for finding primes
- Define max_sieve 10000000 ' 10^7
ReDim Shared As Byte isprime(max_sieve)
' translated the pseudo code to FreeBASIC Sub carmichael3(p1 As Integer)
If isprime(p1) = 0 Then Exit Sub
Dim As Integer h3, d, p2, p3, t1, t2
For h3 = 1 To p1 -1 t1 = (h3 + p1) * (p1 -1) t2 = (-p1 * p1) Mod h3 If t2 < 0 Then t2 = t2 + h3 For d = 1 To h3 + p1 -1 If t1 Mod d = 0 And t2 = (d Mod h3) Then p2 = 1 + (t1 \ d) If isprime(p2) = 0 Then Continue For p3 = 1 + (p1 * p2 \ h3) If isprime(p3) = 0 Or ((p2 * p3) Mod (p1 -1)) <> 1 Then Continue For Print Using "### * #### * #####"; p1; p2; p3 End If Next d Next h3
End Sub
' ------=< MAIN >=------
Dim As UInteger i, j
'set up sieve For i = 3 To max_sieve Step 2
isprime(i) = 1
Next i
isprime(2) = 1 For i = 3 To Sqr(max_sieve) Step 2
If isprime(i) = 1 Then For j = i * i To max_sieve Step i * 2 isprime(j) = 0 Next j End If
Next i
For i = 2 To 61
carmichael3(i)
Next i
' empty keyboard buffer While InKey <> "" : Wend Print : Print "hit any key to end program" Sleep End</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
Go
<lang go>package main
import "fmt"
// Use this rather than % for negative integers func mod(n, m int) int {
return ((n % m) + m) % m
}
func isPrime(n int) bool {
if n < 2 { return false } if n % 2 == 0 { return n == 2 } if n % 3 == 0 { return n == 3 } d := 5 for d * d <= n { if n % d == 0 { return false } d += 2 if n % d == 0 { return false } d += 4 } return true
}
func carmichael(p1 int) {
for h3 := 2; 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 !isPrime(p2) { continue } p3 := 1 + p1 * p2 / h3 if !isPrime(p3) { continue } if p2 * p3 % (p1 - 1) != 1 { continue } c := p1 * p2 * p3 fmt.Printf("%2d %4d %5d %d\n", p1, p2, p3, c) } } }
}
func main() {
fmt.Println("The following are Carmichael munbers for p1 <= 61:\n") fmt.Println("p1 p2 p3 product") fmt.Println("== == == =======")
for p1 := 2; p1 <= 61; p1++ { if isPrime(p1) { carmichael(p1) } }
}</lang>
- Output:
The following are Carmichael munbers for p1 <= 61: p1 p2 p3 product == == == ======= 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 ->
J
<lang J> q =: (,"0 1~ >:@i.@<:@+/"1)&.>@(,&.>"0 1~ >:@i.)&.>@I.@(1&p:@i.)@>: f1 =: (0: = {. | <:@{: * 1&{ + {:) *. ((1&{ | -@*:@{:) = 1&{ | {.) f2 =: 1: = <:@{. | ({: * 1&{) p2 =: 0:`((* 1&p:)@(<.@(1: + <:@{: * {. %~ 1&{ + {:)))@.f1 p3 =: 3:$0:`((* 1&p:)@({: , {. , (<.@>:@(1&{ %~ {. * {:))))@.(*@{.)@(p2 , }.) (-. 3:$0:)@(((*"0 f2)@p3"1)@;@;@q) 61 </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
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.
Kotlin
<lang scala>fun Int.isPrime(): Boolean {
return when { this == 2 -> true this <= 1 || this % 2 == 0 -> false else -> { val max = Math.sqrt(toDouble()).toInt() (3..max step 2) .filter { this % it == 0 } .forEach { return false } true } }
}
fun mod(n: Int, m: Int) = ((n % m) + m) % m
fun main(args: Array<String>) {
for (p1 in 3..61) { if (p1.isPrime()) { for (h3 in 2 until p1) { val g = h3 + p1 for (d in 1 until g) { if ((g * (p1 - 1)) % d == 0 && mod(-p1 * p1, h3) == d % h3) { val q = 1 + (p1 - 1) * g / d if (q.isPrime()) { val r = 1 + (p1 * q / h3) if (r.isPrime() && (q * r) % (p1 - 1) == 1) { println("$p1 x $q x $r") } } } } } } }
}</lang>
- Output:
See D output.
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>func forprimes(a, b, callback) {
for (a = (a-1 -> next_prime); a <= b; a.next_prime!) { callback(a) }
}
forprimes(3, 61, func(p) {
for h3 in (2 ..^ p) { var ph3 = (p + h3) for d in (1 ..^ ph3) { ((-p * p) % h3) != (d % h3) && next ((p-1) * ph3) % d && next var q = 1+((p-1) * ph3 / d) q.is_prime || next var r = 1+((p*q - 1)/h3) r.is_prime || next (q*r) % (p-1) == 1 || next printf("%2d x %5d x %5d = %s\n",p,q,r, p*q*r) } }
})</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>