Jump to content

Chernick's Carmichael numbers

From Rosetta Code
Revision as of 07:56, 15 September 2024 by Zeddicus (talk | contribs) ({{header|REXX}})
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Task
Chernick's Carmichael numbers
You are encouraged to solve this task according to the task description, using any language you may know.

In 1939, Jack Chernick proved that, for n ≥ 3 and m ≥ 1:

   U(n, m) = (6m + 1) * (12m + 1) * Product_{i=1..n-2} (2^i * 9m + 1)

is a Carmichael number if all the factors are primes and, for n > 4, m is a multiple of 2^(n-4).


Example
   U(3, m) = (6m + 1) * (12m + 1) * (18m + 1)
   U(4, m) = U(3, m) * (2^2 * 9m + 1)
   U(5, m) = U(4, m) * (2^3 * 9m + 1)
   ...
   U(n, m) = U(n-1, m) * (2^(n-2) * 9m + 1)
  • The smallest Chernick's Carmichael number with 3 prime factors, is: U(3, 1) = 1729.
  • The smallest Chernick's Carmichael number with 4 prime factors, is: U(4, 1) = 63973.
  • The smallest Chernick's Carmichael number with 5 prime factors, is: U(5, 380) = 26641259752490421121.


For n = 5, the smallest number m that satisfy Chernick's conditions, is m = 380, therefore U(5, 380) is the smallest Chernick's Carmichael number with 5 prime factors.

U(5, 380) is a Chernick's Carmichael number because m = 380 is a multiple of 2^(n-4), where n = 5, and the factors { (6*380 + 1), (12*380 + 1), (18*380 + 1), (36*380 + 1), (72*380 + 1) } are all prime numbers.


Task

For n ≥ 3, let a(n) be the smallest Chernick's Carmichael number with n prime factors.

  • Compute a(n) for n = 3..9.
  • Optional: find a(10).


Note: it's perfectly acceptable to show the terms in factorized form:

 a(3) = 7 * 13 * 19
 a(4) = 7 * 13 * 19 * 37
 a(5) = 2281 * 4561 * 6841 * 13681 * 27361
 ...


See also


Related tasks



C

Library: GMP
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

typedef unsigned long long int u64;

#define TRUE 1
#define FALSE 0

int primality_pretest(u64 k) {
    if (!(k %  3) || !(k %  5) || !(k %  7) || !(k % 11) || !(k % 13) || !(k % 17) || !(k % 19) || !(k % 23)) return (k <= 23);
    return TRUE;
}

int probprime(u64 k, mpz_t n) {
    mpz_set_ui(n, k);
    return mpz_probab_prime_p(n, 0);
}

int is_chernick(int n, u64 m, mpz_t z) {
    u64 t = 9 * m;
    if (primality_pretest(6 * m + 1) == FALSE) return FALSE;
    if (primality_pretest(12 * m + 1) == FALSE) return FALSE;
    for (int i = 1; i <= n - 2; i++) if (primality_pretest((t << i) + 1) == FALSE) return FALSE;
    if (probprime(6 * m + 1, z) == FALSE) return FALSE;
    if (probprime(12 * m + 1, z) == FALSE) return FALSE;
    for (int i = 1; i <= n - 2; i++) if (probprime((t << i) + 1, z) == FALSE) return FALSE;
    return TRUE;
}

int main(int argc, char const *argv[]) {
    mpz_t z;
    mpz_inits(z, NULL);

    for (int n = 3; n <= 10; n ++) {
        u64 multiplier = (n > 4) ? (1 << (n - 4)) : 1;

        if (n > 5) multiplier *= 5;

        for (u64 k = 1; ; k++) {
            u64 m = k * multiplier;

            if (is_chernick(n, m, z) == TRUE) {
                printf("a(%d) has m = %llu\n", n, m);
                break;
            }
        }
    }

    return 0;
}
Output:
a(3) has m = 1
a(4) has m = 1
a(5) has m = 380
a(6) has m = 380
a(7) has m = 780320
a(8) has m = 950560
a(9) has m = 950560
a(10) has m = 3208386195840

C++

Library: GMP
#include <gmp.h>
#include <iostream>

using namespace std;

typedef unsigned long long int u64;

bool primality_pretest(u64 k) {     // for k > 23

    if (!(k %  3) || !(k %  5) || !(k %  7) || !(k % 11) ||
        !(k % 13) || !(k % 17) || !(k % 19) || !(k % 23)
    ) {
        return (k <= 23);
    }

    return true;
}

bool probprime(u64 k, mpz_t n) {
    mpz_set_ui(n, k);
    return mpz_probab_prime_p(n, 0);
}

bool is_chernick(int n, u64 m, mpz_t z) {

    if (!primality_pretest(6 * m + 1)) {
        return false;
    }

    if (!primality_pretest(12 * m + 1)) {
        return false;
    }

    u64 t = 9 * m;

    for (int i = 1; i <= n - 2; i++) {
        if (!primality_pretest((t << i) + 1)) {
            return false;
        }
    }

    if (!probprime(6 * m + 1, z)) {
        return false;
    }

    if (!probprime(12 * m + 1, z)) {
        return false;
    }

    for (int i = 1; i <= n - 2; i++) {
        if (!probprime((t << i) + 1, z)) {
            return false;
        }
    }

    return true;
}

int main() {

    mpz_t z;
    mpz_inits(z, NULL);

    for (int n = 3; n <= 10; n++) {

        // `m` is a multiple of 2^(n-4), for n > 4
        u64 multiplier = (n > 4) ? (1 << (n - 4)) : 1;

        // For n > 5, m is also a multiple of 5
        if (n > 5) {
            multiplier *= 5;
        }

        for (u64 k = 1; ; k++) {

            u64 m = k * multiplier;

            if (is_chernick(n, m, z)) {
                cout << "a(" << n << ") has m = " << m << endl;
                break;
            }
        }
    }

    return 0;
}
Output:
a(3) has m = 1
a(4) has m = 1
a(5) has m = 380
a(6) has m = 380
a(7) has m = 780320
a(8) has m = 950560
a(9) has m = 950560
a(10) has m = 3208386195840

(takes ~3.5 minutes)

F#

