Problem 350: Constraining the Least Greatest and the Greatest Least

Project Euler Problem 350 Solution

EulerSolve provides an optimized solution for Project Euler Problem 350, Constraining the Least Greatest and the Greatest Least, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We must count positive integer $N$-tuples $(x_1,\dots,x_N)$ such that $$\gcd(x_1,\dots,x_N)\ge G,\qquad \operatorname{lcm}(x_1,\dots,x_N)\le L,$$ and return the result modulo $101^4$. A direct search is impossible for the actual parameters, so the solution reorganizes each tuple by separating its common gcd from its normalized lcm structure. Mathematical Approach Let $$F(G,L,N)=\#\left\{(x_1,\dots,x_N)\in \mathbb{Z}_{>0}^N:\gcd(x_1,\dots,x_N)\ge G,\ \operatorname{lcm}(x_1,\dots,x_N)\le L\right\}.$$ The key idea is to classify tuples by their common gcd and by the lcm that remains after dividing that gcd out. Step 1: Normalize by the Common GCD For any admissible tuple define $$g=\gcd(x_1,\dots,x_N),\qquad x_i=g\,y_i.$$ Then the normalized tuple satisfies $$\gcd(y_1,\dots,y_N)=1.$$ Now define $$m=\operatorname{lcm}(y_1,\dots,y_N).$$ Because lcm scales linearly when every coordinate is multiplied by the same factor, we obtain $$\operatorname{lcm}(x_1,\dots,x_N)=g\,m.$$ So once the normalized tuple is fixed, the only remaining freedom is the choice of $g$....

Detailed mathematical approach

Problem Summary

A list of size $n$ is a sequence of $n$ natural numbers, so the order matters: $(2,4,6)$ and $(2,6,4)$ are two different lists. The $\gcd$ of a list is the largest natural number dividing every entry, the $\operatorname{lcm}$ the smallest natural number divisible by every entry; for $(10,6,15,6)$ they are $1$ and $30$. Let $f(G,L,N)$ be the number of lists of size $N$ whose entries satisfy

$$\begin{aligned} \gcd(x_1,\dots,x_N)&\ge G,\\ \operatorname{lcm}(x_1,\dots,x_N)&\le L. \end{aligned}$$

The statement gives $f(10,100,1)=91$, $f(10,100,2)=327$, $f(10,100,3)=1135$ and $f(10,100,1000)\equiv 3286053$ modulo $101^4$, and asks for $f(10^6,10^{12},10^{18})$ modulo $101^4=104060401$. ^[1]

The whole solution in one sentence: divide every list by its own $\gcd$, so that each list becomes a scale factor times a list with no common factor, and count the latter one prime at a time, because for a fixed least common multiple the number of admissible scale factors depends on that multiple alone.

Mathematical Approach

1. Counting a tiny case by hand

Start with $f(4,12,2)$: pairs $(x_1,x_2)$ with $\gcd\ge4$ and $\operatorname{lcm}\le12$. Every entry divides the $\operatorname{lcm}$, so no entry exceeds $12$, and a short search by hand finds exactly three ratio types.

Equal entries $(g,g)$ with $4\le g\le 12$: here $\gcd=\operatorname{lcm}=g$, giving $9$ lists.
Entries in ratio $1:2$, that is $(g,2g)$ or $(2g,g)$: the $\operatorname{lcm}$ is $2g\le12$, so $g\in\{4,5,6\}$ and we get $6$ lists.
Entries in ratio $1:3$: the $\operatorname{lcm}$ is $3g$, and $3g\le12$ with $g\ge4$ forces $g=4$, giving $(4,12)$ and $(12,4)$.

So $f(4,12,2)=9+6+2=17$, which an exhaustive search over all $12^2$ pairs confirms. Two features of this count drive the whole solution: every list was written as a scale factor $g$ times a shape such as $(1,1)$, $(1,2)$ or $(3,1)$ whose entries have no common factor, and the number of admissible scale factors depended only on the $\operatorname{lcm}$ of that shape.

2. Dividing out the greatest common divisor

