Here's an algorithm in Python code for efficiently computing the Carmichael numbers in a given range:
from sympy.ntheory import sieve, isprime, prime
from sympy.core import integer_nthroot
from math import lcm, gcd, isqrt
def find_carmichael_numbers_in_range(A, B):
max_p = 1+((1 + isqrt(8*B + 1)) >> 2)
def _rec(m, lam, lo, k):
if k == 1:
lo = max(lo, A // m + (1 if A % m else 0))
hi = min(B // m + 1, max_p)
u = pow(m, -1, lam)
j,r = divmod(lo - u, lam)
if r: j += 1
u += j*lam
for p in range(u, hi, lam):
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 _rec(m*p, lcm(lam, p - 1), p + 2, k - 1)
def _rec1():
k = 3
l = 3*5*7
while l < B:
yield from _rec(1, 1, 3, k)
k += 1
l *= prime(k+1)
return sorted(_rec1())
assert find_carmichael_numbers_in_range(561, 2000) == [561, 1105, 1729]
assert find_carmichael_numbers_in_range(562, 2000) == [1105, 1729]
assert find_carmichael_numbers_in_range(10585, 75361) == [10585, 15841, 29341, 41041, 46657, 52633, 62745, 63973, 75361]
assert [len(find_carmichael_numbers_in_range(1, 10**n)) for n in range(1, 8)] == [0, 0, 1, 7, 16, 43, 105] # A055553
print(find_carmichael_numbers_in_range(1010, 1010 + 10**8)) #=> [10017089857, 10017782401, 10028287081, 10031651841, 10054063041, 10073377801, 10087556401]