// Generate Chernick's Carmichael numbers. Nigel Galloway: June 1st., 2019
open Open.Numeric.Primes
let fMk m k=Number.IsPrime(6UL*m+1UL) && Number.IsPrime(12UL*m+1UL) && [1..k-2]|>List.forall(fun n->Number.IsPrime(9UL*(pown 2UL n)*m+1UL))
let fX k=Seq.initInfinite(fun n->match k with 3->uint64(n+1) |_->uint64(n+1)*(pown 2UL (k-4)))|>Seq.filter(fun n->fMk n k)
let cherCar k=let m=Seq.head(fX k) in sprintf "m=%d primes-> %A " m ([6UL*m+1UL;12UL*m+1UL]@List.init(k-2)(fun n->9UL*(pown 2UL (n+1))*m+1UL))
[3..9] |> Seq.iter(fun g->printfn $"cherCar %d{g}: %s{cherCar g}")
Output:
cherCar 3: m=1 primes-> [7; 13; 19]
cherCar 4: m=1 primes-> [7; 13; 19; 37]
cherCar 5: m=380 primes-> [2281; 4561; 6841; 13681; 27361]
cherCar 6: m=380 primes-> [2281; 4561; 6841; 13681; 27361; 54721]
cherCar 7: m=780320 primes-> [4681921; 9363841; 14045761; 28091521; 56183041; 112366081; 224732161]
cherCar 8: m=950560 primes-> [5703361; 11406721; 17110081; 34220161; 68440321; 136880641; 273761281;
 547522561]
cherCar 9: m=950560 primes-> [5703361; 11406721; 17110081; 34220161; 68440321; 136880641; 273761281;
 547522561; 1095045121]

FreeBASIC

Basic only

#include "isprime.bas"

Function PrimalityPretest(k As Integer) As Boolean
    Dim As Integer ppp(1 To 8) = {3,5,7,11,13,17,19,23}
    For i As Integer = 1 To Ubound(ppp)
        If k Mod ppp(i) = 0 Then Return (k <= 23)
    Next i
    Return True
End Function

Function isChernick(n As Integer, m As Integer) As Boolean
    Dim As Integer i, t = 9 * m
    
    If Not PrimalityPretest(6 * m + 1) Then Return False
    If Not PrimalityPretest(12 * m + 1) Then Return False
    
    For i = 1 To n-1
        If Not PrimalityPretest(t * (2 ^ i) + 1) Then Return False
    Next i
    
    If Not isPrime(6 * m + 1) Then Return False
    If Not isPrime(12 * m + 1) Then Return False
    
    For i = 1 To n - 2
        If Not isPrime(t * (2 ^ i) + 1) Then Return False
    Next i
    Return True
End Function

Dim As Uinteger multiplier, k, m = 1
For n As Integer = 3 To 9
    multiplier = Iif (n > 4, 2 ^ (n-4), 1)
    If n > 5 Then multiplier *= 5
    
    k = 1
    Do
        m = k * multiplier
        If isChernick(n, m) Then
            Print "a(" & n & ") has m = " & m
            Exit Do
        End If
        k += 1
    Loop
Next n
Sleep

Go

Basic only

package main

import (
    "fmt"
    "math/big"
)

var (
    zero = new(big.Int)
    prod = new(big.Int)
    fact = new(big.Int)
)

func ccFactors(n, m uint64) (*big.Int, bool) {
    prod.SetUint64(6*m + 1)
    if !prod.ProbablyPrime(0) {
        return zero, false
    }
    fact.SetUint64(12*m + 1)
    if !fact.ProbablyPrime(0) { // 100% accurate up to 2 ^ 64
        return zero, false
    }
    prod.Mul(prod, fact)
    for i := uint64(1); i <= n-2; i++ {
        fact.SetUint64((1<<i)*9*m + 1)
        if !fact.ProbablyPrime(0) {
            return zero, false
        }
        prod.Mul(prod, fact)
    }
    return prod, true
}

func ccNumbers(start, end uint64) {
    for n := start; n <= end; n++ {
        m := uint64(1)
        if n > 4 {
            m = 1 << (n - 4)
        }
        for {
            num, ok := ccFactors(n, m)
            if ok {
                fmt.Printf("a(%d) = %d\n", n, num)
                break
            }
            if n <= 4 {
                m++
            } else {
                m += 1 << (n - 4)
            }
        }
    }
}

func main() {
    ccNumbers(3, 9)
}
Output:
a(3) = 1729
a(4) = 63973
a(5) = 26641259752490421121
a(6) = 1457836374916028334162241
a(7) = 24541683183872873851606952966798288052977151461406721
a(8) = 53487697914261966820654105730041031613370337776541835775672321
a(9) = 58571442634534443082821160508299574798027946748324125518533225605795841

Basic plus optional


To reach a(10) in a reasonable time, a much more efficient approach is needed.

The following version takes account of the optimizations referred to in the Talk page and previewed in the C++ entry above.

It also uses a wrapper for the C library, GMP, which despite the overhead of cgo is still much faster than Go's native big.Int library.

The resulting executable is several hundred times faster than before and, even on my modest Celeron @1.6GHZ, reaches a(9) in under 10ms and a(10) in about 22 minutes.

package main

import (
    "fmt"
    big "github.com/ncw/gmp"
)

const (
    min = 3
    max = 10
)

var (
    prod       = new(big.Int)
    fact       = new(big.Int)
    factors    = [max]uint64{}
    bigFactors = [max]*big.Int{}
)

func init() {
    for i := 0; i < max; i++ {
        bigFactors[i] = big.NewInt(0)
    }
}

func isPrimePretest(k uint64) bool {
    if k%3 == 0 || k%5 == 0 || k%7 == 0 || k%11 == 0 ||
        k%13 == 0 || k%17 == 0 || k%19 == 0 || k%23 == 0 {
        return k <= 23
    }
    return true
}

func ccFactors(n, m uint64) bool {
    if !isPrimePretest(6*m + 1) {
        return false
    }
    if !isPrimePretest(12*m + 1) {
        return false
    }
    factors[0] = 6*m + 1
    factors[1] = 12*m + 1
    t := 9 * m
    for i := uint64(1); i <= n-2; i++ {
        tt := (t << i) + 1
        if !isPrimePretest(tt) {
            return false
        }
        factors[i+1] = tt
    }

    for i := 0; i < int(n); i++ {
        fact.SetUint64(factors[i])
        if !fact.ProbablyPrime(0) {
            return false
        }
        bigFactors[i].Set(fact)
    }
    return true
}

func prodFactors(n uint64) *big.Int {
    prod.Set(bigFactors[0])
    for i := 1; i < int(n); i++ {
        prod.Mul(prod, bigFactors[i])
    }
    return prod
}

