I'm trying to solve the problem. I need to find a Carmichael number which is a product of seven prime numbers, each between 10^7 and 10^9. Is there any way to do it?
I tried to solve this task using Chernick's formula:
M(m) = (6m+1)(12m+1)(18m+1)(36m+1)(72m+1)(144m+1)(288*m+1),
on condition that all factors are prime and m is divisible by 8.
Also I used Korselt's criterion: a positive composite integer n is a Carmichael number if and only if n is square-free, and for all prime divisors p of n, it is true that (n-1) is divisible by (p-1). Here is the code:
from gmpy2 import is_prime
from math import prod
m = 1666672
while True:
x = [i * m + 1 for i in [6, 12, 18, 36, 72, 144, 288]]
n = prod(x)
if all(is_prime(i) and (n - 1) % (i - 1) == 0 and i in range(10 ** 7, 10 ** 9 + 1) for i in x):
print(m, n)
break
elif x[6] > 10 ** 9:
print('Nothing found')
break
m += 8
Finally, I got m = 18726360 and n = 112504155402956714375844566364658214310419726067438239203656961.
Prime factors of n are: 112358161, 224716321, 337074481, 674148961, 1348297921, 2696595841, 5393191681.
But only the first four of them are in the required range [10^7, 10^9], so my solution didn't pass :(
Lets first approach this problem by exploring Korselt's Criterion a bit. As you state, n is a Carmichael Number if and only if:
We achieve (1) by making each of the 7 prime factors distinct. That leaves (2). This means n - 1 is a multiple of each p - 1, the smallest such value being the Least Common Denominator (LCM) of the p - 1 values. The probability that a truly random n - 1 is divisible by some other integer K is 1/K. This gives us the intuition that we should try to minimize our LCM, to increase our chances that it divides n - 1. Of course n - 1 will not be truly random, but the heuristic still works.
To minimize the LCM of p - 1 values, we'll generate values that contain intentionally contain many of the same prime factors. In fact, we'll generate 5-smooth numbers, meaning their largest prime factor is at most 5. This way, at worst we a LCM that dominated by three prime powers: 2^x, 3^y, 5^z. In practice, we'll avoid this worst case; most of our p - 1 values will be composed of some "balanced" prime powers, resulting in a smaller LCM.
This alone would be enough to brute force a solution, but we have one more trick we can use to speed up the search further. Let each prime factor p of n be one greater than a multiple of some constant k. Symbolically: p_i = a_i * k + 1. We notice that k divides n - 1:
n = p_0 * p_1 * ... * p_m
n = (a_0 * k + 1) * (a_1 * k + 1) * ... * (a_m * k + 1)
n = k * [...] + 1
n - 1 = k * [...]
With this information, we can actually pick a handful of small, distinct 5-smooth numbers, and some larger common factor k, and multiply them for our p - 1 values. We're guaranteed that k divides n - 1, and the 5-smooth numbers (and their LCM) are much smaller, greatly increasing our chances of their LCM evenly dividing n - 1.
Putting this all into code:
from queue import PriorityQueue
from itertools import combinations
from math import prod
from gmpy2 import is_prime
def gen_5_smooth():
primes = [2, 3, 5]
queue = PriorityQueue()
for i, p in enumerate(primes):
queue.put((p, i))
yield 1
while True:
x, i = queue.get()
yield x
for j, q in enumerate(primes[i:], i):
queue.put((q*x, j))
def solve():
limit_lower = 10**7
limit_upper = 10**9
smooth_lower = 1000
smooth_upper = smooth_lower * 10
mult_lower = limit_lower // smooth_lower + 1
mult_upper = limit_upper // smooth_upper
smooth_factors = []
for fact in gen_5_smooth():
if fact >= smooth_upper:
break
if fact >= smooth_lower:
smooth_factors.append(fact)
for mult in range(mult_lower, mult_upper):
primes = [fact * mult + 1 for fact in smooth_factors if is_prime(fact * mult + 1)]
for comb in combinations(primes, 7):
product = prod(comb)
if all((product - 1) % (x - 1) == 0 for x in comb):
yield comb
I chose the smooth bounds manually based on the number of 5-smooth numbers that fell into range and the size of the corresponding multipliers to push them into the [1e7, 1e9) range. You could adjust that value for potentially faster results. You could also generate combinations in order of increasing LCM with some extra effort, which would also produce solutions faster. However, the above code produces a couple solutions per second on average, so its probably fast enough.
Iterating over the first handful of solutions yields:
(13501351, 14581459, 16201621, 17281729, 58325833, 64806481, 77767777)
(12070801, 17381953, 23175937, 25147501, 45265501, 50295001, 57939841)
(13606651, 14513761, 21770641, 22677751, 36284401, 48983941, 87082561)
(10991161, 14654881, 17585857, 35171713, 36637201, 54955801, 58619521)
(12268801, 15336001, 16562881, 23556097, 26500609, 47112193, 70668289)
(13269761, 17914177, 37321201, 59713921, 62202001, 82936001, 95542273)