Lucas-Carmichael numbers

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.

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.

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

• OEIS A006972: Lucas-Carmichael numbers: squarefree composite numbers k such that p | k => p+1 | k+1

• OEIS A216928: Least Lucas-Carmichael number with n prime factors

• OEIS A216929: Number of Lucas-Carmichael numbers less than 10^n

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.

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