func ccNumbers(start, end uint64) {
    for n := start; n <= end; n++ {
        mult := uint64(1)
        if n > 4 {
            mult = 1 << (n - 4)
        }
        if n > 5 {
            mult *= 5
        }
        m := mult
        for {
            if ccFactors(n, m) {
                num := prodFactors(n)
                fmt.Printf("a(%d) = %d\n", n, num)
                fmt.Printf("m(%d) = %d\n", n, m)
                fmt.Println("Factors:", factors[:n], "\n")
                break
            }
            m += mult
        }
    }
}

func main() {
    ccNumbers(min, max)
}
Output:
a(3) = 1729
m(3) = 1
Factors: [7 13 19] 

a(4) = 63973
m(4) = 1
Factors: [7 13 19 37] 

a(5) = 26641259752490421121
m(5) = 380
Factors: [2281 4561 6841 13681 27361] 

a(6) = 1457836374916028334162241
m(6) = 380
Factors: [2281 4561 6841 13681 27361 54721] 

a(7) = 24541683183872873851606952966798288052977151461406721
m(7) = 780320
Factors: [4681921 9363841 14045761 28091521 56183041 112366081 224732161] 

a(8) = 53487697914261966820654105730041031613370337776541835775672321
m(8) = 950560
Factors: [5703361 11406721 17110081 34220161 68440321 136880641 273761281 547522561] 

a(9) = 58571442634534443082821160508299574798027946748324125518533225605795841
m(9) = 950560
Factors: [5703361 11406721 17110081 34220161 68440321 136880641 273761281 547522561 1095045121] 

a(10) = 24616075028246330441656912428380582403261346369700917629170235674289719437963233744091978433592331048416482649086961226304033068172880278517841921
m(10) = 3208386195840
Factors: [19250317175041 38500634350081 57750951525121 115501903050241 231003806100481 462007612200961 924015224401921 1848030448803841 3696060897607681 7392121795215361] 

J

Brute force:

a=: {{)v
  if.3=y do.1729 return.end.
  m=. z=. 2^y-4
  f=. 6 12,9*2^}.i.y-1
  while.do.
    uf=.1+f*m
    if.*/1 p: uf do. */x:uf return.end.
    m=.m+z
  end.
}}

Task examples:

   a 3
1729
   a 4
63973
   a 5
26641259752490421121
   a 6
1457836374916028334162241
   a 7
24541683183872873851606952966798288052977151461406721
   a 8
53487697914261966820654105730041031613370337776541835775672321
   a 9
58571442634534443082821160508299574798027946748324125518533225605795841

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class ChernicksCarmichaelNumbers {

    public static void main(String[] args) {
        for ( long n = 3 ; n < 10 ; n++ ) {
            long m = 0;
            boolean foundComposite = true;
            List<Long> factors = null;
            while ( foundComposite ) {
                m += (n <= 4 ? 1 : (long) Math.pow(2, n-4) * 5);
                factors = U(n, m);
                foundComposite = false;
                for ( long factor : factors ) {
                    if ( ! isPrime(factor) ) {
                        foundComposite = true;
                        break;
                    }
                }
            }
            System.out.printf("U(%d, %d) = %s = %s %n", n, m, display(factors), multiply(factors));
        }
    }
    
    private static String display(List<Long> factors) {
        return factors.toString().replace("[", "").replace("]", "").replaceAll(", ", " * ");
    }
    
    private static BigInteger multiply(List<Long> factors) {
        BigInteger result = BigInteger.ONE;
        for ( long factor : factors ) {
            result = result.multiply(BigInteger.valueOf(factor));
        }
        return result;
    }
    
    private static List<Long> U(long n, long m) {
        List<Long> factors = new ArrayList<>();
        factors.add(6*m + 1);
        factors.add(12*m + 1);
        for ( int i = 1 ; i <= n-2 ; i++ ) {
            factors.add(((long)Math.pow(2, i)) * 9 * m + 1);
        }
        return factors;
    }

    private static final int MAX = 100_000;
    private static final boolean[] primes = new boolean[MAX];
    private static boolean SIEVE_COMPLETE = false;
    
    private static final boolean isPrimeTrivial(long test) {
        if ( ! SIEVE_COMPLETE ) {
            sieve();
            SIEVE_COMPLETE = true;
        }
        return primes[(int) test];
    }
    
    private static final void sieve() {
        //  primes
        for ( int i = 2 ; i < MAX ; i++ ) {
            primes[i] = true;            
        }
        for ( int i = 2 ; i < MAX ; i++ ) {
            if ( primes[i] ) {
                for ( int j = 2*i ; j < MAX ; j += i ) {
                    primes[j] = false;
                }
            }
        }
    }

    //  See http://primes.utm.edu/glossary/page.php?sort=StrongPRP
    public static final boolean isPrime(long testValue) {
        if ( testValue == 2 ) return true;
        if ( testValue % 2 == 0 ) return false;
        if ( testValue <= MAX ) return isPrimeTrivial(testValue);
        long d = testValue-1;
        int s = 0;
        while ( d % 2 == 0 ) {
            s += 1;
            d /= 2;
        }
        if ( testValue < 1373565L ) {
            if ( ! aSrp(2, s, d, testValue) ) {
                return false;
            }
            if ( ! aSrp(3, s, d, testValue) ) {
                return false;
            }
            return true;
        }
        if ( testValue < 4759123141L ) {
            if ( ! aSrp(2, s, d, testValue) ) {
                return false;
            }
            if ( ! aSrp(7, s, d, testValue) ) {
                return false;
            }
            if ( ! aSrp(61, s, d, testValue) ) {
                return false;
            }
            return true;
        }
        if ( testValue < 10000000000000000L ) {
            if ( ! aSrp(3, s, d, testValue) ) {
                return false;
            }
            if ( ! aSrp(24251, s, d, testValue) ) {
                return false;
            }
            return true;
        }
        //  Try 5 "random" primes
        if ( ! aSrp(37, s, d, testValue) ) {
            return false;
        }
        if ( ! aSrp(47, s, d, testValue) ) {
            return false;
        }
        if ( ! aSrp(61, s, d, testValue) ) {
            return false;
        }
        if ( ! aSrp(73, s, d, testValue) ) {
            return false;
        }
        if ( ! aSrp(83, s, d, testValue) ) {
            return false;
        }
        //throw new RuntimeException("ERROR isPrime:  Value too large = "+testValue);
        return true;
    }

    private static final boolean aSrp(int a, int s, long d, long n) {
        long modPow = modPow(a, d, n);
        //System.out.println("a = "+a+", s = "+s+", d = "+d+", n = "+n+", modpow = "+modPow);
        if ( modPow == 1 ) {
            return true;
        }
        int twoExpR = 1;
        for ( int r = 0 ; r < s ; r++ ) {
            if ( modPow(modPow, twoExpR, n) == n-1 ) {
                return true;
            }
            twoExpR *= 2;
        }
        return false;
    }
    
    private static final long SQRT = (long) Math.sqrt(Long.MAX_VALUE);
    
    public static final long modPow(long base, long exponent, long modulus) {
        long result = 1;
        while ( exponent > 0 ) {
            if ( exponent % 2 == 1 ) {
                if ( result > SQRT || base > SQRT ) {
                    result = multiply(result, base, modulus);
                }
                else {
                    result = (result * base) % modulus;
                }
            }
            exponent >>= 1;
            if ( base > SQRT ) {
                base = multiply(base, base, modulus);
            }
            else {
                base = (base * base) % modulus;
            }
        }
        return result;
    }


    //  Result is a*b % mod, without overflow.
    public static final long multiply(long a, long b, long modulus) {
        long x = 0;
        long y = a % modulus;
        long t;
        while ( b > 0 ) {
            if ( b % 2 == 1 ) {
                t = x + y;
                x = (t > modulus ? t-modulus : t);
            }
            t = y << 1;
            y = (t > modulus ? t-modulus : t);
            b >>= 1;
        }
        return x % modulus;
    }

}
Output:
U(3, 1) = 7 * 13 * 19 = 1729 
U(4, 1) = 7 * 13 * 19 * 37 = 63973 
U(5, 380) = 2281 * 4561 * 6841 * 13681 * 27361 = 26641259752490421121 
U(6, 380) = 2281 * 4561 * 6841 * 13681 * 27361 * 54721 = 1457836374916028334162241 
U(7, 780320) = 4681921 * 9363841 * 14045761 * 28091521 * 56183041 * 112366081 * 224732161 = 24541683183872873851606952966798288052977151461406721 
U(8, 950560) = 5703361 * 11406721 * 17110081 * 34220161 * 68440321 * 136880641 * 273761281 * 547522561 = 53487697914261966820654105730041031613370337776541835775672321 
U(9, 950560) = 5703361 * 11406721 * 17110081 * 34220161 * 68440321 * 136880641 * 273761281 * 547522561 * 1095045121 = 58571442634534443082821160508299574798027946748324125518533225605795841 