Let $x=(x_1,\dots,x_N)$ be any list and $g=\gcd(x)$. Each $x_i$ is a multiple of $g$, so $x_i=g\,y_i$ with positive integers $y_i$, and

$$\gcd(y_1,\dots,y_N)=1,$$

because $\gcd(g\,y_1,\dots,g\,y_N)=g\cdot\gcd(y)$. ^[2] A list with $\gcd=1$ will be called primitive: no prime divides every entry. Its entries need not be pairwise coprime — $(6,10,15)$ is primitive although every pair shares a prime. The $\operatorname{lcm}$ scales the same way. Writing $v_p(t)$ for the exponent of $p$ in $t$, unique factorisation ^[3] gives $v_p(\operatorname{lcm})=\max_i v_p(x_i)$, so

$$v_p\bigl(\operatorname{lcm}(x)\bigr)=v_p(g)+\max_i v_p(y_i),$$

since $v_p(x_i)=v_p(g)+v_p(y_i)$ for every $i$. Doing this for every prime yields

$$\operatorname{lcm}(x)=g\cdot\operatorname{lcm}(y).$$

Conversely a scale factor $g\ge1$ and a primitive list $y$ rebuild $x_i=g\,y_i$, whose $\gcd$ is $g$ again, so lists correspond exactly to pairs of a scale factor and a primitive list: nothing is counted twice.

3. How many scale factors fit

Fix a primitive list $y$ and let $m=\operatorname{lcm}(y)$. By Section 2 the pair $(g,y)$ is admissible exactly when $g\ge G$ and $g\,m\le L$, that is when

$$G\le g\le\lfloor L/m\rfloor .$$

Put $K=\lfloor L/G\rfloor$. This interval contains an integer exactly when $m\le K$: if $m\le L/G$ then $mG\le L$, hence $G\le L/m$ and, $G$ being an integer, $G\le\lfloor L/m\rfloor$; while if $m\ge K+1$ then $mG\gt L$, so $\lfloor L/m\rfloor\lt G$ and no scale factor fits at all. For $m\le K$ the number of choices is therefore

$$w(m)=\lfloor L/m\rfloor-G+1\ge1,$$

which depends only on $m$. Writing

$$H_N(m)=\#\{y:\gcd(y)=1,\ \operatorname{lcm}(y)=m\}$$

for the number of primitive lists of size $N$ with $\operatorname{lcm}$ exactly $m$, grouping the pairs by $m$ gives the master formula

$$f(G,L,N)=\sum_{m=1}^{K}H_N(m)\,w(m).$$

For $(G,L,N)=(4,12,2)$ this reads $K=3$ and $H_2(1)\,w(1)=1\cdot9$, $H_2(2)\,w(2)=2\cdot3$, $H_2(3)\,w(3)=2\cdot1$, total $9+6+2=17$ — exactly the three ratio types found by hand in Section 1.

4. A primitive list is a table of exponents

Every entry divides its $\operatorname{lcm}$, so $\operatorname{lcm}(y)=m$ forces $y_i\mid m$ for all $i$. Let $m=\prod_p p^{e_p}$ be the canonical factorisation. ^[3] Each $y_i$ is then determined by its exponents $a_{i,p}=v_p(y_i)$, which range over $\{0,1,\dots,e_p\}$, and any such table of exponents comes from exactly one list of divisors of $m$. For instance $(2,3,6)$ has $m=6$, with column $(1,0,1)$ at $2$ and $(0,1,1)$ at $3$; both have minimum $0$ and maximum $1$. In general the two requirements hold together exactly when every prime $p$ dividing $m$ satisfies both of

$$\min_i a_{i,p}=0,$$

$$\max_i a_{i,p}=e_p.$$

The first family of equalities says $\gcd(y)=1$, because $v_p(\gcd)=\min_i a_{i,p}$ and a number is $1$ exactly when all its prime exponents vanish; the second says $\operatorname{lcm}(y)=m$, because $v_p(\operatorname{lcm})=\max_i a_{i,p}$. ^[2] Neither family may be split prime by prime — $(2,6)$ has minimum $0$ in its column at $3$ yet $\gcd=2$ — but each single equality involves one prime only, so the columns can be chosen independently.

