This is about the problem that drew controversy in a recent Div. 1 round because it had more accepted solutions than Problem F. That was surprising since G was rated 4000+, while F was intended to be much easier.

This solution path is not a polished textbook derivation. Instead, I brute-forced small cases, noticed patterns in the resulting fractions, and then reverse-engineered a clean formula.

Problem in simple terms

You are given an integer m. For each number x, define a sequence like this: - Start with b(0) = 1 - Then repeatedly compute b(i+1) = x raised to the power b(i), taken modulo m We call x "good" if this sequence eventually becomes 1 and then stays 1 forever. The question: what fraction of all positive integers x are good? (That fraction is called the density of good x.)

Step 1: Start by brute forcing small primes

Begin with a small prime modulus p. For each x = 1, 2, … up to a few hundred or thousand, simulate the sequence: b(0) = 1 b(i+1) = x^b(i) mod p Run it for enough steps and check if it eventually stays at 1. Count how many x work and divide by the total tested to get a fraction. You will see neat fractions: p = 3 → 1/2

p = 5 → 1/2

p = 7 → 5/14

p = 11 → 27/110

p = 13 → 25/78

Step 2: Notice the denominator pattern

In every case, the denominator is the prime p itself times all the distinct prime factors of p − 1. Examples: For p = 11, denominator = 11 × 2 × 5 = 110 For p = 13, denominator = 13 × 2 × 3 = 78 For p = 23, denominator = 23 × 2 × 11 = 506 So the denominator rule is straightforward: p × product of distinct primes dividing p − 1.

for (auto &kv : fac) {
    int p = kv.first, e = kv.second;
    if (p == 2) continue;
    ans = (ans * ((p-1) % MOD)) % MOD;        // multiply by (p-1)
    ans = (ans * modpow(invmod(p), e)) % MOD; // divide by p^e
    ...
}

This *(p-1) and /p^e over every odd prime power p^e in m builds the “modulus and distinct-primes” denominator piece when all parts are multiplied together. (The p==2 case is handled later; see Step 7.)

Step 3: Work out the numerator from p−1

Strip off that denominator; what remains in the numerator (for prime p) matches multiply (q^e + q − 1) for each prime power q^e dividing p−1. the numerators looked like this:

p = 11 → 27

p = 13 → 25

p = 23 → 63

p = 29 → 65

p = 41 → 81

My first instinct was to search for this sequence on OEIS to see if it was known — but it didn’t appear there. So I had to look for the pattern myself.

Factoring p − 1 turned out to be the key. If p − 1 has a prime power q^e, compute (q^e + q − 1). Multiply these together for all prime powers. This exactly matched the numerators.

Note: This is the kind of pattern you may only notice by chance (and with a bit of luck) while playing with enough examples, but once you see it, the formula becomes much clearer.

Examples:

p = 11, p − 1 = 2 × 5 → (2+2−1)(5+5−1) = 3 × 9 = 27

p = 13, p − 1 = 2^2 × 3 → (4+2−1)(3+3−1) = 5 × 5 = 25

p = 23, p − 1 = 2 × 11 → (2+2−1)(11+11−1) = 3 × 21 = 63

map<int,int> pf;
factorize_value(p-1, pf);
for (auto &t : pf) S[t.first] += t.second; // sum v_q(p-1) into S[q]

Later, the second big loop over S turns these exponents into the (q^e + q − 1)-style factors prime-by-prime.

Step 4: Prime powers scale cleanly

Checking prime powers like p^2 or p^3 showed that each higher power simply divides the density by an extra factor of p. So D(p^k) = D(p) / p^(k−1).

Step 5: Combine two primes (and beyond)

For m = p^a * q^b, the denominator becomes the modulus itself ( p^a q^b ) times the distinct primes dividing (p−1)(q−1). For the numerator, add exponents per prime across all (p_i − 1) and apply the same “per prime power” factor rule.

for (auto &t : pf) S[t.first] += t.second; // S[q] = sum of v_q(p_i-1)

Step 6: product of density(x) and x only depends on the prime factors of x

For two primes p and q, tou will find that

D(p^a * q^b) × (p^a * q^b) = D(p * q) × (p * q), for all a ≥ 1 and b ≥ 1.

In other words, once you multiply the density by the modulus, the result depends only on the distinct prime factors of the modulus, not on their exponents.

Step 7: Generalize to many primes

For general m made of several prime powers:

Denominator = m × (product of distinct primes dividing all the (p − 1)).

Numerator = product over prime powers of (q^e + q − 1), where e is the total exponent gathered from all the (p − 1).

Since multiplying by m cancels out the exponents, the whole expression depends only on the prime factors of m.

Step 8: Why look at D(m) × m

Because the denominator always includes m, it’s simpler to study D(m) × m. This shows clearly that the “extra” part of the denominator is just the product of distinct primes dividing the (p − 1), while the numerator comes from the (q^e + q − 1) recipe.

Step 9: Final observation (for odd x)

If x is odd, then multiplying x by any power of 2 does not change the scaled quantity D(·)×(·):

D(2^t · x) × (2^t · x) = D(x) × x for every t ≥ 0.

In other words, once you multiply the density by the modulus, the result depends only on the set of prime factors of the modulus, not on their exponents. So adding more factors of 2 (or raising any odd prime factor to higher powers) doesn’t change D(m)×m—only D(m) itself shrinks by the corresponding prime power while the product stays the same.

Use of AI

GitHub Copilot or cursor will autocomplete 80% of the boilerplate before you’ve even finished thinking it.

Obvious functions like modpow, modinv, sieve (for computing smallest primes), and factorization.

Guarding against overflow errors: e.g. if (x >= MOD) x -= MOD; after modular addition, or using __int128 for modular multiplication.

That separation of labor makes it possible for me to spot the structure without getting bogged down.

A loosely related projecteuler problem from 17 years ago:

Disclaimer: Noticing the right patterns doesn’t always happen, and this simulate→observe→generalize approach may not work on other problems. Still, it’s a nice example of how you can sometimes chip away at very tough-looking problems without “Terence Tao–level math.”