Julia

using Primes

function trial_pretest(k::UInt64)

    if ((k %  3)==0 || (k %  5)==0 || (k %  7)==0 || (k % 11)==0 ||
        (k % 13)==0 || (k % 17)==0 || (k % 19)==0 || (k % 23)==0)
        return (k <= 23)
    end

    return true
end

function gcd_pretest(k::UInt64)

    if (k <= 107)
        return true
    end

    gcd(29*31*37*41*43*47*53*59*61*67, k) == 1 &&
    gcd(71*73*79*83*89*97*101*103*107, k) == 1
end

function is_chernick(n::Int64, m::UInt64)

    t = 9*m

    if (!trial_pretest(6*m + 1))
        return false
    end

    if (!trial_pretest(12*m + 1))
        return false
    end

    for i in 1:n-2
        if (!trial_pretest((t << i) + 1))
            return false
        end
    end

    if (!gcd_pretest(6*m + 1))
        return false
    end

    if (!gcd_pretest(12*m + 1))
        return false
    end

    for i in 1:n-2
        if (!gcd_pretest((t << i) + 1))
            return false
        end
    end

    if (!isprime(6*m + 1))
        return false
    end

    if (!isprime(12*m + 1))
        return false
    end

    for i in 1:n-2
        if (!isprime((t << i) + 1))
            return false
        end
    end

    return true
end

function chernick_carmichael(n::Int64, m::UInt64)
    prod = big(1)

    prod *= 6*m + 1
    prod *= 12*m + 1

    for i in 1:n-2
        prod *= ((big(9)*m)<<i) + 1
    end

    prod
end

function cc_numbers(from, to)

    for n in from:to

        multiplier = 1

        if (n > 4) multiplier = 1 << (n-4) end
        if (n > 5) multiplier *= 5 end

        m = UInt64(multiplier)

        while true

            if (is_chernick(n, m))
                println("a(", n, ") = ", chernick_carmichael(n, m))
                break
            end

            m += multiplier
        end
    end
end

cc_numbers(3, 10)
Output:
a(3) = 1729
a(4) = 63973
a(5) = 26641259752490421121
a(6) = 1457836374916028334162241
a(7) = 24541683183872873851606952966798288052977151461406721
a(8) = 53487697914261966820654105730041031613370337776541835775672321
a(9) = 58571442634534443082821160508299574798027946748324125518533225605795841
a(10) = 24616075028246330441656912428380582403261346369700917629170235674289719437963233744091978433592331048416482649086961226304033068172880278517841921

(takes ~6.5 minutes)

Mathematica / Wolfram Language

ClearAll[PrimeFactorCounts, U]
PrimeFactorCounts[n_Integer] := Total[FactorInteger[n][[All, 2]]]
U[n_, m_] := (6 m + 1) (12 m + 1) Product[2^i 9 m + 1, {i, 1, n - 2}]
FindFirstChernickCarmichaelNumber[n_Integer?Positive] := 
 Module[{step, i, m, formula, value},
  step = Ceiling[2^(n - 4)];
  If[n > 5, step *= 5];
  i = step;
  formula = U[n, m];
  PrintTemporary[Dynamic[i]];
  While[True,
   value = formula /. m -> i;
   If[PrimeFactorCounts[value] == n,
    Break[];
    ];
   i += step
   ];
  {i, value}
  ]
FindFirstChernickCarmichaelNumber[3]
FindFirstChernickCarmichaelNumber[4]
FindFirstChernickCarmichaelNumber[5]
FindFirstChernickCarmichaelNumber[6]
FindFirstChernickCarmichaelNumber[7]
FindFirstChernickCarmichaelNumber[8]
FindFirstChernickCarmichaelNumber[9]
Output:
{1,1729}
{1,63973}
{380,26641259752490421121}
{380,1457836374916028334162241}
{780320,24541683183872873851606952966798288052977151461406721}
{950560,53487697914261966820654105730041031613370337776541835775672321}
{950560,58571442634534443082821160508299574798027946748324125518533225605795841}

Nim