5. One prime power by inclusion and exclusion

Fix a prime power $p^{e}$, $e\ge1$. By Section 4 we must count the vectors $(a_1,\dots,a_N)\in\{0,1,\dots,e\}^N$ with minimum $0$ and maximum $e$. Let $A$ be the vectors avoiding the value $0$ and $B$ those avoiding $e$; we want the vectors in neither, so the complementary form of inclusion–exclusion ^[4] applies:

$$C_e=(e+1)^N-|A|-|B|+|A\cap B|.$$

There are $(e+1)^N$ vectors in all. A vector in $A$ uses only $\{1,\dots,e\}$, so $|A|=e^N$; one in $B$ uses only $\{0,\dots,e-1\}$, so $|B|=e^N$ too; one in both uses only $\{1,\dots,e-1\}$, so $|A\cap B|=(e-1)^N$. Hence

$$C_e=(e+1)^N-2e^N+(e-1)^N,$$

where $0^N=0$ because $N\ge1$. For $e=1$ this is $C_1=2^N-2$: the vectors over $\{0,1\}$ holding both a $0$ and a $1$ are all but the two constant ones. The expression is the second difference of $t\mapsto t^N$, so the small cases are easy to read off; checked for every $e\le50$, it is identically $0$ for $N=1$, identically $2$ for $N=2$ and $6e$ for $N=3$. The value $0$ at $N=1$ is exactly right: a list of size $1$ has $\gcd=\operatorname{lcm}$, so the only primitive one is $(1)$, and the master formula keeps the term $m=1$ alone and returns $L-G+1$ whenever $L\ge G$ — which is $91$ for $(G,L)=(10,100)$.

6. The count for one prime power determines all of them

Because the columns may be chosen independently, the number of primitive lists with $\operatorname{lcm}$ exactly $m$ is the product of the per-column counts:

$$H_N(m)=\prod_{p^{e}\,\|\,m}C_e,$$

the product running over the prime powers dividing $m$ exactly. In particular $H_N(1)=1$ (empty product) and $H_N$ is multiplicative: $H_N(ab)=H_N(a)H_N(b)$ for coprime $a,b$, since coprime arguments contribute disjoint sets of prime powers. ^[5] Two specialisations are worth seeing. For $N=2$ every factor is $2$, so $H_2(m)=2^{\omega(m)}$ with $\omega(m)$ the number of distinct primes dividing $m$; its first ten values $1,2,2,2,2,4,2,2,2,4$ are the first ten terms of the sequence counting unitary divisors. ^[6] For $N=3$ we get $H_3(m)=\prod_p 6e_p$.

7. Worked example: the checkpoints at $(G,L)=(10,100)$

Here $K=\lfloor100/10\rfloor=10$, and the weights $w(m)=\lfloor100/m\rfloor-9$ for $m=1,\dots,10$ are

$$91,41,24,16,11,7,5,3,2,1 .$$

For $N=2$ Section 6 gives $H_2(m)=2^{\omega(m)}$, so the ten values of $H_2$ are $1,2,2,2,2,4,2,2,2,4$ and the products $H_2(m)\,w(m)$ are

$$91,82,48,32,22,28,10,6,4,4,$$

which add up to $327$. For $N=3$ the factors are $C_e=6e$, so $H_3(m)=1,6,6,12,6,36,6,18,12,36$ — note $H_3(4)=C_2=12$, $H_3(6)=C_1C_1=36$ — and the products are

$$91,246,144,192,66,252,30,54,24,36,$$

with sum $1135$. Both agree with the statement, and the same ten weights reproduce $91$ for $N=1$ and $3286053$ for $N=1000$.

8. A second route to $H_N$, used as a check

An independent derivation is useful. Let $\Lambda_N(m)$ count the lists of divisors of $m$ with $\operatorname{lcm}$ exactly $m$ and no condition on the $\gcd$: only the maximum condition of Section 4 survives, so inclusion–exclusion with one bad set gives

$$\Lambda_N(m)=\prod_{p^{e}\,\|\,m}\bigl((e+1)^N-e^N\bigr).$$

