Carmichael numbers in a window
An old library checks primality with the Fermat test: it takes a base a and computes an−1 mod n. Its authors are compiling a list of numbers this test is fooled by for every base coprime to n, the Carmichael numbers, to add them to a blacklist. A composite n is a Carmichael number if an−1 ≡ 1 (mod n) for all a coprime to n. Find all Carmichael numbers in several windows [L, R].
There is no need to check all bases: by Korselt’s criterion a composite n is a Carmichael number exactly when it is not divisible by the square of a prime and for each of its prime divisors p the number p − 1 divides n − 1. For example, 561 = 3 · 11 · 17, and 560 is divisible by 2, 10 and 16.
Input
The first line holds the number of windows q (1 ≤ q ≤ 5). Each of the next q lines holds two integers L and R (1 ≤ L ≤ R ≤ 10¹², R − L < 50,000). Both ends belong to the window.
Output
For each window, one line: first the number of Carmichael numbers in the window, then the numbers themselves in increasing order, separated by spaces.
Examples
Example 1
1 1 3000
5 561 1105 1729 2465 2821
Limits
- Time · Python
- 6000 ms
- Time · Go
- 1500 ms
- Memory
- 256 MB
- Judging
- 7 tests
Walkthrough
Open it step by step: each step gives away a little more. You can jump straight to the solution, but a problem solved with one hint stays with you longer.
1Hint 1
Factoring every number of the window by trial division is expensive: the window lies near 10¹², and the square root is a million. Which fast test, certainly passed by every Carmichael number, can weed out almost the whole window at once?
2Hint 2
Carmichael numbers are odd and pass the Fermat test to base 2. After such a filter a window of 50,000 numbers near 10¹² keeps about 1800 numbers, and almost all of them are prime. How do you discard the primes quickly and without errors?
3Hint 3
Watch the whole criterion. The square of the prime 1093 passes the Fermat test to base 2, and 1092 divides 1093² − 1. But it is not a Carmichael number.
4The idea
A Carmichael number is odd (for an even n the base n − 1 gives (−1)n−1 = −1) and coprime to two, so it passes the Fermat test to base 2. So first the filter: 2n−1 mod n = 1, one fast exponentiation per number. From the survivors throw out the primes with a deterministic Miller–Rabin test: for n < 1.12·10¹² the bases 2, 13, 23 and 1662803 suffice. What remains are composite pseudoprimes to base 2, a handful for the whole window and more often none, and only they are factored by trial division and checked against Korselt’s criterion.
5Walkthrough
Go over the odd n of the window. Raise 2 to the power d, where n − 1 = 2sd, and square s times: the last number of the chain is 2n−1, on which the Fermat filter works, and the whole chain gives the strong test to base 2 at the same time. The reference solution uses this so as not to exponentiate an extra time for a prime.
If n passes the filter and the strong test to bases 2, 13, 23, 1662803, it is prime. Otherwise factor it by trial division: divisors are sought up to √n ≤ 10⁶, and this has to be done rarely. The number qualifies if its factorisation has no repeated primes, there are at least three of them, and each p − 1 divides n − 1.
Why at least three prime divisors? If n = pq and p < q, then n − 1 = p(q − 1) + (p − 1), and from q − 1 | n − 1 it would follow that q − 1 | p − 1, while p − 1 is smaller. Squares of primes are excluded by a separate condition of the criterion.
5+Solution
6Complexity
For each number of the window, one exponentiation, O(log n) multiplications; for each prime, three more rounds. Trial factoring, up to 10⁶ divisions, falls only to pseudoprimes to base 2, and there are only about a hundred thousand of them per trillion numbers up to 10¹². In all, O((R − L) log R) per window.
7Common mistakes
- The criterion without “not divisible by a square”. 1093² = 1,194,649 and 3511² = 12,327,121 are pseudoprimes to base 2, and their p − 1 divides n − 1. The primes 1093 and 3511 are the only known Wieferich primes, for which 2p−1 ≡ 1 even modulo p².
- Pseudoprimes instead of Carmichael numbers. 341 = 11 · 31 passes the Fermat test to base 2, but 3340 ≢ 1 (mod 341).
- The strong test instead of the Fermat test in the filter. Only a few Carmichael numbers pass the strong test to base 2, 43 of 646 below a billion, and the filter would throw out almost all the numbers needed.
- Overflow. In Go the product of two remainders modulo about 10¹² reaches 10²⁴ and does not fit in 64 bits. A 128-bit product is needed,
bits.Mul64andbits.Div64, as in the reference solution.
Algorithms behind it
O(k log³ n)
2 problems
NUM
Algorithm
Fast exponentiation
aⁿ in O(log n) multiplications: squares over the bits of the exponent.
O(log n) multiplications
2 problems