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] = 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]]
    
  4. Kaushik said

    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)

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: