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

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)