Problem 876: Triplet Tricks

Project Euler Problem 876 Solution

EulerSolve provides an optimized solution for Project Euler Problem 876, Triplet Tricks, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Starting from a triple $(a,b,c)$, one move reflects exactly one coordinate: $$ (a,b,c)\to (2(b+c)-a,b,c),\quad (a,b,c)\to (a,2(c+a)-b,c),\quad (a,b,c)\to (a,b,2(a+b)-c). $$ For fixed positive integers $a$ and $b$, the relevant positive values of $c$ are exactly those for which $$ (a+b-c)^2-4ab $$ is a perfect square. If $T(a,b,c)$ denotes the minimum number of moves needed to make at least one coordinate equal to $0$, then $$ F(a,b)=\sum T(a,b,c), $$ where the sum runs over all such positive $c$. The Project Euler task asks for $$ \sum_{k=1}^{18} F(2^k3^k,2^k5^k). $$ A direct graph search is useful only for very small checks. The full solution replaces that search by a divisor-pair parametrization and by counting reflection steps through the Euclidean algorithm. Mathematical Approach Fix $a$ and $b$. The goal is to describe all admissible third coordinates $c$, then evaluate the minimal distance to a zero coordinate for each of them. Step 1: Turn the condition on $c$ into a difference of squares Set $$ u=\lvert a+b-c\rvert. $$ The admissibility condition used by the implementations is that $$ (a+b-c)^2-4ab=s^2 $$ for some integer $s$. With $u=\lvert a+b-c\rvert$, this becomes $$ u^2-s^2=4ab. $$ So every admissible state is encoded by a factorization of $4ab$ as a difference of two squares....

Detailed mathematical approach

Problem Summary

Starting from a triple $(a,b,c)$, one move reflects exactly one coordinate:

$$ (a,b,c)\to (2(b+c)-a,b,c),\quad (a,b,c)\to (a,2(c+a)-b,c),\quad (a,b,c)\to (a,b,2(a+b)-c). $$

For fixed positive integers $a$ and $b$, the relevant positive values of $c$ are exactly those for which

$$ (a+b-c)^2-4ab $$

is a perfect square. If $T(a,b,c)$ denotes the minimum number of moves needed to make at least one coordinate equal to $0$, then

$$ F(a,b)=\sum T(a,b,c), $$

where the sum runs over all such positive $c$. The Project Euler task asks for

$$ \sum_{k=1}^{18} F(2^k3^k,2^k5^k). $$

A direct graph search is useful only for very small checks. The full solution replaces that search by a divisor-pair parametrization and by counting reflection steps through the Euclidean algorithm.

Mathematical Approach

Fix $a$ and $b$. The goal is to describe all admissible third coordinates $c$, then evaluate the minimal distance to a zero coordinate for each of them.

Step 1: Turn the condition on $c$ into a difference of squares

Set

$$ u=\lvert a+b-c\rvert. $$

The admissibility condition used by the implementations is that

$$ (a+b-c)^2-4ab=s^2 $$

for some integer $s$. With $u=\lvert a+b-c\rvert$, this becomes

$$ u^2-s^2=4ab. $$

So every admissible state is encoded by a factorization of $4ab$ as a difference of two squares.

Step 2: Parametrize admissible states by factor pairs

Write

$$ (u-s)(u+s)=4ab. $$

Define

$$ q=u-s,\qquad p=u+s. $$

Then $pq=4ab$, $p\ge q>0$, and $p+q=2u$ is even. Conversely, every factor pair

$$ pq=4ab,\qquad p\ge q,\qquad p+q\text{ even} $$

recovers

$$ u=\frac{p+q}{2},\qquad s=\frac{p-q}{2}. $$

This replaces the original search over $c$ by a finite search over divisor pairs of $4ab$.

Step 3: Each factor pair gives one or two positive values of $c$

Since $u=\lvert a+b-c\rvert$, the corresponding third coordinates are

$$ c_-=a+b+u,\qquad c_+=a+b-u. $$

The larger branch $c_-$ is always positive. The smaller branch $c_+$ is admissible only when

$$ a+b-u>0 \iff \frac{p+q}{2}<a+b \iff p+q<2(a+b). $$

If equality holds, then $c_+=0$, which is excluded because the sum defining $F(a,b)$ only uses positive third coordinates.

Step 4: The move count collapses to Euclidean quotient sums

For positive integers $x\ge y$, run the ordinary Euclidean algorithm

$$ x=q_0y+r_0,\qquad y=q_1r_0+r_1,\qquad \dots,\qquad r_{m-2}=q_m r_{m-1}. $$