On the other hand Section 2 splits any list with $\operatorname{lcm}=m$ into its $\gcd$, which divides $m$, and a primitive list whose $\operatorname{lcm}$ is the quotient, so

$$\Lambda_N(m)=\sum_{d\mid m}H_N(d).$$

Möbius inversion undoes such a sum. Recall $\mu(1)=1$, $\mu(d)=0$ when a prime square divides $d$, and $\mu(d)=(-1)^r$ when $d$ is a product of $r$ distinct primes; then ^[7]

$$H_N(m)=\sum_{d\mid m}\mu(d)\Lambda_N(m/d).$$

For instance $\Lambda_2(4)=9-4=5$ and $\Lambda_2(2)=4-1=3$ give $H_2(4)=5-3=2$, which is $C_2=2$ as it should be. The two descriptions of $H_N$ also agree on prime powers by a telescoping sum,

$$\sum_{j=0}^{e}C_j=(e+1)^N-e^N,$$

since the inner differences cancel in pairs and $C_0$ is read as $H_N(1)=1$. Both identities were confirmed for several $N$ and all $m\le300$, and the inversion is used again below as a full-size cross-check.

9. The size of the target computation

At the asked parameters $K=10^{12}/10^{6}=10^{6}$, so the master formula is a sum of one million terms needing only the factorisation of every $m\le10^6$. The exponents stay small: $2^{19}=524288\le10^6\lt 2^{20}$, so no prime power dividing such an $m$ has exponent above $19$, and twenty values $C_e$ suffice. They need the powers $b^N$ modulo $101^4$ for the $21$ bases $b=1,\dots,21$, and repeated squaring gets each without building the full power: $13=8+4+1$ gives $b^{13}=b^{8}b^{4}b$, with $b^2,b^4,b^8$ three successive squarings. ^[8] Since $N=10^{18}$ has $60$ binary digits, that is $60$ squarings per base. The rest is a sieve, one factorisation pass over $m$ and one weighted-sum pass.

10. Working modulo a prime power

Two cautions about the modulus $101^4=104060401$. It is a prime power, not a prime, so the residues do not form a field: $101\cdot101^3\equiv0$ exhibits zero divisors, and division is not generally available. That costs nothing here, since the master formula only adds and multiplies. Shortening the exponent is the second trap: Euler's theorem replaces $N$ by its residue modulo $\varphi(101^4)=103030100$ only for bases coprime to $101$, true of every base $1,\dots,21$ used here but saving nothing at $21$ exponentiations, so the code uses $N$ in full. Finally the weights are large — $w(1)=L-G+1=999999000001$ at the target — and all three programs reduce each modulo $101^4$ before multiplying by $H_N(m)$; in Java that reduction keeps the product inside the signed 64-bit range, while C++ and Python could also hold it unreduced.

How the Code Works

All three programs implement the master formula of Section 3 in the same four stages, with MOD set to $104060401$. In the C++ file, ^[9] mod_add and mod_sub correct by one conditional addition or subtraction instead of a remainder, mod_mul widens both factors to __int128 before reducing, ^[10] and mod_pow is the right-to-left binary exponentiation of Section 9. ^[8] The widening is a margin rather than a necessity: $\mathtt{MOD}^2=10828567056280801$ is below $2^{63}-1$, so a 64-bit product would not overflow.

build_spf(K) is the first stage. It sets spf[i] = i for all $i\le K$ and then, for each $i$ still equal to spf[i] and hence prime, overwrites the still-unmarked multiples $j=i^2,i^2+i,\dots$. Starting at $i^2$ is enough, ^[8] and writing only where spf[j] == j means the value that survives is the smallest prime factor of $j$.

The second stage builds the local factors of Section 5. A halving loop computes max_exp as the number of binary digits of $K$, which is $20$ for $K=10^6$ — one more than the largest exponent that can occur, a harmless slack. Then pows[b] is set to $b^N$ modulo MOD for $b$ up to max_exp + 1, with pows[0] = 0 because $N\ge1$, and local[a] becomes $C_a=(a+1)^N-2a^N+(a-1)^N$.

