Untitled

#include <bits/stdc++.h>

using namespace std;

using ll = long long;

vector<int> primes = {2, 3, 5, 7, 11};

int modInv(int n, int p) {
    // Return the modular inverse of n mod p.
    //
    // 0 ~= p = q*n + r hence r ~= -q*n hence 1/n ~= -q/r
    int a = 1;
    while (1 < n) {
        // The answer is a*1/n
        a *= -p/n, n = p%n;
    }
    return a;
}

int mod2Equations(int a, int ra, int p, int rp) {
    // Return the smaller nonnegative solution n of the system of two
    // modular equations (n ~= ra [mod a] and n ~= rp [mod p]), where
    // p is a prime number and a is a number coprime with p.

    // n ~= ra [mod a] is equivalent to n = qa*a + ra for some qa.
    // Hence, n ~= rp [mod p] rewrites as qa*a ~= rp - ra [mod p].
    // Since p is prime, this is equivalent to qa ~= (rp - ra)/a [mod
    // p] where 1/a is the modular inverse of a.
    int qa = ((rp - ra)%p + p)%p*modInv(a%p, p)%p;
    return qa*a + ra;
}

int modSystem(const vector<int> &remainders) {
    // Return the smaller nonnegative solution of the system of
    // modular equations (n ~= remainders[i] [mod primes[i]] for each
    // i).
    int n = remainders[0];
    int a = primes[0];
    for (int i = 1; i < 5; i++) {
        n = mod2Equations(a, n, primes[i], remainders[i]), a *= primes[i];
    }
    return n;
}

bool dfs(int i, vector<vector<int>> &remaindersABC,
         const vector<vector<vector<int>>> &possibilities,
         vector<int> &ans) {
    if (i == 5) {
        int a = modSystem(remaindersABC[0]);
        int b = modSystem(remaindersABC[1]);
        int c = modSystem(remaindersABC[2]);
        if (0 < a && a < b && b < c && a*a + b*b == c*c && a + b + c == 1000) {
            // cout << a*a << " + " << b*b << " == " << c*c << "\n";
            ans = {a, b, c};
            return true;
        }
        else {
            // cout << a*a << " + " << b*b << " != "  << c*c << "\n";
            return false;
        }
    }
    else {
        int m = possibilities[i].size();
        for (int j = 0; j < m; j++) {
            // Try the possibility j for the remainders modulo the
            // prime number indexed i
            for (int x = 0; x < 3; x++) {
                remaindersABC[x][i] = possibilities[i][j][x];
            }
            if (dfs(i + 1, remaindersABC, possibilities, ans))
                return true;
        }
        return false;
    }
}

int main() {
    // We search the unique triplet of positive numbers a < b < c such
    // that (1) a^2 + b^2 = c^2 and (2) a + b + c = 1000.
    // 
    // Each number up to 1000 < 2*3*5*7*11 is uniquely defined by its
    // remainders modulo the prime numbers up to 11 as a special case
    // of the Chinese remainders theorem.
    //
    // We enumerate over each possible remainders for a, b and c
    // compatible with equations (1) and (2), and to find the
    // corresponding numbers actually satisfying these equations.
    
    vector<vector<vector<int>>> possibilities(5);
    // possibilities[i]: array with each possible triplets of values
    // {ra, rb, rc} for the remainders of a, b and c mod primes[i],
    // according to equations (1) and (2)
    for (int i = 0; i < 5; i++) {
        // compute possibilities[i]
        for (int ra = 0; ra < primes[i]; ra++) {
            for (int rb = 0; rb < primes[i]; rb++) {
                // ra, rb <= 11 so ra + rb <= 1000
                int rc = (1000 - ra - rb)%primes[i];
                if ((ra*ra + rb*rb)%primes[i] == rc*rc%primes[i])
                    possibilities[i].push_back({ra, rb, rc});
            }
        }
    }

    vector<vector<int>> remaindersABC(3, vector<int>(5));
    // stores the remainders of a, b and c selected during the dfs

    vector<int> ans;
    bool res = dfs(0, remaindersABC, possibilities, ans);
    cout << "a = " << ans[0] << "\nb = " << ans[1] << "\nc = " << ans[2] << "\n";
    cout << "a*b*c = " << ans[0]*ans[1]*ans[2] << "\n";
}