Define the quotient sum

$$ E(x,y)=q_0+q_1+\cdots+q_m. $$

This is exactly the number of subtraction steps compressed by the division form of Euclid's algorithm. For a divisor pair $(p,q)$, the reflection system reduces to two possible Euclidean chains, one aiming to eliminate $a$ and one aiming to eliminate $b$. Their lengths are

$$ E(2a,q)\qquad\text{and}\qquad E(2b,q), $$

so the optimal distance for the larger branch is

$$ f(p,q)=\min\bigl(E(2a,q),E(2b,q)\bigr). $$

The two branches are consecutive along the same reflection chain, because reflecting the third coordinate sends

$$ 2(a+b)-c_-=2(a+b)-(a+b+u)=a+b-u=c_+. $$

Therefore $c_-$ contributes $f(p,q)$, while the smaller branch $c_+$ contributes one less move, namely $f(p,q)-1$, whenever $c_+$ is positive.

Step 5: Closed formula for $F(a,b)$

Summing over all factor pairs gives

$$ F(a,b)= \sum_{\substack{pq=4ab\\ q\le p\\ p+q\text{ even}}} \left( f(p,q)+ \begin{cases} f(p,q)-1, & p+q<2(a+b),\\ 0, & p+q\ge 2(a+b). \end{cases} \right). $$

This is the exact formula evaluated by the implementations.

Step 6: Specialize to the Project Euler input family

For the actual problem,

$$ a=2^k3^k,\qquad b=2^k5^k, $$

$$ 4ab=2^{2k+2}3^k5^k. $$

Hence every divisor $q$ of $4ab$ has the form

$$ q=2^{e_2}3^{e_3}5^{e_5}, \qquad 0\le e_2\le 2k+2, \qquad 0\le e_3,e_5\le k. $$

That is why the implementation can enumerate all candidates by three nested loops over exponents, then recover $p=(4ab)/q$ and apply the filters $q\le p$ and $p+q$ even.

Worked Example: $F(6,10)=17$

Take $k=1$. Then

$$ a=6,\qquad b=10,\qquad 4ab=240. $$

The valid factor pairs with $q\le p$ and $p+q$ even are

$$ (p,q)\in\{(120,2),(60,4),(40,6),(30,8),(24,10),(20,12)\}. $$

Now compute the Euclidean quotient sums:

$$ \begin{aligned} q=2&:& E(12,2)=6,\quad E(20,2)=10,\quad f=6,\\ q=4&:& E(12,4)=3,\quad E(20,4)=5,\quad f=3,\\ q=6&:& E(12,6)=2,\quad E(20,6)=6,\quad f=2,\\ q=8&:& E(12,8)=3,\quad E(20,8)=4,\quad f=3,\\ q=10&:& E(12,10)=6,\quad E(20,10)=2,\quad f=2,\\ q=12&:& E(12,12)=1,\quad E(20,12)=4,\quad f=1. \end{aligned} $$

Here $2(a+b)=32$, and every valid pair has $p+q\ge 32$, so the smaller branch never contributes a positive $c$. Therefore

$$ F(6,10)=6+3+2+3+2+1=17, $$

which matches the sample check used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same arithmetic plan. First they precompute the powers of $2$, $3$, and $5$ needed to build $a$, $b$, and $4ab$ for every $1\le k\le 18$. For each $k$, they enumerate all divisors $q$ of $4ab$ through the exponent ranges above, set $p=(4ab)/q$, and discard pairs with $q>p$ or odd $p+q$.

For each surviving divisor pair they evaluate two Euclidean quotient sums, one for $(2a,q)$ and one for $(2b,q)$, keep the smaller of the two, and add the two branch contributions exactly as in the closed formula. The C++ implementation parallelizes independent values of $k$; the Python and Java implementations perform the same computation sequentially.

Small internal checks confirm the derived formula on tiny cases and on the checkpoints

$$ F(6,10)=17,\qquad F(36,100)=179. $$

Complexity Analysis

For a fixed $k$, the raw divisor enumeration uses

$$ (2k+3)(k+1)^2 $$

choices of $(e_2,e_3,e_5)$. Each surviving choice performs two Euclidean algorithms, and each Euclidean computation takes $O(\log(4ab))$ division steps. So the arithmetic cost for one $k$ is

$$ O\bigl((2k+3)(k+1)^2\log(4ab)\bigr). $$

If one generalizes the upper limit from $18$ to $K$, then the total candidate count is

$$ \sum_{k=1}^{K}(2k+3)(k+1)^2=O(K^4). $$