The third stage fills the table h of length $K+1$ with $H_N$. After h[1] = 1, each $n$ from $2$ to $K$ is factorised by repeatedly reading $p=$ spf[x], dividing it out while it divides, counting the exponent $e$ and multiplying local[e] into a running product — the product of Section 6, with the exponents delivered by the sieve instead of by trial division. The fourth stage is the loop of Section 3: for every $m$ it forms L / m - G + 1, reduces that weight modulo MOD as Section 10 explains, multiplies by h[m] and accumulates.

Only the C++ file carries tests: run_checkpoints calls solve four times and compares against the statement's $91$, $327$, $1135$ and $3286053$, printing which checkpoint broke and exiting with status $2$ on a failure. The flag --skip-checkpoints suppresses them; any other argument is rejected. The Python function solve_350 runs the same four stages with the three-argument pow ^[10] and arbitrary-precision integers, so it needs no overflow analysis. Java keeps the parameters, residues and weights in long but the sieve table, the bound $K$, the indices and the exponents in int; the long arithmetic is safe because every product of two residues stays below $\mathtt{MOD}^2$ and hence inside the signed 64-bit range. ^[10] In C++ the parameters are unsigned 64-bit, the residues signed 64-bit, spf an array of int. One Java line, (pows[a + 1] - 2 * pows[a] % MOD + MOD) % MOD, relies on % binding more tightly than -: it subtracts $(2\,\mathtt{pows}[a])\bmod \mathtt{MOD}$, a value below MOD, so the intermediate lies between $0$ and $2\,\mathtt{MOD}$ and the reduction is correct — confirmed against the exact value for every exponent the program uses. Neither port runs the checkpoints, so a change made only to them would not be caught by a built-in test.

Complexity and Verification

With $K=\lfloor L/G\rfloor$, the sieve costs $O(K\log\log K)$ time ^[8] and performed $2122048$ inner-loop steps at $K=10^6$, finding the $78498$ primes below the bound. Write $\Omega$ for the number of prime factors with multiplicity and $\omega$ for the number of distinct ones. For each $m$ the factorisation stage evaluates $\Omega(m)+\omega(m)$ remainder tests — one per prime factor, plus one failing test per distinct prime — performs $\Omega(m)$ divisions and multiplies in $\omega(m)$ local factors. At the target bound

$$\sum_{m\le10^6}\Omega(m)=3626619,\qquad \sum_{m\le10^6}\omega(m)=2853708,$$

so the stage makes $6480327$ remainder tests, $3626619$ divisions and $2853708$ modular multiplications. The $21$ modular powers — base $0$ is assigned directly — cost $60$ squarings each and are negligible, and the final sum makes $K$ divisions and $K$ multiplications. The factorisation stage obeys the same $O(K\log\log K)$ bound, because $\sum_{m\le K}\Omega(m)$ equals $\sum_{p^a\le K}\lfloor K/p^a\rfloor$, whose terms with $a=1$ are exactly the marking count of the sieve, while the higher powers add less than $K$ because $\sum_p 1/(p(p-1))\lt1$. With the modulus fixed the method is $O(K\log\log K)$ in time and $O(K)$ in memory, the $O(\log K)$ modular powers adding $O(\log K\log N)$ multiplications. In C++ and Java the two tables hold $4(K+1)$ bytes of 32-bit smallest prime factors, about $4$ MB, and $8(K+1)$ bytes of 64-bit residues, about $8$ MB. Measured peak resident memory was about $15$ MB for C++, $58$ MB for Java and $64$ MB for Python, whose lists hold pointers to boxed integers rather than packed words. Measured times on one desktop machine were $0.03$ s, $0.05$ s and $0.45$ s respectively; the checkpoints add nothing measurable, running at $K=10$.

Beyond the four checkpoints in the C++ file, the following checks were carried out for this article.

An exhaustive search over all lists agreed with the master formula for every triple with $1\le L\le12$, $1\le G\le L$, $1\le N\le4$ — $312$ in all — and for larger single cases such as $f(5,30,2)=74$, $f(3,18,3)=148$ and $f(1,4,5)=274$.
$H_N(m)$ counted directly, over all lists of divisors of $m$, matched the product of Section 6 for every $N\le4$ and $m\le30$.
The identities of Section 8 held for all $m\le300$ and several $N$, including $\sum_{d\mid m}\Lambda_N(d)=\tau(m)^N$ with $\tau$ the number of divisors — the statement that every list of divisors of $m$ has some divisor of $m$ as its $\operatorname{lcm}$.
At the full target size three structurally different computations agree: the published product formula; the Möbius inversion of Section 8, built from the local factor $(e+1)^N-e^N$ and summed over the divisors of every $m\le10^6$; and the total re-ordered as a sum over $g$, using prefix sums of $H_N$ and blocks on which $\lfloor L/g\rfloor$ is constant.
The C++, Python and Java programs print the same value.

One sharp edge is worth recording. The code assigns h[1] = 1 right after allocating $K+1$ entries, so a call with $L\lt G$, where $K=0$, writes outside the table: Python raises an IndexError, Java an ArrayIndexOutOfBoundsException and C++ would write past the end of a one-element vector. Every published parameter set has $L\ge G$, so this never arises as the program stands. The opposite extreme is correct but degenerate: $f(10,10,5)=1$, the list $(10,10,10,10,10)$ alone, since every entry must be at least $10$ and divide $10$.

Footnotes and References

The references below supply the standard background used above. The scale-factor reduction, the weight $w(m)$, the local count $C_e$, the two derivations of $H_N$ and all operation counts are derived or measured explicitly here.

Project Euler 350 — Constraining the Least Greatest and the Greatest Least. The official statement defines a list of size n as a sequence of n natural numbers and f(G, L, N) as the number of such lists with gcd at least G and lcm at most L, gives f(10, 100, 1) = 91, f(10, 100, 2) = 327, f(10, 100, 3) = 1135 and f(10, 100, 1000) mod 101⁴ = 3286053, and asks for f(10⁶, 10¹², 10¹⁸) mod 101⁴.
Wikipedia — Greatest common divisor gives the gcd as the product of each prime raised to the minimum of its exponents in the two factorisations, and states gcd(m·a, m·b) = m·gcd(a, b); both are used below. Wikipedia — Least common multiple gives the companion formula with the maximum of the exponents, which turns the lcm conditions into conditions on columns.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a product of primes in exactly one way up to order, and its canonical representation writes it as a product of powers of distinct primes. This makes the correspondence of Section 4 between lists of divisors of m and tables of exponents a bijection.
Wikipedia — Inclusion–exclusion principle. Besides |A ∪ B| = |A| + |B| − |A ∩ B|, the article gives the complementary form, which counts the elements having none of a list of properties by alternately subtracting and adding the sizes of the intersections. Section 5 uses the properties "no exponent is 0" and "no exponent is e".
Wikipedia — Multiplicative function. An arithmetic function is multiplicative when f(1) = 1 and f(ab) = f(a)f(b) for coprime a and b, and is then completely determined by its values at prime powers, as in Section 6.
OEIS A034444 — number of unitary divisors. The sequence equals 2 raised to the number of distinct prime factors, and its first ten terms are 1, 2, 2, 2, 2, 4, 2, 2, 2, 4, which is the table of H₂(m) computed in Sections 6 and 7.
Wikipedia — Möbius inversion formula. If g(n) is the sum of f(d) over the divisors d of n, then f(n) is the sum of μ(d)g(n/d) over those divisors. Section 8 applies it to Λ_N and H_N, and the verification section uses it as a full-size cross-check.
Wikipedia — Modular exponentiation describes the right-to-left binary method, which squares a running base once per binary digit of the exponent, multiplies it into the result for each digit 1 and reduces at every step, in time logarithmic in the exponent; this is mod_pow. Wikipedia — Sieve of Eratosthenes states that marking a prime's multiples may begin at its square, since smaller multiples carry a smaller prime factor, and that the sieve runs in O(n log log n) operations.
C++ — Euler350.cpp. The linked immutable C++ revision contains build_spf, mod_pow, solve and run_checkpoints. It is the implementation record for this article; the extra comparisons reported in the verification section are not packaged in it.
Python — built-in function pow: with a third argument it returns the power modulo that argument, computed more efficiently than taking the remainder of the full power. Java Language Specification §4.2.1: a long holds the integers from −2⁶³ to 2⁶³ − 1, which contains every product of two residues modulo 101⁴. GCC manual — 128-bit integers: the __int128 extension needs a target with a 128-bit-wide integer mode, and mod_mul uses it.

Problem 350 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;

constexpr i64 MOD = 104060401;  // 101^4

i64 mod_add(const i64 a, const i64 b) {
    i64 x = a + b;
    if (x >= MOD) {
        x -= MOD;
    }
    return x;
}

i64 mod_sub(const i64 a, const i64 b) {
    i64 x = a - b;
    if (x < 0) {
        x += MOD;
    }
    return x;
}

i64 mod_mul(const i64 a, const i64 b) {
    return static_cast<i64>((static_cast<__int128>(a) * static_cast<__int128>(b)) % MOD);
}

i64 mod_pow(i64 base, u64 exp) {
    base %= MOD;
    i64 result = 1;
    while (exp > 0) {
        if ((exp & 1ULL) != 0ULL) {
            result = mod_mul(result, base);
        }
        base = mod_mul(base, base);
        exp >>= 1ULL;
    }
    return result;
}

std::vector<int> build_spf(const int n) {
    std::vector<int> spf(static_cast<std::size_t>(n + 1));
    for (int i = 0; i <= n; ++i) {
        spf[static_cast<std::size_t>(i)] = i;
    }

    for (int i = 2; static_cast<i64>(i) * i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] != i) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            if (spf[static_cast<std::size_t>(j)] == j) {
                spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }

    return spf;
}

i64 solve(const u64 G, const u64 L, const u64 N) {
    const int K = static_cast<int>(L / G);
    const std::vector<int> spf = build_spf(K);

    int max_exp = 0;
    {
        int x = K;
        while (x > 0) {
            x /= 2;
            ++max_exp;
        }
    }

    std::vector<i64> pows(static_cast<std::size_t>(max_exp + 2), 0);
    for (int b = 0; b <= max_exp + 1; ++b) {
        if (b == 0) {
            pows[static_cast<std::size_t>(b)] = 0;  // N > 0 here
        } else {
            pows[static_cast<std::size_t>(b)] = mod_pow(b, N);
        }
    }

    std::vector<i64> local(static_cast<std::size_t>(max_exp + 1), 0);
    for (int a = 1; a <= max_exp; ++a) {
        // (a+1)^N - 2*a^N + (a-1)^N
        i64 value = mod_sub(pows[static_cast<std::size_t>(a + 1)], mod_mul(2, pows[static_cast<std::size_t>(a)]));
        value = mod_add(value, pows[static_cast<std::size_t>(a - 1)]);
        local[static_cast<std::size_t>(a)] = value;
    }

    std::vector<i64> h(static_cast<std::size_t>(K + 1), 0);
    h[1] = 1;

    for (int n = 2; n <= K; ++n) {
        int x = n;
        i64 value = 1;

        while (x > 1) {
            const int p = spf[static_cast<std::size_t>(x)];
            int e = 0;
            while (x % p == 0) {
                x /= p;
                ++e;
            }
            value = mod_mul(value, local[static_cast<std::size_t>(e)]);
        }

        h[static_cast<std::size_t>(n)] = value;
    }

    i64 answer = 0;
    for (int m = 1; m <= K; ++m) {
        const u64 weight = L / static_cast<u64>(m) - G + 1ULL;
        const i64 w = static_cast<i64>(weight % static_cast<u64>(MOD));
        answer = mod_add(answer, mod_mul(h[static_cast<std::size_t>(m)], w));
    }

    return answer;
}

