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