The memory usage is $O(1)$ beyond the small power tables and loop variables.

Footnotes and References

Problem page: https://projecteuler.net/problem=876
Euclidean algorithm: Wikipedia — Euclidean algorithm
Difference of two squares: Wikipedia — Difference of two squares
Divisor function: Wikipedia — Divisor function

Problem 876 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdlib>
#include <cstdint>
#include <deque>
#include <iostream>
#include <thread>
#include <unordered_set>
#include <utility>
#include <vector>

using namespace std;

namespace {

using u128 = unsigned __int128;

u128 euclid_sum(u128 a, u128 b) {
    if (a < b) swap(a, b);
    u128 sum = 0;
    while (b != 0) {
        sum += a / b;
        u128 r = a % b;
        a = b;
        b = r;
    }
    return sum;
}

string to_string_u128(u128 v) {
    if (v == 0) return "0";
    string s;
    while (v > 0) {
        int digit = static_cast<int>(v % 10);
        s.push_back(static_cast<char>('0' + digit));
        v /= 10;
    }
    reverse(s.begin(), s.end());
    return s;
}

struct Triple {
    int a;
    int b;
    int c;
};

struct TripleHash {
    size_t operator()(const Triple& t) const noexcept {
        size_t h1 = std::hash<int>{}(t.a);
        size_t h2 = std::hash<int>{}(t.b);
        size_t h3 = std::hash<int>{}(t.c);
        size_t h = h1;
        h ^= h2 + 0x9e3779b9 + (h << 6) + (h >> 2);
        h ^= h3 + 0x9e3779b9 + (h << 6) + (h >> 2);
        return h;
    }
};

bool operator==(const Triple& x, const Triple& y) {
    return x.a == y.a && x.b == y.b && x.c == y.c;
}

int min_steps_bfs(int a, int b, int c, int bound, int depth_limit) {
    if (a == 0 || b == 0 || c == 0) return 0;
    deque<pair<Triple, int>> q;
    unordered_set<Triple, TripleHash> seen;

    q.push_back({{a, b, c}, 0});
    seen.insert({a, b, c});

    while (!q.empty()) {
        auto [cur, dist] = q.front();
        q.pop_front();
        if (dist >= depth_limit) continue;

        const int aa = cur.a;
        const int bb = cur.b;
        const int cc = cur.c;

        const Triple nexts[3] = {
            {2 * (bb + cc) - aa, bb, cc},
            {aa, 2 * (cc + aa) - bb, cc},
            {aa, bb, 2 * (aa + bb) - cc}
        };

        for (const auto& nxt : nexts) {
            if (nxt.a == 0 || nxt.b == 0 || nxt.c == 0) return dist + 1;
            if (max({abs(nxt.a), abs(nxt.b), abs(nxt.c)}) > bound) continue;
            if (seen.insert(nxt).second) {
                q.push_back({nxt, dist + 1});
            }
        }
    }
    return 0;
}

u128 compute_F_bruteforce(int a, int b) {
    const int N = 4 * a * b;
    unordered_set<int> candidates;
    for (int d = 1; d * d <= N; ++d) {
        if (N % d != 0) continue;
        const int p = d;
        const int q = N / d;
        if ((p + q) & 1) continue;
        const int u = (p + q) / 2;
        const int c_pos = a + b - u;
        if (c_pos > 0) candidates.insert(c_pos);
        const int c_neg = a + b + u;
        if (c_neg > 0) candidates.insert(c_neg);
    }

    u128 total = 0;
    for (int c : candidates) {
        int f = min_steps_bfs(a, b, c, 500, 12);
        total += static_cast<u128>(f);
    }
    return total;
}

u128 compute_F_power(int k,
                     const vector<u128>& pow2,
                     const vector<u128>& pow3,
                     const vector<u128>& pow5) {
    const u128 a = pow2[k] * pow3[k];
    const u128 b = pow2[k] * pow5[k];
    const u128 two_a = a * 2;
    const u128 two_b = b * 2;
    const u128 two_ab = (a + b) * 2;

    const u128 N = pow2[2 * k + 2] * pow3[k] * pow5[k];

    u128 total = 0;

    for (int e2 = 0; e2 <= 2 * k + 2; ++e2) {
        const u128 v2 = pow2[e2];
        for (int e3 = 0; e3 <= k; ++e3) {
            const u128 v23 = v2 * pow3[e3];
            for (int e5 = 0; e5 <= k; ++e5) {
                const u128 q = v23 * pow5[e5];
                const u128 p = N / q;
                if (q > p) continue;
                if ((p + q) & 1) continue;

                // For each divisor pair p*q=4ab, q is the smaller factor.
                // The step count to zero a or b equals the sum of Euclidean
                // quotients for (2a,q) or (2b,q); take the minimum. For the
                // positive-u branch (c <= a+b), we need one fewer step.
                const u128 s1 = euclid_sum(max(two_a, q), min(two_a, q));
                const u128 s2 = euclid_sum(max(two_b, q), min(two_b, q));
                const u128 f = min(s1, s2);

                total += f; // u negative branch is always valid.
                if (p + q < two_ab) {
                    total += f - 1; // u positive branch, if c>0.
                }
            }
        }
    }
    return total;
}

void run_validations(const vector<u128>& pow2,
                     const vector<u128>& pow3,
                     const vector<u128>& pow5) {
    const u128 f1 = compute_F_power(1, pow2, pow3, pow5);
    const u128 f2 = compute_F_power(2, pow2, pow3, pow5);
    if (f1 != 17 || f2 != 179) {
        cerr << "Validation failed for given samples: F(6,10)="
             << to_string_u128(f1) << ", F(36,100)="
             << to_string_u128(f2) << "\n";
        exit(1);
    }

    const vector<pair<int, int>> small_cases = {
        {2, 3},
        {3, 5},
        {4, 6}
    };

    for (const auto& [a, b] : small_cases) {
        const u128 brute = compute_F_bruteforce(a, b);
        // Compute formula result via brute divisor enumeration (small N).
        const int N = 4 * a * b;
        u128 formula = 0;
        for (int d = 1; d * d <= N; ++d) {
            if (N % d != 0) continue;
            const int q = d;
            const int p = N / d;
            if ((p + q) & 1) continue;
            const u128 s1 = euclid_sum(max<u128>(2ULL * a, q), min<u128>(2ULL * a, q));
            const u128 s2 = euclid_sum(max<u128>(2ULL * b, q), min<u128>(2ULL * b, q));
            const u128 f = min(s1, s2);
            formula += f;
            if (p + q < 2 * (a + b)) {
                formula += f - 1;
            }
        }
        if (brute != formula) {
            cerr << "Validation failed for small case (" << a << "," << b
                 << "): brute=" << to_string_u128(brute)
                 << ", formula=" << to_string_u128(formula) << "\n";
            exit(1);
        }
    }
}

} // namespace

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int maxK = 18;
    unsigned threads = thread::hardware_concurrency();
    bool validate = true;

    // Optional CLI: ./a.out [maxK] [threads] [validate(0/1)]
    if (argc >= 2) maxK = stoi(argv[1]);
    if (argc >= 3) threads = static_cast<unsigned>(stoul(argv[2]));
    if (argc >= 4) validate = (stoi(argv[3]) != 0);

    if (maxK <= 0) {
        cout << "0\n";
        return 0;
    }

    const int max_pow2 = 2 * maxK + 2;
    vector<u128> pow2(max_pow2 + 1, 1);
    for (int i = 1; i <= max_pow2; ++i) pow2[i] = pow2[i - 1] * 2;

    vector<u128> pow3(maxK + 1, 1);
    vector<u128> pow5(maxK + 1, 1);
    for (int i = 1; i <= maxK; ++i) {
        pow3[i] = pow3[i - 1] * 3;
        pow5[i] = pow5[i - 1] * 5;
    }

    if (validate) {
        run_validations(pow2, pow3, pow5);
    }

    if (threads == 0) threads = 1;
    threads = min<unsigned>(threads, static_cast<unsigned>(maxK));

    vector<u128> results(maxK + 1, 0);
    atomic<int> nextK(1);

    auto worker = [&]() {
        while (true) {
            int k = nextK.fetch_add(1);
            if (k > maxK) break;
            results[k] = compute_F_power(k, pow2, pow3, pow5);
        }
    };

    vector<thread> pool;
    pool.reserve(threads);
    for (unsigned i = 0; i < threads; ++i) {
        pool.emplace_back(worker);
    }
    for (auto& th : pool) th.join();

    u128 total = 0;
    for (int k = 1; k <= maxK; ++k) {
        total += results[k];
    }

    cout << to_string_u128(total) << "\n";
    return 0;
}

