Rosetta Code

Lucas-Carmichael numbers

A Lucas-Carmichael number, is a squarefree composite integer n such that p+1 divides n+1 for all the prime factors p of n.

Lucas-Carmichael numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Task

Write a function that generates all the Lucas-Carmichael numbers with n prime factors in a given range [A, B].

  • Show the smallest Lucas-Carmichael with n prime factors, for n = 3..12.
  • Show the count of Lucas-Carmichael numbers less than 10^n, for n = 1..10.


Algorithm
Algorithm for generating Lucas-Carmichael numbers with n prime factors
Algorithm for generating Lucas-Carmichael numbers with n prime factors


See also


Julia 

using Primes

const BIG = false       # true to use big integers

function big_prod(arr)
    BIG || return prod(arr)
    r = big(1)
    for n in arr
        r *= n
    end
    return r
end

function lucas_carmichael(A, B, k)

    LC = []
    max_p = isqrt(B+1)-1
    A = max(A, fld(big_prod(primes(prime(k+1))), 2))

    F = function(m, L, lo, k)

        hi = round(Int64, fld(B, m)^(1/k))

        if (lo > hi)
            return nothing
        end

        if (k == 1)

            hi = min(hi, max_p)
            lo = round(Int64, max(lo, cld(A, m)))
            lo > hi && return nothing

            t = L - invmod(m, L)
            t > hi && return nothing

            while (t < lo)
                t += L
            end

            for p in t:L:hi
                if (isprime(p))
                    n = m*p
                    if ((n+1) % (p+1) == 0)
                        push!(LC, n)
                    end
                end
            end

            return nothing
        end

        for p in primes(lo, hi)
            if (gcd(m, p+1) == 1)
                F(m*p, lcm(L, p+1), p+1, k-1)
            end
        end
    end

    F((BIG ? big(1) : 1), (BIG ? big(1) : 1), 3, k)

    return sort(LC)
end

function LC_with_n_primes(n)

    n < 3 && return nothing

    x = big_prod(primes(prime(n + 1))) >> 1
    y = 2 * x

    while true
        LC = lucas_carmichael(x, y, n)
        if (length(LC) >= 1)
            return LC[1]
        end
        x = y + 1
        y = 2 * x
    end
end

println("Least Lucas-Carmichael number with n prime factors:")

for n in 3:12
    println([n, LC_with_n_primes(n)])
end

function LC_count(A, B)
    count = 0
    for k in 3:10^6
        if (big_prod(primes(prime(k+1)))/2 > B)
            break
        end
        count += length(lucas_carmichael(A, B, k))
    end
    return count
end

println("\nNumber of Lucas-Carmichael numbers less than 10^n:")

for n in 1:10
    println([n, LC_count(1, 10^n)])
end
Output:
Least Lucas-Carmichael number with n prime factors:
[3, 399]
[4, 8855]
[5, 588455]
[6, 139501439]
[7, 3512071871]
[8, 199195047359]
[9, 14563696180319]
[10, 989565001538399]
[11, 20576473996736735]
[12, 4049149795181043839]

Number of Lucas-Carmichael numbers less than 10^n:
[1, 0]
[2, 0]
[3, 2]
[4, 8]
[5, 26]
[6, 60]
[7, 135]
[8, 323]
[9, 791]
[10, 1840]

PARI/GP 

lucas_carmichael(A, B, k) = {
    A=max(A, vecprod(primes(k+1))\2);
    my(max_p=sqrtint(B+1)-1);
    (f(m, l, lo, k) =
        my(list=List());
        my(hi=sqrtnint(B\m, k));
        if(lo > hi, return(list));
        if(k==1,
            hi = min(hi, max_p);
            lo=max(lo, ceil(A/m));
            my(t=lift(-1/Mod(m,l)));
            while(t < lo, t += l);
            forstep(p=t, hi, l,
                if(isprime(p),
                    my(n=m*p);
                    if((n+1)%(p+1) == 0, listput(list, n));
                );
            ),
        \\ else
            forprime(p=lo, hi,
                if(gcd(m, p+1) == 1,
                    list=concat(list, f(m*p, lcm(l, p+1), p+1, k-1))
                );
            );
        );
        list;
    );
    vecsort(Vec(f(1, 1, 3, k)));
};

