7
\$\begingroup\$

Write a program to print Carmichael numbers. The challenge is to print the maximum number of them in less than 5 secs. The numbers should be in increasing order. Please do not use built in functions. Also please don't store the carmichael numbers in any form. To keep a common time measurement use Ideone. The program to print the maximum number of Carmichael number wins.

Sample Output
561
1105
1729
2465
2821
6601
8911
.
.
.

OEIS A002997

\$\endgroup\$
6
  • 3
    \$\begingroup\$ @JB: I for one usually accept the shortest answer only after a week or so. Also I doubt that this site will be that acceptance-heavy as SO. \$\endgroup\$
    – Joey
    Feb 12, 2011 at 18:05
  • \$\begingroup\$ 10 seconds seems an awfully short time for this problem. Unless you hardcode the answer... \$\endgroup\$ Feb 14, 2011 at 3:18
  • \$\begingroup\$ Is multi-threading allowed? \$\endgroup\$ Dec 22, 2016 at 19:56
  • \$\begingroup\$ I’m voting to close this question because online services aren't valid for fastest code due to connection and server differences. \$\endgroup\$
    – emanresu A
    Mar 30, 2022 at 4:05
  • \$\begingroup\$ @emanresuA This contradicts codegolf.meta.stackexchange.com/q/12707/46076, which says that they aren't a good choice for timing, but doesn't explicitly disallow them. (And I don't see any distinguishing factor between ideone and TIO) \$\endgroup\$ Mar 30, 2022 at 4:24

6 Answers 6

2
\$\begingroup\$

Python

This generates 1332 Carmichael numbers in under 5 seconds. Remove the print to actually get ideone to run to completion.

N=180000

maxprime = 18*N+1
prime = [1]*maxprime
prime[0] = prime[1] = 0
for i in xrange(maxprime):
  if prime[i]:
    for j in xrange(2*i,maxprime,i):prime[j]=0

for i in xrange(N):
  a=6*i+1
  b=12*i+1
  c=18*i+1
  if prime[a] and prime[b] and prime[c]:print a*b*c

Note that you could make this a whole lot faster in C.

\$\endgroup\$
2
  • \$\begingroup\$ You are printing Carmichael numbers only of a particular type. \$\endgroup\$
    – fR0DDY
    Feb 15, 2011 at 4:06
  • 1
    \$\begingroup\$ Indeed I am.... \$\endgroup\$ Feb 15, 2011 at 4:56
1
\$\begingroup\$

unGolf3D dear Java :) maxCount=1193221 (at my local) and 552721 (at ideone) in 4.8s

    public class Main{
    public static void main(String[] args){
          long start=System.currentTimeMillis();
          int N=1200000;
          int list[]=new int[N];
            for (int i=0;i<N;i++) {
                list[i]=(i+1);
            }
            for (int i=1;i<N;i++) {
                if (Math.sqrt(list[i])-Math.floor(Math.sqrt(list[i]))==0) {
                    for (int j=i;j<N;j+=(i+1)){
                        list[j]=-1;
                    }
                }
            }
        boolean Carmichael=false;
        int chemCount=0;
        int product=1;
          for(int i=561  ;(System.currentTimeMillis()-start)/1000.0<=5;i++){
              if(!isPrime(i)&&list[i-1]==i){
                  for(int div=2;div<Math.sqrt(i);div++){
                      if(i%div==0 && isPrime(div)){
                          if((i-1)%(div-1)==0){
                              Carmichael= true;
                              product*=div;
                              chemCount++;
                              if (product==i) break;
                              }
                          else{
                              Carmichael=false;
                              break;
                              }
                      }
                  }
                  if (Carmichael && chemCount>=3 && product==i)
                      System.out.println(i);
                  Carmichael=false;
                   chemCount=0;
                   product=1;
              }
          }
  }

        static boolean isPrime(long number){
            if (number==2)
                return true;
            if (number==1)
                return false;
            for (long i=2;i<=Math.sqrt(number);i++) {
                if (number%i==0)
                    return false;
            }
            return true;
        }
    }