Python

import sys

def euclid_sum(a, b):
    if a < b:
        a, b = b, a
    total = 0
    while b != 0:
        total += a // b
        r = a % b
        a = b
        b = r
    return total

def compute_F_power(k, pow2, pow3, pow5):
    a = pow2[k] * pow3[k]
    b = pow2[k] * pow5[k]
    two_a = a * 2
    two_b = b * 2
    two_ab = (a + b) * 2

    N = pow2[2 * k + 2] * pow3[k] * pow5[k]

    total = 0

    for e2 in range(2 * k + 3):
        v2 = pow2[e2]
        for e3 in range(k + 1):
            v23 = v2 * pow3[e3]
            for e5 in range(k + 1):
                q = v23 * pow5[e5]
                p = N // q
                if q > p:
                    continue
                if (p + q) & 1:
                    continue

                s1 = euclid_sum(max(two_a, q), min(two_a, q))
                s2 = euclid_sum(max(two_b, q), min(two_b, q))
                f = min(s1, s2)

                total += f
                if p + q < two_ab:
                    total += f - 1

    return total

def solve():
    maxK = 18
    max_pow2 = 2 * maxK + 2
    
    pow2 = [1] * (max_pow2 + 1)
    for i in range(1, max_pow2 + 1):
        pow2[i] = pow2[i - 1] * 2

    pow3 = [1] * (maxK + 1)
    pow5 = [1] * (maxK + 1)
    for i in range(1, maxK + 1):
        pow3[i] = pow3[i - 1] * 3
        pow5[i] = pow5[i - 1] * 5

    total = 0
    for k in range(1, maxK + 1):
        total += compute_F_power(k, pow2, pow3, pow5)

    return str(total)

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

