experimental-projects/oeis-research/Tim Johannes Ohrtmann/Number of Lucas-Carmichael numbers less than 10^n/prog.sf

#!/usr/bin/ruby

# Number of Lucas-Carmichael numbers less than 10^n.
# https://oeis.org/A216929

# Known terms:
#   0, 0, 2, 8, 26, 60, 135, 323, 791, 1840, 4216, 9967

# New terms (12 November 2023):
#   a(13) = 23070  (took 5 minutes)
#   a(14) = 54356  (took ~1 hour)
#   a(15) = 129125 (took ~5 hours)

#`{

# PARI/GP programs:

# Generate Lucas-Carmichael numbers <= n (slow)

isok(n) = my(f=factor(n)); for(k=1, #f~, if ((n+1) % (f[k,1]+1) != 0 || f[k,2] >= 2, return(0))); #f~ > 1;
upto(n) = my(A=List()); forprime(p=3, sqrtint(n+1)-1, my(t = p*(p+1), u = p+t); while(u <= n, if(isok(u), listput(A, u)); u += t)); vecsort(Set(A));

# Count the number of Lucas-Carmichael numbers <= n (slow)

islc(n) = my(f=factor(n)); for(k=1, #f~, if ((n+1) % (f[k,1]+1) != 0 || f[k,2] >= 2, return(0))); #f~ > 1;
a(n) = my(A=List()); forprime(p=3, sqrtint(n+1)-1, my(t = p*(p+1), u = p+t); while(u <= n, if(islc(u), listput(A, u)); u += t)); #Set(A);

# Count the number of Lucas-Carmichael numbers <= n (fast)

lucas_carmichael(A, B, k) = A=max(A, vecprod(primes(k+1))\2); (f(m, l, lo, k) = my(list=List()); my(hi=min(sqrtint(B+1)-1, sqrtnint(B\m, k))); if(lo > hi, return(list)); if(k==1, 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)))), 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); f(1, 1, 3, k);
a(n) = my(N=10^n); my(count=0); for(k=3, oo, if(vecprod(primes(k+1))\2 > N, break); count += #lucas_carmichael(1, N, k)); count;

# Generate Lucas-Carmichael numbers <= n (fast)

lucas_carmichael(A, B, k) = A=max(A, vecprod(primes(k+1))\2); (f(m, l, lo, k) = my(list=List()); my(hi=min(sqrtint(B+1)-1, sqrtnint(B\m, k))); if(lo > hi, return(list)); if(k==1, 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)))), 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); f(1, 1, 3, k);
upto(n) = my(list=List()); for(k=3, oo, if(vecprod(primes(k+1))\2 > n, break); list=concat(list, lucas_carmichael(1, n, k))); vecsort(Vec(list));

}

var min = 1
var max = 10

var n = 0
var count = 0

loop {

    for k in (3..Inf) {
        if (pn_primorial(k+1)/2 > max) {
            break
        }
        count += k.lucas_carmichael(min, max).len
    }

    say "a(#{++n}) = #{count}"

    (min, max) = (max, max*10)
}

__END__
a(1) = 0
a(2) = 0
a(3) = 2
a(4) = 8
a(5) = 26
a(6) = 60
a(7) = 135
a(8) = 323
a(9) = 791
a(10) = 1840
a(11) = 4216
a(12) = 9967
a(13) = 23070
a(14) = 54356
a(15) = 129125