## Lucas-Carmichael Numbers

### January 6, 2015

We iterate through the odd numbers from 3 to a million, factoring each, check that the number is composite and square-free, then check that all its factors meet the divisibility requirement:

```(define (lucas-carmichael limit)   (list-of n   (n range 3 limit 2)   (fs is (factors n))   (< 1 (length fs)) ; composite   (no-dups? fs) ; square-free   (all? (lambda (f) (divides? (+ f 1) (+ n 1))) fs)))```

We use auxiliary functions to check for duplicates and divisibility:

```(define (no-dups? xs)   (or (null? xs)       (and (not (member (car xs) (cdr xs)))            (no-dups? (cdr xs)))))```

`(define (divides? d n) (zero? (modulo n d)))`

For numbers as small as a million, the easiest factoring algorithm is a simple wheel:

```(define (factors n) ; 2,3,5-wheel   (let ((wheel (vector 1 2 2 4 2 4 2 4 6 2 6)) (next 0))     (let loop ((n n) (f 2) (fs (list)))       (if (= n 1) (reverse fs)         (if (< n (* f f)) (reverse (cons n fs))           (if (zero? (modulo n f))               (loop (/ n f) f (cons f fs))               (let ((f (+ f (vector-ref wheel next))))                 (set! next (if (= next 10) 3 (+ next 1)))                 (loop n f fs))))))))```

And that’s it. Here’s the output:

```> (lucas-carmichael 1000000) (399 935 2015 2915 4991 5719 7055 8855 12719 18095 20705  20999 22847 29315 31535 46079 51359 60059 63503 67199 73535  76751 80189 81719 88559 90287 104663 117215 120581 147455  152279 155819 162687 191807 194327 196559 214199 218735  230159 265895 357599 388079 390335 482143 588455 653939  663679 676799 709019 741311 760655 761039 776567 798215  880319 895679 913031 966239 966779 973559)```

It’s nine hundred years to the next Lucas-Carmichael year, so enjoy this one while you can. We used list comprehensions and the `all?` function from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/AXSKaB7M.

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