Lucas-Carmichael Numbers

January 6, 2015

We start the new year with a simple task from number theory. A Lucas-Carmichael number is a positive composite integer n, odd and square-free, such that p + 1 divides n + 1 for all prime factors p of n. For instance, 2015 is a Lucas-Carmichael number because 2015 = 5 × 13 × 31 and 5 + 1 = 6, 13 + 1 = 14, and 31 + 2 = 32 all divide 2015 + 1 = 2016. The restriction that n be square-free prevents the cubes of prime numbers, such as 8 or 27, from being considered as Lucas-Carmichael numbers.

Your task is to write a program to find all the Lucas-Carmichael numbers less than a million. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

5 Responses to “Lucas-Carmichael Numbers”

  1. danaj said

    Perl, 5s for 1e7 (135 L-C numbers), <1.5min for 1e8 (323 L-C numbers). Perl:

    use warnings;
    use strict;
    use ntheory qw/:all/;
    foroddcomposites {
    my $n=$_;
    say unless moebius($n) == 0 || vecsum( map { ($n+1) % ($_+1) != 0 } factor($n) );
    } 1e6;

  2. Rutger said
    def factorize(number):
        # returns a list of all prime factors of number
        factors = []
        while not (number % 2):
            number = number / 2
            factors.append(2)
        i = 3
        while i <= number**0.5 and number-1:
            if not (number % i):
                factors.append(i)
                number = number / i
            else:
                i += 2
        if number != 1 and i >= number**0.5:
            factors.append(number)
        return factors
    
    def lucas_carmichael(number, factors):
        if len(factors) == 1:
            # prime
            return False
        if any((factors.count(f) > 1 for f in factors)):
            # not square free
            return False
        for f in factors:
            if (number+1) % (f+1):
                return False
        return True
    
    for n in range(3,1000000,2):
        factors = factorize(n)
        if lucas_carmichael(n, factors):
            print n, factors
    
  3. Mike said

    I decided to use a sieve.

    For p = 3, write out the multiples of p and of (p+1):

    3 6 9 12 15 18 21 24 27
    4 8 12 16 20 24 28

    The only values of n for which (n+1)%(p+1) == 0 are n = 15, 27, 39, …. Similarly, for p = 5, the only values are n = 35, 65, 95, …. In general, for a prime ‘p’, the only values of ‘n’ for which (n+1)%(p+1) == 0 are given by: p + k*p*(p+1) for k = 1,2,3,…

    Initialize the sieve to all 1s.

    For each odd prime ‘p’:
    for k in 1, 2, 3, …:
    sieve[p+k*p*(p+1)] *= p

    return a list of all n where sieve[n] == n

    The code below sieves the primes I parallel with sieving the LCNs.
    The sieves also only contain odd n, so the index for a value of n is given by n//2.

    Finds 323 LCNs under 1e8 in about 12 seconds.

    import math
    
    def lcn(limit=1000000):
        oddprimes = list(range(1,math.ceil(math.sqrt(limit)),2))
        oddprimes[0] = 0
    
        sieve = [1]*(limit//2)
        sieve[0] = 0
    
        for p in oddprimes:
            if p:
                # sieve the primes
                for j in range(p*p//2, len(oddprimes), p):
                    oddprimes[j] = 0
        
                # accumulate factors for possible LCNs
                start = p // 2
                step = p * (p + 1) // 2
                for k in range(start+step, len(sieve), step):
                    sieve[k] *= p
    
        # n is an LCN if product of prime factors == n
        return [p for p in range(1, limit, 2) if p==sieve[p//2]]
    
  4. Kaushik said

    Haskell

    import Data.List

    factors :: Int -> [Int]
    factors n = factors’ n 2
    where
    factors’ n f
    | f*f > n = [n]
    | n `mod` f == 0 = f : factors’ (n `div` f) f
    | otherwise = factors’ n (f+1)

    noDupes :: [Int] -> Bool
    noDupes xs = length(nub xs) == length xs

    lcn :: [Int] -> [Int]
    lcn xs = filter (\x ->
    let fs = factors x
    in
    length(fs) > 1 && noDupes(fs) && (all (\f -> (x + 1) `mod` (f + 1) == 0) fs) ) xs

    limit = 1000000
    n = [3,5..limit]

    main = print (lcn n)

  5. Dan said
    function getPrimeFactors(n){
    	j = [];
    	k = n;
    	for(i=2;i<=k;i++){
    		while(k%i==0){
    			j.push(i);
    			k /= i;
    		}
    		if(Math.sqrt(k)<i && k>1){
    			j.push(k);
    			break;
    		}
    	}
    	return j;
    }
    function isLucasCarmichael(n){
    	q = getPrimeFactors(n);
    	if(q.length==1) return false;
    	if(q[q.length-1]==q[0]) return false;
    	for(p=0;p<q.length;p++){
    		if((n+1)%(q[p]+1)!=0) return false;
    		if(i!=0) if(q[p]==q[p-1]) return false;
    	}
    	return true;
    }
    function run(){
    	for(r=2;r<1000000;r++) if(isLucasCarmichael(r)) console.log(r);
    }
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: