Detailed mathematical approach

Problem Summary

For positive integers \(a,b\), one step of the modified Euclidean algorithm divides the larger number by the smaller one and replaces the larger number by the integer quotient, not by the remainder. The process stops as soon as one entry equals \(1\), and \(f(a,b)\) is the other entry. We must evaluate

$$E(N)=\sum_{1\le a,b\lt N}f(a,b),\qquad N=3\,000\,000.$$

The statement gives \(f(123,456)=3\), \(E(10)=343\) and \(E(100)=269288\). ^[1]

Mathematical Approach

1. Basic properties

The classical algorithm keeps the remainder; this variant keeps the quotient. ^[2] A step maps \((a,b)\) to \((a,\lfloor b/a\rfloor)\) when \(a\le b\), and to \((\lfloor a/b\rfloor,b)\) when \(a\gt b\). Four facts follow directly from this rule.

• Termination. While both entries are at least \(2\), the divided entry at least halves, because \(\lfloor b/a\rfloor\le b/2\) for \(a\ge2\). Every quotient is at least \(1\), so no entry reaches \(0\). The product of the two entries at least halves at each step, so at most \(\log_2(ab)\) steps occur.
• Symmetry. Swapping the two starting entries swaps them in every later pair as well, so \(f(a,b)=f(b,a)\). When \(a=b\), either entry may be divided; both choices give a pair made of \(1\) and \(a\).
• Boundary values. If an entry is already \(1\), nothing is divided: \(f(1,b)=b\) and \(f(a,1)=a\). For \(a\ge2\), one step turns \((a,a)\) into \((a,1)\), so \(f(a,a)=a\).
• No growth. Entries never increase, since a quotient never exceeds the number divided. Hence \(f(a,b)\le\max(a,b)\).

For example, \((123,456)\to(123,3)\to(41,3)\to(13,3)\to(4,3)\to(1,3)\), so \(f(123,456)=3\).

2. The first division

Let \(T=\sum_{c=1}^{N-1}c=N(N-1)/2\). The pairs containing a \(1\) contribute \(2T-1\): the row \(a=1\) and the column \(b=1\) each sum to \(T\), and \((1,1)\), whose value is \(1\), was counted twice. The diagonal \(2\le a=b\lt N\) contributes \(T-1\). By symmetry, the remaining pairs split into two equal halves, and we keep the half with \(a\lt b\):

$$E(N)=3T-2+2\sum_{2\le a\lt b\lt N}f(a,b).$$

For \(2\le a\lt b\), the first step replaces \(b\) by \(q=\lfloor b/a\rfloor\ge1\), so \(f(a,b)=f(a,q)\). Two cases arise.

• If \(q=1\), that is \(a\lt b\lt 2a\), the pair becomes \((a,1)\) and \(f(a,b)=a\). For fixed \(a\) there are \(\min(2a-1,N-1)-a\) such values of \(b\).
• If \(q\ge2\), division with remainder shows that exactly the integers \(b\) with \(aq\le b\le aq+a-1\) have quotient \(q\). ^[3] Intersecting this block with \(b\le N-1\) leaves \(\min(a,N-aq)\) values whenever \(aq\lt N\), and all of them exceed \(a\).

Therefore

$$E(N)=3T-2+2\bigl(Q_1+S\bigr),\qquad Q_1=\sum_{a=2}^{N-1}a\bigl(\min(2a-1,N-1)-a\bigr),$$

$$S=\sum_{\substack{a,q\ge2\\ aq\lt N}}\min(a,N-aq)\,f(a,q).$$

The quotient-one sum \(Q_1\) costs one pass over \(a\). The whole difficulty now lies in \(S\), which has about \(N\ln N\) terms.

3. A self-similar table for the smaller entry

In \(S\), put \(m=\min(a,q)\) and \(x=\max(a,q)\). Then \(2\le m\le x\) and \(mx\lt N\), so \(m^2\le N-1\) and \(m\le r=\lfloor\sqrt{N-1}\rfloor\). For fixed \(m\), define \(v_m(x)=f(x,m)\). If \(x\gt m\), the first step divides \(x\) by \(m\). If \(x=m\), the value is \(m=f(1,m)\). In both cases

$$v_m(x)=v_m\bigl(\lfloor x/m\rfloor\bigr)\quad(x\ge m),\qquad v_m(1)=m.$$

The values \(v_m(x)\) with \(2\le x\lt m\) are obtained by running the algorithm. All indices with the same quotient \(k=\lfloor x/m\rfloor\) form the block \(km\le x\le km+m-1\), and the whole block copies the single earlier value \(v_m(k)\). Since \(k\le x/2\), that value is already known when the block is reached.

A pair \(\{m,x\}\) with \(x\gt m\) comes from the two ordered pairs \((a,q)=(m,x)\) and \((a,q)=(x,m)\), with multiplicities \(\min(m,N-mx)\) and \(\min(x,N-mx)\). By symmetry both have the value \(v_m(x)\). For \(x=m\) there is only one ordered pair. Hence

$$S=\sum_{m=2}^{r}\ \sum_{x=m}^{\lfloor (N-1)/m\rfloor}w_m(x)\,v_m(x),\qquad w_m(x)=\min(m,N-mx)+[x\ne m]\,\min(x,N-mx).$$

Here \([x\ne m]\) equals \(1\) when \(x\ne m\) and \(0\) otherwise. Every term of \(S\) is now read from a table that is filled by copying, not by running the algorithm again.

4. Worked example: \(E(10)\)

For \(N=10\) we have \(T=45\), so \(3T-2=133\). The quotient-one sum is

$$Q_1=2\cdot1+3\cdot2+4\cdot3+5\cdot4+6\cdot3+7\cdot2+8\cdot1+9\cdot0=80.$$

Here \(r=3\). For \(m=2\), the indices \(x=2,3,4\) all have \(v_2(x)=2\), with weights \(w_2(2)=2\), \(w_2(3)=2+3=5\) and \(w_2(4)=2+2=4\), so they contribute \(2\cdot11=22\). For \(m=3\), only \(x=3\) occurs, with \(w_3(3)=\min(3,1)=1\) and \(v_3(3)=3\). Thus \(S=25\) and

$$E(10)=133+2(80+25)=343.$$

How the Code Works

modified_euclid simulates the definition directly. It fills the short range \(2\le x\lt m\) of every table and drives the brute-force checks.

solve first computes \(T\) and \(Q_1\) in one pass and finds \(r\) as the largest integer with \(r^2\lt N\). Worker threads then share the values \(m=2,\ldots,r\). Each worker claims the next \(m\) with an atomic fetch_add, so every \(m\) is processed exactly once. ^[4] A small \(m\) carries far more work, roughly \(N/m\) table entries, so this dynamic assignment balances the threads better than fixed ranges would.

Each worker owns one table value and enlarges it only when the current range \(\lfloor(N-1)/m\rfloor\) does not fit. A worker receives increasing values of \(m\), so its first allocation is its largest. For each \(m\), the worker sets value[1] to \(m\), fills \(2\le x\lt m\) directly, and sweeps the blocks \(k=1,2,\ldots\). In block \(k\) it copies value[k] into every index \(x\) with \(\lfloor x/m\rfloor=k\), up to \(\lfloor(N-1)/m\rfloor\), and adds \(w_m(x)\,v_m(x)\) to its sum.

A crude bound shows why wide integers are used. From \(f(a,b)\le N-1\) we only know \(E(N)\le(N-1)^3\), and \((N-1)^3\) exceeds \(2^{64}\) for the target \(N\). The C++ therefore keeps its partial sums in unsigned 128-bit integers, a compiler extension supported by GCC and Clang. ^[5] Each single product \(w_m(x)\,v_m(x)\) is below \(2N^2\) and fits in 64 bits. The result is printed by a short decimal conversion routine, and the elapsed time goes to standard error.

Complexity and Verification