My JVM gets this to 1193221 number in 5 sec... I am sure this is slower than any other scripting languages and algo!! :( ... OTOH, http://www.ideone.com/clone/oNJin fails to create this big process.. calculates upto 552721 in 4.8 sec and breaks after that throwing http://en.wikipedia.org/wiki/SIGXCPU :)

\$\endgroup\$
3
  • \$\begingroup\$ High complexity involved, Less Time, and on top Java!! \$\endgroup\$ Feb 14, 2011 at 21:47
  • \$\begingroup\$ Please mention the number of Carmichael number it generates, rather than the maximum Carmichael number you found. \$\endgroup\$
    – fR0DDY
    Feb 15, 2011 at 3:30
  • \$\begingroup\$ Oops that's only 49 :(... I thought that we need to quote the max "challenge is to print maximum of them in less than 5 secs" \$\endgroup\$ Feb 15, 2011 at 4:11
1
\$\begingroup\$

C++

Machine i5-760 @ 2.8Ghz, only one core used

Numbers found in 5 seconds: 129

#define WIN32_MEAN_AND_LEAN
#include <windows.h>
#include <iostream>
using namespace std;

bool IsPrime (int v)
{
  bool
    is_prime = true;

  for (int i = 3, i2 = 9 ; i2 < v ; ++i, i2 += i + i + 1)
  {
    if (v % i == 0)
    {
      is_prime = false;
      break;
    }
  }

  return is_prime;
}

int Prime (int i, int &i2)
{
  static int
    primes [10000] = {2, 3, 5},
    primes2 [10000] = {4, 9, 25},
    num_primes = 3;

  while (i >= num_primes)
  {
    int 
      v = primes [num_primes - 1];

    do
    {
      v += 2;
    } while (!IsPrime (v));

    primes [num_primes] = v;
    primes2 [num_primes] = v * v;
    ++num_primes;
  }

  i2 = primes2 [i];
  return primes [i];
}

bool IsCarmichael (int n) 
{
   if (n < 2) 
   {
      return false;
   }

   int
     k = n--,
     i2 = 4;

   for (int i = 1, p = 2 ; i2 <= k ; p = Prime (i++, i2))
   {
      if (k % p == 0) 
      {
         k /= p;

         if (k % p == 0 || n % (p - 1) != 0)
         {
            return false;
         }

         i = 0;
      }
   }

   return --k != n && n % k == 0;
}

int main ()
{
  DWORD
    end = GetTickCount () + 5000;

  int
    count = 0,
    n = 561;

  while (GetTickCount () < end)
  {
    for (int i = 0 ; i < 1000 ; ++i)
    {
      if (IsCarmichael (n))
      {
        ++count;
        cout << n << endl;
      }
      ++n;
    }
  }

  cout << count << " Carmichael numbers found out of " << (n - 561) << " values tested." << endl;
}
\$\endgroup\$
1
  • \$\begingroup\$ Can you fix it to run on ideone.com? e.g replace DWORD with unsigned. \$\endgroup\$
    – Alexandru
    Feb 25, 2011 at 22:12
0
\$\begingroup\$

D Language, 139 numbers, 18900973 max

Implemented Korselt's criterion to extract Carmichael numbers, considering only odd numbers greater than 3 since an even number cannot satisfy the relation.

Code

import std.stdio, std.math;

bool isCarmichael(uint n) pure nothrow @nogc {
   auto maxValue = ceil(sqrt(float(n)));
   uint cn = n;

   for (uint i = 3; i < maxValue; i += 2) {
      if (cn % i == 0) {
         cn /= i;

         if (cn % i == 0)
            return false;

         if ((n-1) % (i-1))
            return false;
      }
   }

   if (cn > 1)
      return false;

   return true;
}

void main() {
   for (uint n = 3; ; n += 2) {
      if (isCarmichael(n))
         writeln(n);
   }
}

Output

561
1105
1729
2465
2821
6601
8911
10585
15841
29341
41041
46657
52633
62745
63973
75361
101101
115921
126217
162401
172081
188461
252601
278545
294409
314821
334153
340561
399001
410041
449065
488881
512461
530881
552721
656601
658801
670033
748657
825265
838201
852841
997633
1024651
1033669
1050985
1082809
1152271
1193221
1461241
1569457
1615681
1773289
1857241
1909001
2100901
2113921
2433601
2455921
2508013
2531845
2628073
2704801
3057601
3146221
3224065
3581761
3664585
3828001
4335241
4463641
4767841
4903921
4909177
5031181
5049001
5148001
5310721
5444489
5481451
5632705
5968873
6049681
6054985
6189121
6313681
6733693
6840001
6868261
7207201
7519441
7995169
8134561
8341201
8355841
8719309
8719921
8830801
8927101
9439201
9494101
9582145
9585541
9613297
9890881
10024561
10267951
10402561
10606681
10837321
10877581
11119105
11205601
11921001
11972017
12261061
12262321
12490201
12945745
13187665
13696033
13992265
14469841
14676481
14913991
15247621
15403285
15829633
15888313
16046641
16778881
17098369
17236801
17316001
17586361
17812081
18162001
18307381
18900973

The code runs on a single core, without parallelism. No luck making it work on IDEONE though; time is checked on my machine running the process with timeout 5s UNIX command.

\$\endgroup\$
0
\$\begingroup\$

Rust

Naive brute force. It can be optimized more.

Run it on Rust Playground!

extern crate rand;

fn gcd(a: i128, b: i128) -> i128 {
    if b == 0 {
        a
    } else {
        gcd(b, a % b)
    }
}

fn power(mut a: i128, mut b: i128, m: i128) -> i128 {
    let mut result = 1;
    a %= m;
    while b > 0 {
        if b % 2 != 0 {
            result = (result * a) % m;
        }
        b >>= 1;
        a = (a * a) % m;
    }
    result
}

fn is_prime_miller_rabin(n: i128) -> bool {
    if n < 2 {
        return false;
    }
    if n != 2 && n % 2 == 0 {
        return false;
    }
    let s = ((n - 1) as u128).trailing_zeros() as i128;
    let d = (n - 1) / 2i128.pow(s as u32);

    let witnesses: [i128; 12] = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37];
    'next_witness: for &a in witnesses.iter().take_while(|&&x| x <= n - 1) {
        let mut x = power(a, d, n);
        if x == 1 || x == n - 1 {
            continue;
        }
        for _r in 1..s {
            x = power(x, 2, n);
            if x == n - 1 {
                continue 'next_witness;
            }
        }
        return false;
    }
    true
}