Java

import java.math.BigInteger;

public class Euler876 {

    static BigInteger euclidSum(BigInteger a, BigInteger b) {
        if (a.compareTo(b) < 0) {
            BigInteger temp = a;
            a = b;
            b = temp;
        }
        BigInteger sum = BigInteger.ZERO;
        while (b.compareTo(BigInteger.ZERO) > 0) {
            BigInteger[] divAndRem = a.divideAndRemainder(b);
            sum = sum.add(divAndRem[0]);
            a = b;
            b = divAndRem[1];
        }
        return sum;
    }

    static BigInteger computeFPower(int k, BigInteger[] pow2, BigInteger[] pow3, BigInteger[] pow5) {
        BigInteger a = pow2[k].multiply(pow3[k]);
        BigInteger b = pow2[k].multiply(pow5[k]);
        BigInteger twoA = a.shiftLeft(1);
        BigInteger twoB = b.shiftLeft(1);
        BigInteger twoAb = a.add(b).shiftLeft(1);

        BigInteger N = pow2[2 * k + 2].multiply(pow3[k]).multiply(pow5[k]);
        BigInteger total = BigInteger.ZERO;

        for (int e2 = 0; e2 <= 2 * k + 2; ++e2) {
            BigInteger v2 = pow2[e2];
            for (int e3 = 0; e3 <= k; ++e3) {
                BigInteger v23 = v2.multiply(pow3[e3]);
                for (int e5 = 0; e5 <= k; ++e5) {
                    BigInteger q = v23.multiply(pow5[e5]);
                    BigInteger p = N.divide(q);

                    if (q.compareTo(p) > 0)
                        continue;
                    if (p.add(q).testBit(0))
                        continue;

                    BigInteger s1 = euclidSum(twoA.max(q), twoA.min(q));
                    BigInteger s2 = euclidSum(twoB.max(q), twoB.min(q));
                    BigInteger f = s1.min(s2);

                    total = total.add(f);
                    if (p.add(q).compareTo(twoAb) < 0) {
                        total = total.add(f.subtract(BigInteger.ONE));
                    }
                }
            }
        }
        return total;
    }

    public static String solve() {
        int maxK = 18;
        int maxPow2 = 2 * maxK + 2;

        BigInteger[] pow2 = new BigInteger[maxPow2 + 1];
        pow2[0] = BigInteger.ONE;
        for (int i = 1; i <= maxPow2; ++i) {
            pow2[i] = pow2[i - 1].shiftLeft(1);
        }

        BigInteger[] pow3 = new BigInteger[maxK + 1];
        BigInteger[] pow5 = new BigInteger[maxK + 1];
        pow3[0] = BigInteger.ONE;
        pow5[0] = BigInteger.ONE;

        BigInteger three = BigInteger.valueOf(3);
        BigInteger five = BigInteger.valueOf(5);

        for (int i = 1; i <= maxK; ++i) {
            pow3[i] = pow3[i - 1].multiply(three);
            pow5[i] = pow5[i - 1].multiply(five);
        }

        BigInteger total = BigInteger.ZERO;
        for (int k = 1; k <= maxK; ++k) {
            total = total.add(computeFPower(k, pow2, pow3, pow5));
        }

        return total.toString();
    }

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