## 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

sieve = *(limit//2)
sieve = 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

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) 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);
}
```