Problem 842: Irregular Star Polygons

Problem 842 source code

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <thread>
#include <unordered_map>
#include <vector>
#include <boost/multiprecision/cpp_dec_float.hpp>

using namespace std;

namespace {

constexpr int64_t MOD = 1'000'000'007LL;
using BigFloat = boost::multiprecision::cpp_dec_float_50;
const BigFloat SCALE_BF("1e18");

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

int64_t mod_pow(int64_t a, int64_t e) {
    int64_t r = 1 % MOD;
    int64_t x = (a % MOD + MOD) % MOD;
    while (e > 0) {
        if (e & 1) r = mod_mul(r, x);
        x = mod_mul(x, x);
        e >>= 1;
    }
    return r;
}

struct Precomp {
    vector<int64_t> fact;
    vector<int64_t> invfact;
    vector<int64_t> pow2;
    int64_t inv2;
};

Precomp precompute(int max_n) {
    Precomp p;
    p.fact.assign(max_n + 1, 1);
    p.invfact.assign(max_n + 1, 1);
    p.pow2.assign(max_n + 1, 1);
    for (int i = 1; i <= max_n; ++i) {
        p.fact[i] = mod_mul(p.fact[i - 1], i);
        p.pow2[i] = mod_mul(p.pow2[i - 1], 2);
    }
    p.invfact[max_n] = mod_pow(p.fact[max_n], MOD - 2);
    for (int i = max_n; i >= 1; --i) {
        p.invfact[i - 1] = mod_mul(p.invfact[i], i);
    }
    p.inv2 = mod_pow(2, MOD - 2);
    return p;
}

int64_t comb(int n, int k, const Precomp& p) {
    if (k < 0 || k > n) return 0;
    return mod_mul(p.fact[n], mod_mul(p.invfact[k], p.invfact[n - k]));
}

int64_t cycles_with_t_edges(int n, int t, const Precomp& p) {
    if (t == 0) return mod_mul(p.fact[n - 1], p.inv2);
    return mod_mul(p.fact[n - t - 1], p.pow2[t - 1]);
}

int64_t count_at_least_two(int n, int m, const Precomp& p) {
    int64_t total = cycles_with_t_edges(n, 0, p);
    int64_t none = 0;
    for (int t = 0; t <= m; ++t) {
        int64_t term = mod_mul(comb(m, t, p), cycles_with_t_edges(n, t, p));
        if (t & 1) {
            none = (none - term) % MOD;
        } else {
            none = (none + term) % MOD;
        }
    }
    if (none < 0) none += MOD;

    int64_t exactly_one = 0;
    for (int t = 0; t <= m - 1; ++t) {
        int64_t term = mod_mul(comb(m - 1, t, p), cycles_with_t_edges(n, t + 1, p));
        if (t & 1) {
            exactly_one = (exactly_one - term) % MOD;
        } else {
            exactly_one = (exactly_one + term) % MOD;
        }
    }
    if (exactly_one < 0) exactly_one += MOD;
    exactly_one = mod_mul(exactly_one, m % MOD);

    int64_t result = (total - none - exactly_one) % MOD;
    if (result < 0) result += MOD;
    return result;
}

struct PointKey {
    int64_t x;
    int64_t y;