fn fermat_test(n: i128) -> bool {
    for a in 2..n {
        if gcd(a, n) == 1 {
            if power(a, n - 1, n) != 1 {
                return false;
            }
        }
    }
    true
}

fn is_carmichael(n: i128) -> bool {
    !is_prime_miller_rabin(n) && fermat_test(n)
}

fn main() {
    let start = 561;
    let end = 100000;  // Adjust end limit based on time constraint
    for n in start..end {
        if is_carmichael(n) {
            println!("{}", n);
        }
    }
}
\$\endgroup\$
0
\$\begingroup\$

Python 3, 1184 numbers

Somewhat simplified implementation of Richard Pinch's methods (Carmichael numbers up to 10^21, Carmichael numbers up to 10^16).

The idea is to enumerate all the Carmichael numbers less than \$N\$. We can prove that none of them will have prime factors \$>\sqrt{N}\$, so we split into two cases based on the size of the largest prime factor (specifically if it is larger or smaller than \$N^{0.3}\$).

We'll let \$n\$ be the number, and \$q\$ it's largest factor. If there is a factor greater than \$N^{0.3}\$, then since \$\frac{n}{q}\equiv1\mod{q-1}\$, enumerating all possible values of \$\frac{n}{q}\$ given \$q\$ is \$\mathcal{O}\left(N^{0.4}\right)\$. So enumerate \$q\$ from \$N^{0.3}\$ to \$N^{0.5}\$, and test with a Fermat test to base 2 and Korselt's criterion.

In the other case, we use a depth first search on products of prime factors, where we can limit the branching factor to \$N^{0.3}\$.

Try it on ideone!

import math
from collections import deque
import itertools
# from gmpy2 import powmod as pow

def main():
    x = find_carmichael(5_000_000_000)
    for idx, i in enumerate(sorted(x), 1):# + len(s)):
        print(idx, i)

def find_carmichael(X):
    if X = q: return False
        x = i
        f.append(i)
    if len(f) >= 2:
        n1 = P*q-1
        return all(n1%(p-1)==0 for p in f)

def dfs(X, primes, n=1, t=(), i=0, lcm=1):
    if math.gcd(n, lcm) > 1: return
    if lcm * n * len(primes) > X:
        for Q in range(pow(n, -1, lcm), X//n + 1, lcm):
            if math.gcd(n, Q) > 1: continue
            N = n*Q
            if pow(2, N, Q) == 2:
                if is_carmichael0(n, Q, t):
                    yield N
        return
    for i1 in range(i, len(primes)):
        p = primes[i1]
        n1 = n*p
        if n1 > X: break
        if len(t) >= 2 and (n-1)%(p-1) == 0 and n1%lcm == 1:
            yield n1
        yield from dfs(X, primes, n1, t+(p,), i1+1, math.lcm(lcm, p-1))

def is_carmichael0(n, Q, t):
    x = 0
    f = list(t)
    for i in factor(Q):
        if i == x: return False
        if i in t: return False
        x = i
        f.append(i)
    if len(f) >= 3:
        n1 = n*Q-1
        return all(n1%(p-1)==0 for p in f)

class LazyList:
    def __init__(self, I):
        self.I = I
        self.l = []
    def __getitem__(self, i):
        try:
            return self.l[i]
        except IndexError:
            self.l.extend(itertools.islice(self.I, i+1 - len(self.l)))
            return self.l[i]
    def range(self, a, b=None):
        if b is None:
            a, b = 0, a
        for i in self:
            if i = b: break
            yield i

def primes_():
    yield 2
    l = [3,5,7,11,13,17,19,23]
    yield from l
    p2s = deque([(p, p*p) for p in l[2:]])
    n = 23
    stuff = [(0, 3), (4, 5)]
    while True:
        p2, next_p2 = p2s.popleft()
        R = range(n+4, next_p2, 2)
        z = [[] for _ in R]
        for i, p in stuff:
            z[i].append(p)
        stuff.clear()
        stuff.append((p2-1, p2))
        for i, (l, n) in enumerate(zip(z, R)):
            # print(i, l, n)
            if l:
                for p in l:
                    i1 = i + p
                    if i1 >= len(R):
                        stuff.append((i1 - len(R) - 1, p))
                    else:
                        z[i1].append(p)
            else:
                yield n
                p2s.append((n, n*n))

primes = LazyList(primes_())

def factor(n):
    for p in primes:
        while n%p == 0:
            yield p
            n //= p
        if p*p > n:
            if n > 1: yield n
            return

if __name__ == '__main__': main()
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.