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.
Programming Praxis 
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;
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, factorsI 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]]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)
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); }