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.
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;
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.
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)