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, it is a poor man’s, non-mathematician programmer’s way: brute force small cases, notice patterns in the resulting fractions, and then reverse-engineer 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 ^ b(i) % 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.

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.”