    bool operator==(const PointKey& other) const noexcept {
        return x == other.x && y == other.y;
    }
};

struct PointKeyHash {
    size_t operator()(const PointKey& p) const noexcept {
        size_t h1 = std::hash<int64_t>{}(p.x);
        size_t h2 = std::hash<int64_t>{}(p.y);
        return h1 ^ (h2 + 0x9e3779b97f4a7c15ULL + (h1 << 6) + (h1 >> 2));
    }
};

int64_t choose4(int n) {
    if (n < 4) return 0;
    return static_cast<int64_t>(n) * (n - 1) * (n - 2) * (n - 3) / 24;
}

int64_t round_to_int64(const BigFloat& v) {
    BigFloat scaled = v * SCALE_BF;
    if (scaled >= 0) {
        return static_cast<int64_t>(scaled + BigFloat("0.5"));
    }
    return -static_cast<int64_t>(-scaled + BigFloat("0.5"));
}

vector<int64_t> intersection_multiplicities(int n) {
    using boost::multiprecision::acos;
    using boost::multiprecision::cos;
    using boost::multiprecision::sin;
    vector<pair<BigFloat, BigFloat>> pts(n);
    const BigFloat pi = acos(BigFloat(-1));
    for (int i = 0; i < n; ++i) {
        BigFloat ang = BigFloat(2) * pi * i / n;
        pts[i] = {cos(ang), sin(ang)};
    }

    unordered_map<PointKey, int, PointKeyHash> counts;
    counts.reserve(static_cast<size_t>(choose4(n) + 1));
    counts.max_load_factor(0.7f);

    for (int a = 0; a < n; ++a) {
        for (int b = a + 1; b < n; ++b) {
            for (int c = b + 1; c < n; ++c) {
                for (int d = c + 1; d < n; ++d) {
                    const auto& A = pts[a];
                    const auto& C = pts[c];
                    const auto& B = pts[b];
                    const auto& D = pts[d];

                    BigFloat x1 = A.first, y1 = A.second;
                    BigFloat x2 = C.first, y2 = C.second;
                    BigFloat x3 = B.first, y3 = B.second;
                    BigFloat x4 = D.first, y4 = D.second;

                    BigFloat den = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);
                    BigFloat det1 = x1 * y2 - y1 * x2;
                    BigFloat det2 = x3 * y4 - y3 * x4;
                    BigFloat px = (det1 * (x3 - x4) - (x1 - x2) * det2) / den;
                    BigFloat py = (det1 * (y3 - y4) - (y1 - y2) * det2) / den;

                    // Quantize high-precision coordinates for stable hashing.
                    PointKey key{
                        round_to_int64(px),
                        round_to_int64(py)
                    };
                    ++counts[key];
                }
            }
        }
    }

    vector<int64_t> Nm(n / 2 + 2, 0);
    for (const auto& kv : counts) {
        int q = kv.second;
        int64_t D = 1 + 8LL * q;
        int64_t s = static_cast<int64_t>(llroundl(sqrtl(static_cast<long double>(D))));
        while (s * s < D) ++s;
        while (s * s > D) --s;
        if (s * s != D) {
            cerr << "Validation failure: non-triangular count for n=" << n << '\n';
            continue;
        }
        int m = static_cast<int>((1 + s) / 2);
        if (m >= static_cast<int>(Nm.size())) Nm.resize(m + 1, 0);
        Nm[m] += 1;
    }
    return Nm;
}

int64_t compute_T(int n, const Precomp& p) {
    vector<int64_t> Nm = intersection_multiplicities(n);
    int64_t total = 0;
    for (int m = 2; m < static_cast<int>(Nm.size()); ++m) {
        if (Nm[m] == 0) continue;
        int64_t ways = count_at_least_two(n, m, p);
        total = (total + mod_mul(Nm[m] % MOD, ways)) % MOD;
    }
    return total;
}

}  // namespace

int main() {
    const int N_MIN = 3;
    const int N_MAX = 60;
    Precomp pre = precompute(N_MAX);

    vector<int64_t> Tvals(N_MAX + 1, 0);
    int thread_count = static_cast<int>(thread::hardware_concurrency());
    if (thread_count <= 0) thread_count = 1;
    vector<int64_t> partial(thread_count, 0);

    atomic<int> next_n(N_MIN);
    vector<thread> workers;
    workers.reserve(thread_count);
    for (int t = 0; t < thread_count; ++t) {
        workers.emplace_back([&, t]() {
            for (;;) {
                int n = next_n.fetch_add(1);
                if (n > N_MAX) break;
                int64_t Tn = compute_T(n, pre);
                Tvals[n] = Tn;
                partial[t] = (partial[t] + Tn) % MOD;
            }
        });
    }
    for (auto& th : workers) th.join();

    if (Tvals[5] != 20 || Tvals[8] != 14640) {
        cerr << "Validation failure: T(5) or T(8) mismatch\n";
        return 1;
    }

    int64_t answer = 0;
    for (int t = 0; t < thread_count; ++t) {
        answer = (answer + partial[t]) % MOD;
    }
    cout << answer << '\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

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler842 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler842.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler842_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler842 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler842 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler842 C++ bridge produced empty output.");
        }
        return parsed;
    }

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