bool run_checkpoints() {
    if (solve(10ULL, 100ULL, 1ULL) != 91) {
        std::cerr << "Checkpoint failed: f(10,100,1)\n";
        return false;
    }
    if (solve(10ULL, 100ULL, 2ULL) != 327) {
        std::cerr << "Checkpoint failed: f(10,100,2)\n";
        return false;
    }
    if (solve(10ULL, 100ULL, 3ULL) != 1135) {
        std::cerr << "Checkpoint failed: f(10,100,3)\n";
        return false;
    }
    if (solve(10ULL, 100ULL, 1000ULL) != 3286053) {
        std::cerr << "Checkpoint failed: f(10,100,1000) mod 101^4\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const i64 answer = solve(1000000ULL, 1000000000000ULL, 1000000000000000000ULL);
    std::cout << answer << '\n';
    return 0;
}

Python

MOD = 104060401

def mod_pow(base, exp):
    return pow(base, exp, MOD)

def build_spf(n):
    spf = list(range(n + 1))
    for i in range(2, int(n ** 0.5) + 1):
        if spf[i] == i:
            for j in range(i * i, n + 1, i):
                if spf[j] == j:
                    spf[j] = i
    return spf

def solve_350(G, L, N):
    K = L // G
    spf = build_spf(K)
    
    max_exp = 0
    x = K
    while x > 0:
        x //= 2
        max_exp += 1
        
    pows = [0] * (max_exp + 2)
    for b in range(max_exp + 2):
        if b == 0:
            pows[b] = 0
        else:
            pows[b] = mod_pow(b, N)
            
    local = [0] * (max_exp + 1)
    for a in range(1, max_exp + 1):
        val = (pows[a + 1] - 2 * pows[a] + pows[a - 1]) % MOD
        local[a] = (val + MOD) % MOD
        
    h = [0] * (K + 1)
    h[1] = 1
    
    for n in range(2, K + 1):
        x = n
        value = 1
        while x > 1:
            p = spf[x]
            e = 0
            while x % p == 0:
                x //= p
                e += 1
            value = (value * local[e]) % MOD
        h[n] = value
        
    answer = 0
    for m in range(1, K + 1):
        weight = L // m - G + 1
        w = weight % MOD
        answer = (answer + h[m] * w) % MOD
        
    return answer

def solve():
    G = 1000000
    L = 1000000000000
    N = 1000000000000000000
    ans = solve_350(G, L, N)
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;

public class Euler350 {

    static final long MOD = 104060401L;

    static long modPow(long base, long exp) {
        long result = 1;
        base %= MOD;
        while (exp > 0) {
            if ((exp & 1) == 1) {
                result = (result * base) % MOD;
            }
            base = (base * base) % MOD;
            exp >>= 1;
        }
        return result;
    }

    static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        for (int i = 0; i <= n; i++)
            spf[i] = i;

        for (int i = 2; i * i <= n; i++) {
            if (spf[i] == i) {
                for (int j = i * i; j <= n; j += i) {
                    if (spf[j] == j) {
                        spf[j] = i;
                    }
                }
            }
        }
        return spf;
    }

    static long solve(long G, long L, long N) {
        int K = (int) (L / G);
        int[] spf = buildSpf(K);

        int maxExp = 0;
        int x = K;
        while (x > 0) {
            x /= 2;
            maxExp++;
        }

        long[] pows = new long[maxExp + 2];
        for (int b = 0; b <= maxExp + 1; b++) {
            if (b == 0) {
                pows[b] = 0;
            } else {
                pows[b] = modPow(b, N);
            }
        }

        long[] local = new long[maxExp + 1];
        for (int a = 1; a <= maxExp; a++) {
            long value = (pows[a + 1] - 2 * pows[a] % MOD + MOD) % MOD;
            value = (value + pows[a - 1]) % MOD;
            local[a] = value;
        }

        long[] h = new long[K + 1];
        h[1] = 1;

        for (int n = 2; n <= K; n++) {
            int cx = n;
            long value = 1;

            while (cx > 1) {
                int p = spf[cx];
                int e = 0;
                while (cx % p == 0) {
                    cx /= p;
                    e++;
                }
                value = (value * local[e]) % MOD;
            }
            h[n] = value;
        }

        long answer = 0;
        for (int m = 1; m <= K; m++) {
            long weight = L / m - G + 1;
            long w = weight % MOD;
            answer = (answer + h[m] * w) % MOD;
        }

        return answer;
    }

    public static String solve() {
        long ans = solve(1000000L, 1000000000000L, 1000000000000000000L);
        return String.valueOf(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}