Solution

Answer: 84664213

Let

$$f(G,L,N)=#{(a_1,\dots,a_N)\in \mathbb N^N:\gcd(a_i)\ge G,\ \operatorname{lcm}(a_i)\le L}.$$

We must compute

$$f(10^6,10^{12},10^{18}) \pmod{101^4}.$$

The modulus is

$$101^4=104060401.$$

1. Mathematical analysis

Step 1: Normalize by the gcd

Let

$$g=\gcd(a_1,\dots,a_N).$$

Write

$$a_i=g b_i.$$

Then:

$\gcd(b_1,\dots,b_N)=1$,
$\operatorname{lcm}(a_i)=g\operatorname{lcm}(b_i)$.

Define

$$d=\operatorname{lcm}(b_i).$$

Then every valid list corresponds uniquely to:

an integer $g\ge G$,
an integer $d\ge 1$,
a list $(b_i)$ with

$$\gcd(b_i)=1,\qquad \operatorname{lcm}(b_i)=d,$$

and the constraint

$$gd\le L.$$

Hence for fixed $d$, the number of allowed $g$ is

$$\left\lfloor \frac{L}{d}\right\rfloor-G+1,$$

provided $d\le L/G$.

Therefore

$$f(G,L,N) = \sum_{d\le L/G} h_N(d)\left(\left\lfloor\frac{L}{d}\right\rfloor-G+1\right),$$

where $h_N(d)$ is the number of length-$N$ lists with:

$$\gcd=1,\qquad \operatorname{lcm}=d.$$

So the whole problem reduces to computing $h_N(d)$.

Step 2: Compute $h_N(d)$

Factor

$$d=\prod p^{e_p}.$$

Because gcd/lcm conditions separate prime-by-prime, we only need to study one prime power.

Suppose for a prime $p$, the exponents in the $N$ entries are

$$x_1,\dots,x_N\in{0,1,\dots,e}.$$

To have:

lcm exponent $=e$: at least one $x_i=e$,
gcd exponent $=0$: at least one $x_i=0$.

Thus we count exponent-tuples with:

$$\max x_i=e,\qquad \min x_i=0.$$

Total tuples in $[0,e]$:

$$(e+1)^N.$$

Subtract tuples with no zero:

$$e^N.$$

Subtract tuples with no $e$:

$$e^N.$$

Add back tuples with neither $0$ nor $e$:

$$(e-1)^N.$$

Therefore the contribution of exponent $e$ is

$$A_e=(e+1)^N-2e^N+(e-1)^N.$$

Since different primes are independent:

$$h_N(d)=\prod_{p^e\parallel d} A_e.$$

This is multiplicative.

Step 3: Final summation

Let

$$M=\left\lfloor \frac{L}{G}\right\rfloor.$$

Then

$$f(G,L,N) = \sum_{d=1}^{M} \left( \prod_{p^e\parallel d} \big((e+1)^N-2e^N+(e-1)^N\big) \right) \left( \left\lfloor \frac{L}{d}\right\rfloor-G+1 \right).$$

For this problem:

$$G=10^6,\quad L=10^{12},\quad N=10^{18},$$

$$M=10^6.$$

The maximum exponent appearing in numbers $\le10^6$ is only $19$, since

$$2^{19}=524288,\qquad 2^{20}>10^6.$$

Thus we precompute $A_e$ for $e=0,\dots,19$ modulo $101^4$, then sieve multiplicatively up to $10^6$.

Complexity:

sieve: $O(M\log\log M)$,
evaluation: $O(M)$,

which is easily feasible.

2. Python implementation

MOD = 101**4

G = 10**6
L = 10**12
N = 10**18

M = L // G   # = 10^6

# ------------------------------------------------------------
# Precompute:
# A[e] = (e+1)^N - 2*e^N + (e-1)^N  (mod MOD)
# ------------------------------------------------------------

MAX_E = 20

A = [0] * (MAX_E + 1)

for e in range(MAX_E + 1):
    A[e] = (
        pow(e + 1, N, MOD)
        - 2 * pow(e, N, MOD)
        + pow(e - 1, N, MOD)
    ) % MOD

# ------------------------------------------------------------
# Smallest prime factor sieve
# ------------------------------------------------------------

spf = list(range(M + 1))

for i in range(2, int(M**0.5) + 1):
    if spf[i] == i:  # i is prime
        for j in range(i * i, M + 1, i):
            if spf[j] == j:
                spf[j] = i

# ------------------------------------------------------------
# Compute h[d] multiplicatively
#
# h[d] = product over p^e || d of A[e]
# ------------------------------------------------------------

h = [0] * (M + 1)
h[1] = 1

for n in range(2, M + 1):

    p = spf[n]

    e = 0
    t = n

    while t % p == 0:
        t //= p
        e += 1

    h[n] = (h[t] * A[e]) % MOD

# ------------------------------------------------------------
# Final summation
# ------------------------------------------------------------

answer = 0

for d in range(1, M + 1):

    weight = (L // d) - G + 1

    answer = (answer + h[d] * weight) % MOD

print(answer)

3. Code walkthrough

Precomputing $A_e$

A[e] = (e+1)^N - 2e^N + (e-1)^N

This is exactly the inclusion–exclusion count for exponent tuples whose minimum is $0$ and maximum is $e$.

We compute all powers modulo $101^4$.

Smallest-prime-factor sieve

The spf array stores the smallest prime factor of every number up to $10^6$.

This lets us factor each number very quickly.

Building $h[n]$

$$n = p^e \cdot t$$

with $p\nmid t$, then multiplicativity gives

$$h(n)=h(t)\cdot A_e.$$

So each value is computed in constant amortized time.

Final sum

For each $d$:

weight = floor(L/d) - G + 1

counts the valid gcd choices $g$.

We accumulate

$$h(d)\times \text{weight}.$$

4. Verification on the examples

The program reproduces all sample values:

$$f(10,100,1)=91$$

$$f(10,100,2)=327$$

$$f(10,100,3)=1135$$

and

$$f(10,100,1000)\bmod 101^4 = 3286053.$$

So the derivation is consistent.

5. Final computation

Running the algorithm for

$$(G,L,N)=(10^6,10^{12},10^{18})$$

gives

$$84664213.$$

Answer: 84664213

Project Euler Problem 350