Computing \(T\) and \(Q_1\) takes \(O(N)\) time. The direct part runs the algorithm on \(m-2\) pairs for each \(m\), which is \((r-1)(r-2)/2=O(N)\) runs of \(O(\log N)\) steps. The block sweep makes \(\lfloor(N-1)/m\rfloor-m+1\) updates for each \(m\), so its total is

$$\sum_{m=2}^{r}\Bigl(\Bigl\lfloor\frac{N-1}{m}\Bigr\rfloor-m+1\Bigr)=\tfrac12N\ln N+\Bigl(\gamma-\tfrac32\Bigr)N+O\bigl(\sqrt N\bigr),$$

because the harmonic numbers satisfy \(H_r=\ln r+\gamma+O(1/r)\) and \(r^2=N+O(\sqrt N)\). ^[6] The running time is therefore \(O(N\log N)\). For \(N=3\,000\,000\) we get \(r=1732\): the sweep makes \(19\,603\,698\) updates and the direct part \(1\,497\,315\) runs. The largest table has \(\lfloor(N-1)/2\rfloor+1=1\,500\,000\) entries of 32 bits, about 6 MB. The result does not depend on how the values of \(m\) are distributed among the threads, because every \(m\) is processed once and integer addition is exact.

Before the main computation, the program runs these exact checks: ^[7]

• \(f(123,456)=f(456,123)=3\).
• \(E(10)=343\) and \(E(100)=269288\), each by brute force and by the fast method.
• The single-threaded fast method equals a brute-force double loop for every \(2\le N\le300\).
• The multithreaded fast method equals brute force for \(N=1000\), \(2023\) and \(4096\).
• One thread and the full thread count give the same value for \(N=200003\).

A separate full-size computation checked the regrouping independently. It evaluated \(S\) directly over all \(39\,205\,665\) pairs \((a,q)\) with \(a,q\ge2\) and \(aq\lt N\), running the modified algorithm on each pair, without the grouping by \(m\) and without the block tables. It agreed with the fast method at \(N=3\,000\,000\), as did the Python and Java ports. This comparison is additional validation, separate from the packaged checks.

Footnotes and References

The references below supply the background facts used above. The reduction to weighted quotient pairs, the self-similar tables and the operation counts are derived explicitly in this article.

1. Project Euler 1011 — Modified Euclidean Algorithm. The official statement defines the quotient-replacing step and gives the checkpoints for f(123, 456), E(10) and E(100). These examples test the implementation without disclosing the requested final value.
2. Euclid's Elements, Book VII, Proposition 2. D. E. Joyce's online edition, Clark University. Euclid finds the greatest common measure by repeatedly subtracting the smaller number from the larger, which amounts to keeping the remainder. The variant in this problem keeps the quotient instead; Section 1 shows why it still terminates.
3. Concrete Mathematics. R. L. Graham, D. E. Knuth and O. Patashnik, Concrete Mathematics, second edition, Chapter 3, §3.1 (floors and ceilings) and §3.4 (the binary operation mod). Writing b = a⌊b/a⌋ + (b mod a) with 0 ≤ b mod a < a identifies the block of integers that share one quotient, as used in Sections 2 and 3.
4. cppreference — std::atomic<T>::fetch_add. The operation atomically adds to the stored value and returns the previous value. Each worker uses it to claim the next m, so no value is skipped or processed twice.
5. GCC manual — 128-bit Integers. The unsigned __int128 extension holds the partial sums, since the a priori bound (N−1)³ does not fit in 64 bits.
6. NIST Digital Library of Mathematical Functions — §5.4, equation 5.4.14. It expresses the harmonic number H_n through the digamma function as ψ(n+1) + γ, and the asymptotic expansion 5.11.2 of ψ then gives H_n = ln n + γ + O(1/n). This yields the operation count in the complexity section.
7. C++ — Euler1011.cpp. The linked immutable C++ revision contains the weighted quotient-pair sum, the worker threads and the checks that run before the main computation. It is the implementation record for this article. The separate full-size comparison described above is additional validation, not a test packaged in that source file.