Library: bignum

Until a(9) a simple primality test using divisions by odd numbers is sufficient. But for a(10), it is necessary to improve the test. We have used here some optimizations found in other solutions:

– eliminating multiples of 3, 5, 7, 11, 13, 17, 19, 23;
– using a probability test which implies to use big integers; so, we have to convert the tested number to a big integer;
– for n >= 5, checking only values of m which are multiple of 5 (in fact, we check only the multiples of 5 × 2^(n-4).

With these optimizations, the program executes in 4-5 minutes.

import strutils, sequtils
import bignum

const
  Max = 10
  Factors: array[3..Max, int] = [1, 1, 2, 4, 8, 16, 32, 64]   # 1 for n=3 then 2^(n-4).
  FirstPrimes = [3, 5, 7, 11, 13, 17, 19, 23]

#---------------------------------------------------------------------------------------------------

iterator factors(n, m: Natural): Natural =
  ## Yield the factors of U(n, m).

  yield 6 * m + 1
  yield 12 * m + 1
  var k = 2
  for _ in 1..(n - 2):
    yield 9 * k * m + 1
    inc k, k

#---------------------------------------------------------------------------------------------------

proc mayBePrime(n: int): bool =
  ## First primality test.

  if n < 23: return true

  for p in FirstPrimes:
    if n mod p == 0:
      return false

  result = true

#---------------------------------------------------------------------------------------------------

proc isChernick(n, m: Natural): bool =
  ## Check if U(N, m) if a Chernick-Carmichael number.

  # Use the first and quick test.
  for factor in factors(n, m):
    if not factor.mayBePrime():
      return false

  # Use the slow probability test (need to use a big int).
  for factor in factors(n, m):
    if probablyPrime(newInt(factor), 25) == 0:
      return false

  result = true

#---------------------------------------------------------------------------------------------------

proc a(n: Natural): tuple[m: Natural, factors: seq[Natural]] =
  ## For a given "n", find the smallest Charnick-Carmichael number.

  var m: Natural = 0
  var incr = (if n >= 5: 5 else: 1) * Factors[n]  # For n >= 5, a(n) is a multiple of 5.

  while true:
    inc m, incr
    if isChernick(n, m):
      return (m, toSeq(factors(n, m)))

#———————————————————————————————————————————————————————————————————————————————————————————————————

import strformat

for n in 3..Max:
  let (m, factors) = a(n)

  stdout.write fmt"a({n}) = U({n}, {m}) = "
  var s = ""
  for factor in factors:
    s.addSep(" × ")
    s.add($factor)
  stdout.write s, '\n'
Output:
a(3) = U(3, 1) = 7 × 13 × 19
a(4) = U(4, 1) = 7 × 13 × 19 × 37
a(5) = U(5, 380) = 2281 × 4561 × 6841 × 13681 × 27361
a(6) = U(6, 380) = 2281 × 4561 × 6841 × 13681 × 27361 × 54721
a(7) = U(7, 780320) = 4681921 × 9363841 × 14045761 × 28091521 × 56183041 × 112366081 × 224732161
a(8) = U(8, 950560) = 5703361 × 11406721 × 17110081 × 34220161 × 68440321 × 136880641 × 273761281 × 547522561
a(9) = U(9, 950560) = 5703361 × 11406721 × 17110081 × 34220161 × 68440321 × 136880641 × 273761281 × 547522561 × 1095045121
a(10) = U(10, 3208386195840) = 19250317175041 × 38500634350081 × 57750951525121 × 115501903050241 × 231003806100481 × 462007612200961 × 924015224401921 × 1848030448803841 × 3696060897607681 × 7392121795215361

PARI/GP

cherCar(n)={
  my(C=vector(n));C[1]=6; C[2]=12; for(g=3,n,C[g]=2^(g-2)*9);
  my(i=1); my(N(g)=while(i<=n&ispseudoprime(g*C[i]+1),i=i+1); return(i>n));
     i=1;  my(G(g)=while(i<=n&isprime(g*C[i]+1),i=i+1); return(i>n));
  i=1; if(n>4,i=2^(n-4)); if(n>5,i=i*5); my(m=i); while(!(N(m)&G(m)),m=m+i);
  printf("cherCar(%d): m = %d\n",n,m)}
for(x=3,9,cherCar(x))
Output:
cherCar(3): m = 1
cherCar(4): m = 1
cherCar(5): m = 380
cherCar(6): m = 380
cherCar(7): m = 780320
cherCar(8): m = 950560
cherCar(9): m = 950560
cherCar(10): m = 3208386195840

Perl

Library: ntheory
use 5.020;
use warnings;
use ntheory qw/:all/;
use experimental qw/signatures/;

sub chernick_carmichael_factors ($n, $m) {
    (6*$m + 1, 12*$m + 1, (map { (1 << $_) * 9*$m + 1 } 1 .. $n-2));
}

sub chernick_carmichael_number ($n, $callback) {

    my $multiplier = ($n > 4) ? (1 << ($n-4)) : 1;

    for (my $m = 1 ; ; ++$m) {
        my @f = chernick_carmichael_factors($n, $m * $multiplier);
        next if not vecall { is_prime($_) } @f;
        $callback->(@f);
        last;
    }
}

foreach my $n (3..9) {
    chernick_carmichael_number($n, sub (@f) { say "a($n) = ", vecprod(@f) });
}
Output:
a(3) = 1729
a(4) = 63973
a(5) = 26641259752490421121
a(6) = 1457836374916028334162241
a(7) = 24541683183872873851606952966798288052977151461406721
a(8) = 53487697914261966820654105730041031613370337776541835775672321
a(9) = 58571442634534443082821160508299574798027946748324125518533225605795841

Phix

Library: Phix/mpfr
Translation of: Sidef
with javascript_semantics
function chernick_carmichael_factors(integer n, m)
    sequence res = {6*m + 1, 12*m + 1}
    for i=1 to n-2 do
        res &= power(2,i) * 9*m + 1
    end for
    return res
end function
 
include mpfr.e
mpz p = mpz_init()
 
function m_prime(atom a)
    mpz_set_d(p,a)
    return mpz_prime(p)
end function
 
function is_chernick_carmichael(integer n, m)
    return iff(n==2 ? m_prime(6*m + 1) and m_prime(12*m + 1)
                    : m_prime(power(2,n-2) * 9*m + 1) and 
                      is_chernick_carmichael(n-1, m))
end function
 
function chernick_carmichael_number(integer n)
    integer m = iff(n>4 ? power(2,n-4) : 1), mm = m
    while not is_chernick_carmichael(n, mm) do mm += m end while
    return {chernick_carmichael_factors(n, mm),mm}
end function
 
for n=3 to 9 do
    {sequence f, integer m} = chernick_carmichael_number(n)
    mpz_set_si(p,1)
    for i=1 to length(f) do
        mpz_mul_d(p,p,f[i])
        f[i] = sprintf("%d",f[i])
    end for
    printf(1,"U(%d,%d): %s = %s\n",{n,m,mpz_get_str(p),join(f," * ")})
end for
Output:
U(3,1): 1729 = 7 * 13 * 19
U(4,1): 63973 = 7 * 13 * 19 * 37
U(5,380): 26641259752490421121 = 2281 * 4561 * 6841 * 13681 * 27361
U(6,380): 1457836374916028334162241 = 2281 * 4561 * 6841 * 13681 * 27361 * 54721
U(7,780320): 24541683183872873851606952966798288052977151461406721 = 4681921 * 9363841 * 14045761 * 28091521 * 56183041 * 112366081 * 224732161
U(8,950560): 53487697914261966820654105730041031613370337776541835775672321 = 5703361 * 11406721 * 17110081 * 34220161 * 68440321 * 136880641 * 273761281 * 547522561
U(9,950560): 58571442634534443082821160508299574798027946748324125518533225605795841 = 5703361 * 11406721 * 17110081 * 34220161 * 68440321 * 136880641 * 273761281 * 547522561 * 1095045121

Pleasingly fast, note however that a(10) remains well out of reach / would probably need a complete rewrite.

with cheat

Translation of: C

with added cheat for the a(10) case - I found a nice big prime factor of k and added that on each iteration instead of 1.

You could also use the sequence {1,1,1,1,19,19,4877,457,457,12564169}, if you know a way to build that, and then it wouldn't be cheating anymore...

with javascript_semantics
include mpfr.e
sequence ppp = {3,5,7,11,13,17,19,23}
function primality_pretest(atom k)
    for i=1 to length(ppp) do
        if remainder(k,ppp[i])=0 then return (k<=23) end if
    end for
    return true
end function
 
function probprime(atom k, mpz n)
    mpz_set_d(n, k)
    return mpz_prime(n)
end function
 
function is_chernick(integer n, atom m, mpz z)
    atom t = 9 * m;
    if primality_pretest(6 * m + 1) == false then return false end if
    if primality_pretest(12 * m + 1) == false then return false end if
    for i=1 to n-3 do
        if primality_pretest(t*power(2,i) + 1) == false then return false end if
    end for
    if probprime(6 * m + 1, z) == false then return false end if
    if probprime(12 * m + 1, z) == false then return false end if
    for i=1 to n-2 do
        if probprime(t*power(2,i) + 1, z) == false then return false end if
    end for
    return true
end function
 
procedure main()
    atom t0 = time()
    mpz z = mpz_init(0)
 
    for n=3 to 10 do
        atom multiplier = iff(n>4 ? power(2,n-4) : 1), k = 1
        if n>5 then multiplier *= 5 end if
 
        while true do
            if n=10 then k += 12564168 end if   -- cheat!
            atom m = k * multiplier;
            if is_chernick(n, m, z) then
                printf(1,"a(%d) has m = %d\n", {n, m})
                exit
            end if
            k += 1
        end while
    end for
    ?elapsed(time()-t0)
end procedure
main()
Output:
a(3) has m = 1
a(4) has m = 1
a(5) has m = 380
a(6) has m = 380
a(7) has m = 780320
a(8) has m = 950560
a(9) has m = 950560
a(10) has m = 3208386195840
"0.1s"

Prolog

SWI Prolog is too slow to solve for a(10), even with optimizing by increasing the multiplier and implementing a trial division check. (actually, my implementation of Miller-Rabin in Prolog already starts with a trial division by small primes.)

?- use_module(library(primality)).

u(3, M, A * B * C) :-
    A is 6*M + 1, B is 12*M + 1, C is 18*M + 1, !.
u(N, M, U0 * D) :-
    succ(Pn, N), u(Pn, M, U0),
    D is 9*(1 << (N - 2))*M + 1.

prime_factorization(A*B) :- prime(B), prime_factorization(A), !.
prime_factorization(A) :- prime(A).

step(N, 1) :- N < 5, !.
step(5, 2) :- !.
step(N, K) :- K is 5*(1 << (N - 4)).

a(N, Factors) :- % due to backtracking nature of Prolog, a(n) will return all chernick-carmichael numbers.
    N > 2, !,
    step(N, I),
    between(1, infinite, J), M is I * J,
    u(N, M, Factors),
    prime_factorization(Factors).

main :-
    forall(
        (between(3, 9, K), once(a(K, Factorization)), N is Factorization),
        format("~w: ~w = ~w~n", [K, Factorization, N])),
    halt.

?- main.

isprime predicate:

prime(N) :-
    integer(N),
    N > 1,
    divcheck(
        N,
        [  2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31, 
          37,  41,  43,  47,  53,  59,  61,  67,  71,  73,  79,
          83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137,
         139, 149],
        Result),
    ((Result = prime, !); miller_rabin_primality_test(N)).

divcheck(_, [],    unknown) :- !.
divcheck(N, [P|_], prime) :- P*P > N, !.
divcheck(N, [P|Ps], State) :- N mod P =\= 0, divcheck(N, Ps, State).

miller_rabin_primality_test(N) :-
    bases(Bases, N),
    forall(member(A, Bases), strong_fermat_pseudoprime(N, A)).

miller_rabin_precision(16).

bases([31, 73], N) :- N < 9_080_191, !.
bases([2, 7, 61], N) :- N < 4_759_123_141, !.
bases([2, 325, 9_375, 28_178, 450_775, 9_780_504, 1_795_265_022], N) :-
    N < 18_446_744_073_709_551_616, !. % 2^64
bases(Bases, N) :-
    miller_rabin_precision(T), RndLimit is N - 2,
    length(Bases, T), maplist(random_between(2, RndLimit), Bases).

strong_fermat_pseudoprime(N, A) :-  % miller-rabin strong pseudoprime test with base A.
    succ(Pn, N), factor_2s(Pn, S, D),
    X is powm(A, D, N),
    ((X =:= 1, !); \+ composite_witness(N, S, X)).

composite_witness(_, 0, _) :- !.
composite_witness(N, K, X) :-
    X =\= N-1,
    succ(Pk, K), X2 is (X*X) mod N, composite_witness(N, Pk, X2).

factor_2s(N, S, D) :- factor_2s(0, N, S, D).
factor_2s(S, D, S, D) :- D /\ 1 =\= 0, !.
factor_2s(S0, D0, S, D) :-
    succ(S0, S1), D1 is D0 >> 1,
    factor_2s(S1, D1, S, D).
Output:
3: 7*13*19 = 1729
4: 7*13*19*37 = 63973
5: 2281*4561*6841*13681*27361 = 26641259752490421121
6: 2281*4561*6841*13681*27361*54721 = 1457836374916028334162241
7: 4681921*9363841*14045761*28091521*56183041*112366081*224732161 = 24541683183872873851606952966798288052977151461406721
8: 5703361*11406721*17110081*34220161*68440321*136880641*273761281*547522561 = 53487697914261966820654105730041031613370337776541835775672321
9: 5703361*11406721*17110081*34220161*68440321*136880641*273761281*547522561*1095045121 = 58571442634534443082821160508299574798027946748324125518533225605795841

Python

"""

Python implementation of
http://rosettacode.org/wiki/Chernick%27s_Carmichael_numbers

"""

# use sympy for prime test

from sympy import isprime

# based on C version

def primality_pretest(k):
    if not (k % 3) or not (k % 5) or not (k % 7) or not (k % 11) or not(k % 13) or not (k % 17) or not (k % 19) or not (k % 23):
        return (k <= 23)
        
    return True

def is_chernick(n, m):

    t = 9 * m
    
    if not primality_pretest(6 * m + 1):
        return False
        
    if not primality_pretest(12 * m + 1):
        return False
        
    for i in range(1,n-1):
        if not primality_pretest((t << i) + 1):
            return False
        
    if not isprime(6 * m + 1):
        return False
        
    if not isprime(12 * m + 1):
        return False
        
    for i in range(1,n - 1):
        if not isprime((t << i) + 1):
            return False
        
    return True
    
for n in range(3,10):

    if n > 4:
        multiplier = 1 << (n - 4)
    else:
        multiplier = 1
    
    if n > 5:
        multiplier *= 5
        
        
    k = 1
    
    while True:
        m = k * multiplier
        
        if is_chernick(n, m): 
            print("a("+str(n)+") has m = "+str(m))
            break
            
        k += 1
Output:
a(3) has m = 1
a(4) has m = 1
a(5) has m = 380
a(6) has m = 380
a(7) has m = 780320
a(8) has m = 950560
a(9) has m = 950560

Raku

(formerly Perl 6)

Use the ntheory library from Perl for primality testing since it is much, much faster than Raku's built-in .is-prime method.

Translation of: Perl
Library: ntheory
use Inline::Perl5;
use ntheory:from<Perl5> <:all>;

sub chernick-factors ($n, $m) {
    6×$m + 1, 12×$m + 1, |((1 .. $n-2).map: { (1 +< $_) × 9×$m + 1 } )
}

sub chernick-carmichael-number ($n) {

    my $multiplier = 1 +< (($n-4) max 0);
    my $iterator   = $n < 5 ?? (1 .. *) !! (1 .. *).map: * × 5;

    $multiplier × $iterator.first: -> $m {
        [&&] chernick-factors($n, $m × $multiplier).map: { is_prime($_) }
    }

}

for 3 .. 9 -> $n {
    my $m = chernick-carmichael-number($n);
    my @f = chernick-factors($n, $m);
    say "U($n, $m): {[×] @f} = {@f.join(' ⨉ ')}";
}
Output:
U(3, 1): 1729 = 7 ⨉ 13 ⨉ 19
U(4, 1): 63973 = 7 ⨉ 13 ⨉ 19 ⨉ 37
U(5, 380): 26641259752490421121 = 2281 ⨉ 4561 ⨉ 6841 ⨉ 13681 ⨉ 27361
U(6, 380): 1457836374916028334162241 = 2281 ⨉ 4561 ⨉ 6841 ⨉ 13681 ⨉ 27361 ⨉ 54721
U(7, 780320): 24541683183872873851606952966798288052977151461406721 = 4681921 ⨉ 9363841 ⨉ 14045761 ⨉ 28091521 ⨉ 56183041 ⨉ 112366081 ⨉ 224732161
U(8, 950560): 53487697914261966820654105730041031613370337776541835775672321 = 5703361 ⨉ 11406721 ⨉ 17110081 ⨉ 34220161 ⨉ 68440321 ⨉ 136880641 ⨉ 273761281 ⨉ 547522561
U(9, 950560): 58571442634534443082821160508299574798027946748324125518533225605795841 = 5703361 ⨉ 11406721 ⨉ 17110081 ⨉ 34220161 ⨉ 68440321 ⨉ 136880641 ⨉ 273761281 ⨉ 547522561 ⨉ 1095045121

REXX

Libraries: How to use
Library: Functions
Library: Numbers
Library: Settings
Library: Abend

Below program uses the improvements as given in Discussion and several other entries. The procedure Prime() is in library Numbers (Miller-Rabin with many tricks).

include Settings

say version; say 'Chernick''s Carmichael numbers'; say
numeric digits 80
say Copies('-',80)
say 'n   m(n) a(n)'
say Copies('-',80)
do n = 3 to 9
   mp = 1
   if n > 4 then
      mp = 2**(n-4)
   if n > 5 then
      mp = mp*5
   k = 0
   do x = 1
      k = k+1; m = mp*k; f.1 = 6*m+1
      if \ IsPrime(f.1) then
         iterate x
      f.2 = 12*m+1
      if \ IsPrime(f.2) then
         iterate x
      f = 2
      do i = 1 to n-2
         f = f+1; f.f = 2**i*9*m+1
         if \ IsPrime(f.f) then
            iterate x
      end
      a = 1
      do i = 1 to f
         a = a*f.i
      end
      say n Right(m,6) a
      leave x
   end
end
say Copies('-',80)
say Format(Time('e'),3,3) 'seconds'
say
exit

include Numbers
include Functions
include Abend
Output:
REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022
Chernick's Carmichael numbers

--------------------------------------------------------------------------------
n   m(n) a(n)
--------------------------------------------------------------------------------
3      1 1729
4      1 63973
5    380 26641259752490421121
6    380 1457836374916028334162241
7 780320 24541683183872873851606952966798288052977151461406721
8 950560 53487697914261966820654105730041031613370337776541835775672321
9 950560 58571442634534443082821160508299574798027946748324125518533225605795841
--------------------------------------------------------------------------------
  1.406 seconds

a(10) was not reachable with REXX. It would run for weeks.

Sidef

func chernick_carmichael_factors (n, m) {
    [6*m + 1, 12*m + 1, {|i| 2**i * 9*m + 1 }.map(1 .. n-2)...]
}

func is_chernick_carmichael (n, m) {
    (n == 2) ? (is_prime(6*m + 1) && is_prime(12*m + 1))
             : (is_prime(2**(n-2) * 9*m + 1) && __FUNC__(n-1, m))
}

func chernick_carmichael_number(n, callback) {
    var multiplier = (n>4 ? 2**(n-4) : 1)
    var m = (1..Inf -> first {|m| is_chernick_carmichael(n, m * multiplier) })
    var f = chernick_carmichael_factors(n, m * multiplier)
    callback(f...)
}

for n in (3..9) {
    chernick_carmichael_number(n, {|*f| say "a(#{n}) = #{f.join(' * ')}" })
}
Output:
a(3) = 7 * 13 * 19
a(4) = 7 * 13 * 19 * 37
a(5) = 2281 * 4561 * 6841 * 13681 * 27361
a(6) = 2281 * 4561 * 6841 * 13681 * 27361 * 54721
a(7) = 4681921 * 9363841 * 14045761 * 28091521 * 56183041 * 112366081 * 224732161
a(8) = 5703361 * 11406721 * 17110081 * 34220161 * 68440321 * 136880641 * 273761281 * 547522561
a(9) = 5703361 * 11406721 * 17110081 * 34220161 * 68440321 * 136880641 * 273761281 * 547522561 * 1095045121

Wren

Translation of: Go
Library: Wren-big
Library: Wren-fmt

Based on Go's 'more efficient' version. Reaches a(9) in just over 0.1 seconds but a(10) would still be out of reasonable reach for Wren so I've had to be content with that.

import "./big" for BigInt, BigInts
import "./fmt" for Fmt

var min = 3
var max = 9
var prod = BigInt.zero
var fact = BigInt.zero
var factors = List.filled(max, 0)
var bigFactors = List.filled(max, null)

var init = Fn.new {
    for (i in 0...max) bigFactors[i] = BigInt.zero
}

var isPrimePretest = Fn.new { |k|
    if (k%3 == 0 || k%5 == 0 || k%7 == 0 || k%11 == 0 ||
       (k%13 == 0) || k%17 == 0 || k%19 == 0 || k%23 == 0) return k <= 23
    return true
}

var ccFactors = Fn.new { |n, m|
    if (!isPrimePretest.call(6*m + 1)) return false
    if (!isPrimePretest.call(12*m + 1)) return false
    factors[0] = 6*m + 1
    factors[1] = 12*m + 1
    var t = 9 * m
    var i = 1
    while (i <= n-2) {
        var tt = (t << i) + 1
        if (!isPrimePretest.call(tt)) return false
        factors[i+1] = tt
        i = i + 1
    }
    for (i in 0...n) {
        fact = BigInt.new(factors[i])
        if (!fact.isProbablePrime(1)) return false
        bigFactors[i] = fact
    }
    return true
}

var ccNumbers = Fn.new { |start, end|
    for (n in start..end) {
        var mult = 1
        if (n > 4) mult = 1 << (n - 4)
        if (n > 5) mult = mult * 5
        var m = mult
        while (true) {
            if (ccFactors.call(n, m)) {
                var num = BigInts.prod(bigFactors.take(n))
                Fmt.print("a($d) = $i", n, num)
                Fmt.print("m($d) = $d", n, m)
                Fmt.print("Factors: $n\n", factors[0...n])
                break
            }
            m = m + mult
        }
    }
}

init.call()
ccNumbers.call(min, max)
Output:
a(3) = 1729
m(3) = 1
Factors: [7, 13, 19]

a(4) = 63973
m(4) = 1
Factors: [7, 13, 19, 37]

a(5) = 26641259752490421121
m(5) = 380
Factors: [2281, 4561, 6841, 13681, 27361]

a(6) = 1457836374916028334162241
m(6) = 380
Factors: [2281, 4561, 6841, 13681, 27361, 54721]

a(7) = 24541683183872873851606952966798288052977151461406721
m(7) = 780320
Factors: [4681921, 9363841, 14045761, 28091521, 56183041, 112366081, 224732161]

a(8) = 53487697914261966820654105730041031613370337776541835775672321
m(8) = 950560
Factors: [5703361, 11406721, 17110081, 34220161, 68440321, 136880641, 273761281, 547522561]

a(9) = 58571442634534443082821160508299574798027946748324125518533225605795841
m(9) = 950560
Factors: [5703361, 11406721, 17110081, 34220161, 68440321, 136880641, 273761281, 547522561, 1095045121]

zkl

Translation of: Go
Library: GMP

GNU Multiple Precision Arithmetic Library

Using GMP (probabilistic primes), because it is easy and fast to check primeness.

var [const] BI=Import("zklBigNum");  // libGMP

fcn ccFactors(n,m){	// not re-entrant
   prod:=BI(6*m + 1);
   if(not prod.probablyPrime())    return(False);
   fact:=BI(12*m + 1);
   if(not fact.probablyPrime())    return(False);
   prod.mul(fact);
   foreach i in ([1..n-2]){
      fact.set((2).pow(i) *9*m + 1);
      if(not fact.probablyPrime()) return(False);
      prod.mul(fact);
   }
   prod
}

fcn ccNumbers(start,end){
   foreach n in ([start..end]){
      a,m := ( if(n<=4) 1  else (2).pow(n - 4) ), a;
      while(1){
	 if(num := ccFactors(n,m)){
	    println("a(%d) = %,d".fmt(n,num));
	    break;
	 }
	 m+=a;
      }
   }
}
ccNumbers(3,9);
Output:
a(3) = 1,729
a(4) = 63,973
a(5) = 26,641,259,752,490,421,121
a(6) = 1,457,836,374,916,028,334,162,241
a(7) = 24,541,683,183,872,873,851,606,952,966,798,288,052,977,151,461,406,721
a(8) = 53,487,697,914,261,966,820,654,105,730,041,031,613,370,337,776,541,835,775,672,321
a(9) = 58,571,442,634,534,443,082,821,160,508,299,574,798,027,946,748,324,125,518,533,225,605,795,841
Cookies help us deliver our services. By using our services, you agree to our use of cookies.