LC_count(n) = {
    my(count=0);
    for(k=3, oo,
        if(vecprod(primes(k+1))\2 > n, break);
        count += #lucas_carmichael(1, n, k)
    );
    count;
};

LC_with_n_primes(n) = {
    if(n < 3, return());
    my(x=vecprod(primes(n+1))\2, y=2*x);
    while(1,
        my(v=lucas_carmichael(x,y,n));
        if(#v >= 1, return(v[1]));
        x=y+1; y=2*x;
    );
};

print("Least Lucas-Carmichael number with n prime factors:")
for(n=3, 12, print([n, LC_with_n_primes(n)]));

print("\nNumber of Lucas-Carmichael numbers less than 10^n:")
for(n=1, 10, print([n, LC_count(10^n)]));
Output:
Least Lucas-Carmichael number with n prime factors:
[3, 399]
[4, 8855]
[5, 588455]
[6, 139501439]
[7, 3512071871]
[8, 199195047359]
[9, 14563696180319]
[10, 989565001538399]
[11, 20576473996736735]
[12, 4049149795181043839]

Number of Lucas-Carmichael numbers less than 10^n:
[1, 0]
[2, 0]
[3, 2]
[4, 8]
[5, 26]
[6, 60]
[7, 135]
[8, 323]
[9, 791]
[10, 1840]

Perl 

use 5.020;
use ntheory      qw(:all);
use experimental qw(signatures);

sub divceil ($x, $y) {    # ceil(x/y)
    (($x % $y == 0) ? 0 : 1) + int($x / $y);
}

sub LC_in_range ($A, $B, $k) {

    my @LC;
    my $max_p = sqrtint($B + 1) - 1;

    $A = vecmax(pn_primorial($k + 1) >> 1, $A);

    sub ($m, $L, $lo, $k) {

        my $hi = rootint(int($B / $m), $k);

        return if ($lo > $hi);

        if ($k == 1) {

            $hi = $max_p if ($hi > $max_p);
            $lo = vecmax($lo, divceil($A, $m));
            $lo > $hi && return;

            my $t = $L - invmod($m, $L);
            $t += $L while ($t < $lo);

            for (my $p = $t ; $p <= $hi ; $p += $L) {
                if (is_prime($p)) {
                    my $n = $m * $p;
                    if (($n + 1) % ($p + 1) == 0) {
                        push @LC, $n;
                    }
                }
            }

            return;
        }

        foreach my $p (@{primes($lo, $hi)}) {
            if (gcd($m, $p + 1) == 1) {
                __SUB__->($m * $p, lcm($L, $p + 1), $p + 1, $k - 1);
            }
        }
      }
      ->(1, 1, 3, $k);

    return sort { $a <=> $b } @LC;
}

sub LC_with_n_primes ($n) {
    return if ($n < 3);

    my $x = pn_primorial($n + 1) >> 1;
    my $y = 2 * $x;

    for (; ;) {
        my @LC = LC_in_range($x, $y, $n);
        @LC and return $LC[0];
        $x = $y + 1;
        $y = 2 * $x;
    }
}

sub LC_count ($A, $B) {
    my $count = 0;
    for (my $k = 3 ; ; ++$k) {
        if (pn_primorial($k + 1) / 2 > $B) {
            last;
        }
        my @LC = LC_in_range($A, $B, $k);
        $count += @LC;
    }
    return $count;
}

say "Least Lucas-Carmichael number with n prime factors:";

foreach my $n (3 .. 12) {
    printf("%2d: %s\n", $n, LC_with_n_primes($n));
}

say "\nNumber of Lucas-Carmichael numbers less than 10^n:";

my $acc = 0;
foreach my $n (1 .. 10) {
    my $c = LC_count(10**($n - 1), 10**$n);
    printf("%2d: %s\n", $n, $acc + $c);
    $acc += $c;
}
Output:
Least Lucas-Carmichael number with n prime factors:
 3: 399
 4: 8855
 5: 588455
 6: 139501439
 7: 3512071871
 8: 199195047359
 9: 14563696180319
10: 989565001538399
11: 20576473996736735
12: 4049149795181043839

Number of Lucas-Carmichael numbers less than 10^n:
 1: 0
 2: 0
 3: 2
 4: 8
 5: 26
 6: 60
 7: 135
 8: 323
 9: 791
10: 1840

Phix

Translation of: Julia

mul_inv() copied from Modular_inverse#Phix 

requires(64,true)

function mul_inv(integer a, n)
    if n<0 then n = -n end if
    if a<0 then a = n - mod(-a,n) end if
    integer t = 0,  nt = 1,
            r = n,  nr = a;   
    while nr!=0 do
        integer q = floor(r/nr)
        {t, nt} = {nt, t-q*nt}
        {r, nr} = {nr, r-q*nr}
    end while
    if r>1 then return "a is not invertible" end if
    if t<0 then t += n end if
    return t
end function
 
function F(sequence LC, atom A, B, integer max_p, m, L, lo, k)

    integer hi = round(power(floor(B/m),1/k))

    if lo<=hi then
        if k==1 then
            hi = min(hi, max_p)
            lo = round(max(lo, ceil(A/m)))
            if lo <= hi then
                integer imL = mul_inv(m, L),
                          t = L - imL
                if t <= hi then
                    while t<lo do
                        t += L
                    end while
                    for p=t to hi by L do
                        if is_prime(p) then
                            atom n = m*p
                            if remainder(n+1,p+1)==0 then
                                LC &= n
                            end if
                        end if
                    end for
                end if
            end if
        else
            sequence plohi = get_primes_le(hi)
            lo = abs(binary_search(lo,plohi))
            plohi = plohi[lo..$]
            for p in plohi do
                if gcd(m, p+1) == 1 then
                    LC = F(LC, A, B, max_p, m*p, lcm(L, p+1), p+1, k-1)
                end if
            end for
        end if
    end if
    return LC
end function

function lucas_carmichael(atom A, B, integer k)
    atom max_p = floor(sqrt(B+1))-1,
          fppp = floor(product(get_primes(-k-1))/2)
    return sort(F({}, max(A,fppp), B, max_p, 1, 1, 3, k))
end function

function LC_with_n_primes(integer n)
    atom x = product(get_primes(-n-1))/2,
         y = 2 * x
    while true do
        sequence LC = lucas_carmichael(x, y, n)
        if length(LC)>=1 then
            return LC[1]
        end if
        x = y + 1
        y = 2 * x
    end while
end function

printf(1,"Least Lucas-Carmichael number with n prime factors:\n")

for n=3 to 12 do
    printf(1,"%d: %,d\n",{n, LC_with_n_primes(n)})
end for

function LC_count(atom A, B)
    integer count = 0
    for k=3 to 1e6 do
        if product(get_primes(-k-1))/2 > B then
            exit
        end if
        count += length(lucas_carmichael(A, B, k))
    end for
    return count
end function

printf(1,"\nNumber of Lucas-Carmichael numbers less than 10^n:\n")

for n=1 to 10 do
    printf(1,"%d: %,d\n",{n, LC_count(1, power(10,n))})
end for
Output:
Least Lucas-Carmichael number with n prime factors:
3: 399
4: 8,855
5: 588,455
6: 139,501,439
7: 3,512,071,871
8: 199,195,047,359
9: 14,563,696,180,319
10: 989,565,001,538,399
11: 20,576,473,996,736,735
12: 4,049,149,795,181,043,839

Number of Lucas-Carmichael numbers less than 10^n:
1: 0
2: 0
3: 2
4: 8
5: 26
6: 60
7: 135
8: 323
9: 791
10: 1,840

Python

Uses the SymPy library. 

from sympy.ntheory import sieve, isprime, prime
from sympy.core import integer_nthroot
from math import lcm, gcd, isqrt

def lucas_carmichael(A, B, n):
    max_p = isqrt(B)+1

    def f(m, l, lo, k):

        if k == 1:

            lo = max(lo, A // m + (1 if A % m else 0))
            hi = min(B // m + 1, max_p)

            u = l - pow(m, -1, l)
            while u < lo: u += l
            if u > hi: return

            for p in range(u, hi, l):
                if (m*p+1) % (p+1) == 0 and isprime(p):
                    yield m*p

        else:
            hi = (integer_nthroot(B // m, k))[0]+1
            for p in sieve.primerange(lo, hi):
                if gcd(m, p+1) == 1:
                    yield from f(m*p, lcm(l, p+1), p + 2, k - 1)

    return sorted(f(1, 1, 3, n))

def LC_with_n_primes(n):
    x = 2
    y = 2*x
    while True:
        LC = lucas_carmichael(x, y, n)
        if len(LC) >= 1: return LC[0]
        x = y+1
        y = 2*x

def LC_count(A, B):
    k = 3
    l = 3*5*7
    count = 0
    while l < B:
        count += len(lucas_carmichael(A, B, k))
        k += 1
        l *= prime(k+1)
    return count

print("Least Lucas-Carmichael number with n prime factors:")

for n in range(3, 12+1):
    print("%2d: %d" % (n, LC_with_n_primes(n)))

print("\nNumber of Lucas-Carmichael numbers less than 10^n:")

for n in range(1, 10+1):
    print("%2d: %d" % (n, LC_count(1, 10**n)))
Output:
Least Lucas-Carmichael number with n prime factors:
 3: 399
 4: 8855
 5: 588455
 6: 139501439
 7: 3512071871
 8: 199195047359
 9: 14563696180319
10: 989565001538399
11: 20576473996736735
12: 4049149795181043839

Number of Lucas-Carmichael numbers less than 10^n:
 1: 0
 2: 0
 3: 2
 4: 8
 5: 26
 6: 60
 7: 135
 8: 323
 9: 791
10: 1840

Raku

Translation of: Perl
# 20231224 Raku programming solution

sub primorial($n) { [*] ((2..*).grep: *.is-prime)[^$n] }
sub divceil($x, $y) { ($x %% $y ?? 0 !! 1) + $x div $y } # ceil(x/y)

sub LC_in_range ($A is copy, $B, $k) {

   my ($max_p, @LC) = Int(sqrt($B + 1) - 1);
   $A = max(primorial($k + 1) +> 1, $A);

   sub SUB($m, $L, $lo is copy, $k) {

      if $k == 1 {

         my $hi = min($B div $m, $max_p);
         $lo = max($lo, divceil($A, $m));

         my $t = $L - expmod($m, -1, $L);
         $t += $L while $t < $lo;
         return if $t > $hi;

         for $t, $t+$L ... $hi -> $p {
            if $p.is-prime {
               my $n = $m * $p;
               @LC.push($n) if ($n + 1) %% ($p + 1);
            }
         }
         return;
      }

      for $lo .. Int(($B div $m)**(1/$k)) -> $p {
         if $p.is-prime and ($m gcd ($p + 1)) == 1 {
            SUB($m * $p, ($L lcm ($p + 1)), $p + 1, $k - 1)
         }
      }
   };

   SUB(1, 1, 3, $k);

   return @LC.sort;
}

sub LC_with_n_primes ($n) {
   return if $n < 3;

   my $y = 2 * ( my $x = primorial($n + 1) +> 1);

   loop {
      my @LC = LC_in_range($x, $y, $n);
      return @LC[0] if @LC.Bool;
      $x = $y + 1;
      $y = 2 * $x;
   }
}

sub LC_count ($A, $B) {
   my $count = 0;
   for 3 .. Inf -> $k {
      last if primorial($k + 1) / 2 > $B;
      $count += LC_in_range($A, $B, $k).elems
   }
   return $count;
}

say "Least Lucas-Carmichael number with n prime factors:";

for 3 .. 12 -> $n { printf("%2d: %d\n", $n, LC_with_n_primes($n)) }

say "\nNumber of Lucas-Carmichael numbers less than 10^n:";

my $acc = 0;
for 1 .. 10 -> $n {
   my $c = LC_count(10**($n - 1), 10**$n);
   printf("%2d: %s\n", $n, $acc + $c);
   $acc += $c
}

Output is the same as Perl.

REXX

Include: How to use
Include: Source code
Several procedures from the field 'number theory' are used: Prime and SquareFree from module numbers and Ufactors from module Sequences. 

-- 8 May 2025
include Settings
numeric digits 10

say 'LUCAS-CARMICHAEL NUMBERS'
say version
say
call GetNumbers 1,1e5
call ShowNumbers 1,1e5
call Timer
exit

GetNumbers:
procedure expose luca. ufac.
arg xx,yy
say 'Get LuCa numbers...'
xx = xx+Even(xx)
luca. = 0; n = 0
do i = xx by 2 to yy
   if Prime(i) then
      iterate i
   if \ Squarefree(i) then
      iterate i
   i1 = i+1
   a = Ufactors(i)
   do j = 1 to a
      u1 = ufac.j+1
      if i1//u1 > 0 then
         iterate i
   end
   n = n+1; luca.n = i
end
luca.0 = n
say n 'numbers found'
say
return

ShowNumbers:
procedure expose luca.
arg xx,yy
say 'Lucas-Carmichael numbers between' xx 'and' yy'...'
do i = 1 to luca.0
   call charout ,right(luca.i,7)
   if i//10 = 0 then
      say
end
say; say
return

include Numbers
include Sequences
include Functions
include Helper
include Abend
Output:
LUCAS-CARMICHAEL NUMBERS
REXX-Regina_3.9.7(MT) 5.00 18 Mar 2025

Get LuCa numbers...
27 numbers found

Lucas-Carmichael numbers between 1 and 1E5...
      1    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

3.802 seconds

Sidef

The function is also built-in as n.lucas_carmichael(A,B). 

func LC_in_range(A, B, k) {

    var LC = []

    func (m, L, lo, k) {

        var hi = idiv(B,m).iroot(k)

        if (lo > hi) {
            return nil
        }

        if (k == 1) {

            hi = min(B.isqrt, hi)
            lo = max(lo, idiv_ceil(A, m))
            lo > hi && return nil

            var t = mulmod(m.invmod(L), -1, L)

            t > hi && return nil
            t += L while (t < lo)

            for p in (range(t, hi, L).primes) {
                with (m*p) {|n|
                    LC << n if (p+1 `divides` n+1)
                }
            }

            return nil
        }

        each_prime(lo, hi, {|p|
            __FUNC__(m*p, lcm(L, p+1), p+1, k-1) if m.is_coprime(p+1)
        })
    }(1, 1, 3, k)

    return LC.sort
}

func LC_with_n_primes(n) {
    return nil if (n < 3)

    var x = pn_primorial(n+1)/2
    var y = 2*x

    loop {
        var arr = LC_in_range(x,y,n)
        arr && return arr[0]
        x = y+1
        y = 2*x
    }
}

func LC_count(A, B) {
    var count = 0
    for k in (3..Inf) {
        if (pn_primorial(k+1)/2 > B) {
            break
        }
        count += LC_in_range(A,B,k).len
    }
    return count
}

say "Least Lucas-Carmichael number with n prime factors:"

for n in (3..12) {
    say "#{'%2d'%n}: #{LC_with_n_primes(n)}"
}

say "\nNumber of Lucas-Carmichael numbers less than 10^n:"

var acc = 0
for n in (1..10) {
    var c = LC_count(10**(n-1), 10**n)
    say "#{'%2d'%n}: #{c + acc}"
    acc += c
}
Output:
Least Lucas-Carmichael number with n prime factors:
 3: 399
 4: 8855
 5: 588455
 6: 139501439
 7: 3512071871
 8: 199195047359
 9: 14563696180319
10: 989565001538399
11: 20576473996736735
12: 4049149795181043839

Number of Lucas-Carmichael numbers less than 10^n:
 1: 0
 2: 0
 3: 2
 4: 8
 5: 26
 6: 60
 7: 135
 8: 323
 9: 791
10: 1840

Wren

Translation of: Julia
Library: Wren-psieve
Library: Wren-math
Library: Wren-big
Library: Wren-fmt
import "./psieve" for Primes
import "./math" for Int, Nums
import "./big" for BigInt
import "./fmt" for Fmt

var lucasCarmichael = Fn.new { |A, B, k|
    var LC = []
    var maxP = (B+1).sqrt.floor - 1
    var p = Primes.nthAfter(k+1, 1)
    var primes = Primes.between(2, p)
    var x = (Nums.prod(primes)/2).floor
    A = A.max(x)
    var M = 1
    var P = 1
    var F
    F = Fn.new { |m, L, lo, k|
        var hi = (B/m).floor.pow(1/k).round
        if (lo > hi) return []
        if (k == 1) {
            hi = hi.min(maxP)
            lo = lo.max((A/m).ceil).round
            if (lo > hi) return []
            var t
            if (m <= Num.maxSafeInteger) {
                t = L - Int.modInv(m, L)
            } else {
                var bm = BigInt.new(M) * P
                t = L - bm.modInv(L).toNum
            }
            if (t > hi) return []
            while (t < lo) t = t + L
            p = t
            while (p <= hi) {
                if (Int.isPrime(p)) {
                    var n = m * p
                    if (n + 1 <= Num.maxSafeInteger) {
                        if ((n+1) % (p+1) == 0) LC.add(n)
                    } else {
                        var bn = BigInt.new(M) * P * p
                        if ((bn+1) % (p+1) == 0) LC.add(bn)
                    }
                }
                p = p + L
            }
            return []
        }
        for (p in Primes.between(lo, hi)) {
            if (Int.gcd(m, p+1) == 1) {
                M = m
                P = p
                F.call(m*p, Int.lcm(L, p+1), p+1, k-1)
            }
        }
    }
    F.call(1, 1, 3, k)
    if (LC.any { |e| e is BigInt }) {
        for (i in 0...LC.count) {
            if (LC[i] is Num) LC[i] = BigInt.new(LC[i])
        }
    }
    return LC.sort()
}

var lcWithNPrimes = Fn.new { |n|
    if (n < 3) return []
    var p = Primes.nthAfter(n+1, 1)
    var primes = Primes.between(2, p)
    var x = (Nums.prod(primes)/2).floor
    var y = 2 * x
    while (true) {
        var LC = lucasCarmichael.call(x, y, n)
        if (LC.count >= 1) return LC[0]
        x = y + 1
        y = 2 * x
    }
}

var lcCount = Fn.new { |A, B|
    var count = 0
    for (k in 3..1e6) {
        var p = Primes.nthAfter(k+1, 1)
        var primes = Primes.between(2, p)
        var x = (Nums.prod(primes)/2).floor
        if (x > B) break
        count = count + lucasCarmichael.call(A, B, k).count
    }
    return count
}

System.print("Least Lucas-Carmichael number with n prime factors:")   
for (n in 3..12) {
    Fmt.print("$2d: $i", n, lcWithNPrimes.call(n))
}

System.print("\nNumber of Lucas-Carmichael numbers less than 10^n:")
for (n in 1..10) {
    Fmt.print("$2d: $d", n, lcCount.call(1, 10.pow(n)))
}
Output:
Least Lucas-Carmichael number with n prime factors:
 3: 399
 4: 8855
 5: 588455
 6: 139501439
 7: 3512071871
 8: 199195047359
 9: 14563696180319
10: 989565001538399
11: 20576473996736735
12: 4049149795181043839

Number of Lucas-Carmichael numbers less than 10^n:
 1: 0
 2: 0
 3: 2
 4: 8
 5: 26
 6: 60
 7: 135
 8: 323
 9: 791
10: 1840
Retrieved from "https://rosettacode.org/wiki/Lucas-Carmichael_numbers?oldid=380770"