Problem 557: Cutting Triangles

Problem 557 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <set>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

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

// For a valid cut, with integer piece areas (a,b,c,d) and b<=c, one can solve for the
// (generally rational) areas X,Y of triangles BCP and CAP in terms of (a,b,c). The
// total area is:
//   T = a*(a+b)*(a+c) / (a^2 - b*c)
// and the cut exists with integer (a,b,c,d) iff (a^2 - b*c) > 0 and T is an integer
// (then d = T-a-b-c is automatically a positive integer).
static u128 S_range(int n, int a_begin, int a_end) {
    u128 total = 0;
    for (int a = a_begin; a <= a_end; ++a) {
        const int a_room = n - a;
        if (a_room < 3) continue;  // need b,c,d >= 1

        // Since b<=c and d>=1: a + 2b + 1 <= n.
        int bmax = std::min(a - 1, (a_room - 1) / 2);
        if (bmax < 1) continue;

        // T is increasing in c. The minimum for given (a,b) happens at c=b:
        //   T_min = a*(a+b)/(a-b).
        // If T_min > n then no c works, which gives an additional b upper bound.
        const u64 bound_num = static_cast<u64>(a) * static_cast<u64>(n - a);
        const u64 bound_den = static_cast<u64>(a + n);
        const int bmax_by_Tmin = static_cast<int>(bound_num / bound_den);
        bmax = std::min(bmax, bmax_by_Tmin);
        if (bmax < 1) continue;

        const u64 a2 = static_cast<u64>(a) * static_cast<u64>(a);

        for (int b = 1; b <= bmax; ++b) {
            const int rem = n - a - b;
            if (rem <= 2) continue;  // need c>=b and d>=1

            // From T<=n:
            //   c * (a(a+b) + n b) <= a^2 (n-a-b)
            const u64 rhs = a2 * static_cast<u64>(rem);
            const u64 lhs_coeff = static_cast<u64>(a) * static_cast<u64>(a + b) + static_cast<u64>(n) * static_cast<u64>(b);
            const u64 cmax_T = rhs / lhs_coeff;

            // From bc < a^2 and c integer:
            const u64 cmax_bc = (a2 - 1) / static_cast<u64>(b);

            // From a+b+c+d<=n and d>=1:
            const u64 cmax_sum = static_cast<u64>(n - a - b - 1);

            u64 cmax = std::min({cmax_T, cmax_bc, cmax_sum});
            if (cmax < static_cast<u64>(b)) continue;

            const u64 base = static_cast<u64>(a) * static_cast<u64>(a + b);  // fits in u64 for n<=10000
            for (int c = b; c <= static_cast<int>(cmax); ++c) {
                const u64 bc = static_cast<u64>(b) * static_cast<u64>(c);
                const u64 denom = a2 - bc;  // >0 by construction

                const u64 numer = base * static_cast<u64>(a + c);
                if (numer % denom != 0) continue;

                const u64 T = numer / denom;
                if (T > static_cast<u64>(n)) continue;

                const int d = static_cast<int>(T) - a - b - c;
                if (d <= 0) continue;

                total += T;
            }
        }
    }

    return total;
}

static u128 S(int n) {
    unsigned threads = std::max(1U, std::thread::hardware_concurrency());
    if (threads <= 1U || n < 200) {
        return S_range(n, 1, n);
    }

    threads = std::min<unsigned>(threads, static_cast<unsigned>(n));
    const int chunk = (n + static_cast<int>(threads) - 1) / static_cast<int>(threads);

    std::vector<u128> partial(static_cast<std::size_t>(threads), 0U);
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(threads));

    for (unsigned t = 0; t < threads; ++t) {
        const int begin = static_cast<int>(t) * chunk + 1;
        const int end = std::min(n, begin + chunk - 1);
        workers.emplace_back([&, t, begin, end]() {
            if (begin > end) {
                return;
            }
            partial[static_cast<std::size_t>(t)] = S_range(n, begin, end);
        });
    }
    for (auto& th : workers) {
        th.join();
    }

    u128 total = 0U;
    for (u128 v : partial) {
        total += v;
    }
    return total;
}

static std::vector<std::array<int, 4>> quads_with_total_area(int target_total) {
    std::vector<std::array<int, 4>> quads;
    for (int a = 1; a <= target_total - 3; ++a) {
        for (int b = 1; b <= std::min(a - 1, target_total - a - 2); ++b) {
            for (int c = b; c <= target_total - a - b - 1; ++c) {
                if (1LL * b * c >= 1LL * a * a) break;
                const u64 a2 = static_cast<u64>(a) * static_cast<u64>(a);
                const u64 denom = a2 - static_cast<u64>(b) * static_cast<u64>(c);
                const u128 numer = static_cast<u128>(a) * static_cast<u128>(a + b) * static_cast<u128>(a + c);
                if (numer % denom != 0) continue;
                const u64 T = static_cast<u64>(numer / denom);
                if (static_cast<int>(T) != target_total) continue;
                const int d = target_total - a - b - c;
                if (d <= 0) continue;
                quads.push_back({a, b, c, d});
            }
        }
    }
    std::sort(quads.begin(), quads.end());
    return quads;
}

}  // namespace

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    assert(to_string_u128(S(20)) == "259");

    {
        const auto quads = quads_with_total_area(55);
        const std::vector<std::array<int, 4>> expected = {
            {20, 2, 24, 9},
            {22, 8, 11, 14},
        };
        assert(quads == expected);
    }

    std::cout << to_string_u128(S(10000)) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

public class Euler557 {
    public static String solve() {
        int n = 10000;
        long total = 0;
        for (int a = 1; a <= n - 3; a++) {
            if (n - a < 3)
                continue;
            int bmax = Math.min(a - 1, (n - a - 1) / 2);
            long bn = (long) a * (n - a), bd = a + n;
            bmax = Math.min(bmax, (int) (bn / bd));
            if (bmax < 1)
                continue;
            long a2 = (long) a * a;
            for (int b = 1; b <= bmax; b++) {
                int rem = n - a - b;
                if (rem <= 2)
                    continue;
                long rhs = a2 * rem;
                long lc = (long) a * (a + b) + (long) n * b;
                long cmaxT = rhs / lc;
                long cmaxBc = (a2 - 1) / b;
                long cmaxSum = n - a - b - 1;
                long cmax = Math.min(cmaxT, Math.min(cmaxBc, cmaxSum));
                if (cmax < b)
                    continue;
                long base = (long) a * (a + b);
                for (long c = b; c <= cmax; c++) {
                    long bc = (long) b * c;
                    long denom = a2 - bc;
                    if (denom <= 0)
                        break;
                    long numer = base * (a + c);
                    if (numer % denom != 0)
                        continue;
                    long T = numer / denom;
                    if (T > n)
                        continue;
                    long d = T - a - b - c;
                    if (d <= 0)
                        continue;
                    total += T;
                }
            }
        }
        return String.valueOf(total);
    }

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