Problem 471: Triangle Inscribed in Ellipse

Problem 471 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <sstream>
#include <string>
#include <thread>
#include <vector>

namespace {

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

constexpr u64 kDefaultN = 100'000'000'000ULL;
constexpr u64 kDefaultValidationN = 2'000ULL;
constexpr u64 kDirectHarmonicLimit = 5'000'000ULL;

// From the geometry of the problem (t = a / b), one gets:
// r(a,b) = b * (t - 2) / (t - 1) = b * (a - 2b) / (a - b).
long double r_formula(u64 a, u64 b) {
    return static_cast<long double>(b) * static_cast<long double>(a - 2 * b) /
           static_cast<long double>(a - b);
}

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

unsigned choose_thread_count(bool allow_multithreading,
                             unsigned requested_threads,
                             std::size_t workload) {
    if (!allow_multithreading || workload < 2'000'000ULL) {
        return 1;
    }

    unsigned threads = requested_threads;
    if (threads == 0) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0) {
            threads = 1;
        }
    }

    return std::max(1u, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}

long double harmonic_asymptotic(u64 n) {
    if (n == 0) {
        return 0.0L;
    }

    // H_n = log(n) + gamma + 1/(2n) - 1/(12n^2) + 1/(120n^4) - ...
    constexpr long double kEulerGamma = 0.577215664901532860606512090082402431L;

    const long double x = static_cast<long double>(n);
    const long double inv = 1.0L / x;
    const long double inv2 = inv * inv;
    const long double inv4 = inv2 * inv2;
    const long double inv6 = inv4 * inv2;
    const long double inv8 = inv4 * inv4;
    const long double inv10 = inv8 * inv2;

    return std::log(x) + kEulerGamma + inv / 2.0L - inv2 / 12.0L + inv4 / 120.0L -
           inv6 / 252.0L + inv8 / 240.0L - 5.0L * inv10 / 660.0L;
}

long double harmonic_range(u64 left, u64 right) {
    if (left > right) {
        return 0.0L;
    }

    const u64 width = right - left + 1ULL;
    if (right <= kDirectHarmonicLimit || width <= kDirectHarmonicLimit) {
        long double sum = 0.0L;
        for (u64 k = left; k <= right; ++k) {
            sum += 1.0L / static_cast<long double>(k);
        }
        return sum;
    }

    return harmonic_asymptotic(right) - harmonic_asymptotic(left - 1ULL);
}

long double compute_g_fast(u64 n) {
    if (n < 3) {
        return 0.0L;
    }

    const u64 q = n / 2ULL;
    const u64 m = (n - 1ULL) / 2ULL;
    const u64 left = n - m;
    const u64 right = n - 1ULL;

    const long double qd = static_cast<long double>(q);
    const long double md = static_cast<long double>(m);
    const long double nd = static_cast<long double>(n);
    const long double ld = static_cast<long double>(left);
    const long double rd = static_cast<long double>(right);

    // Split by c = a - b:
    // c <= floor(n/2): closed polynomial contribution.
    const long double part_low =
        qd * (2.0L * qd * qd + 3.0L * qd - 5.0L) / 36.0L;

    // c > floor(n/2), with x = n - c in [1..m].
    const long double part_high_poly = md * (md + 1.0L) * (md + 2.0L) / 6.0L;

    const long double harmonic = harmonic_range(left, right);

    const long double acoef = nd * (2.0L * nd + 1.0L) * (nd + 1.0L);
    const long double bcoef = -(6.0L * nd * nd + 6.0L * nd + 1.0L);
    const long double ccoef = 6.0L * nd + 3.0L;

    const long double sum_k = (ld + rd) * md / 2.0L;
    const long double sum_k2 =
        (rd * (rd + 1.0L) * (2.0L * rd + 1.0L) -
         (ld - 1.0L) * ld * (2.0L * ld - 1.0L)) /
        6.0L;

    const long double rational_tail =
        acoef * harmonic + bcoef * md + ccoef * sum_k - 2.0L * sum_k2;

    return part_low + part_high_poly - rational_tail / 6.0L;
}

long double compute_g_direct_parallel(u64 n,
                                      bool allow_multithreading,
                                      unsigned requested_threads) {
    if (n < 3) {
        return 0.0L;
    }

    const long double pair_estimate =
        static_cast<long double>(n) * static_cast<long double>(n) / 4.0L;
    const std::size_t workload =
        (pair_estimate >= static_cast<long double>(std::numeric_limits<std::size_t>::max()))
            ? std::numeric_limits<std::size_t>::max()
            : static_cast<std::size_t>(pair_estimate);

    const unsigned threads = choose_thread_count(allow_multithreading, requested_threads, workload);
    if (threads == 1) {
        long double total = 0.0L;
        for (u64 a = 3; a <= n; ++a) {
            const u64 bmax = (a - 1ULL) / 2ULL;
            for (u64 b = 1; b <= bmax; ++b) {
                total += r_formula(a, b);
            }
        }
        return total;
    }

    const u64 a_count = n - 2ULL;
    std::vector<std::thread> workers;
    std::vector<long double> partial(threads, 0.0L);
    workers.reserve(threads);

    for (unsigned t = 0; t < threads; ++t) {
        const u64 a_begin = 3ULL + static_cast<u64>((static_cast<u128>(a_count) * t) / threads);
        const u64 a_end =
            3ULL + static_cast<u64>((static_cast<u128>(a_count) * (t + 1ULL)) / threads);

        workers.emplace_back([&, t, a_begin, a_end]() {
            long double local = 0.0L;
            for (u64 a = a_begin; a < a_end; ++a) {
                const u64 bmax = (a - 1ULL) / 2ULL;
                for (u64 b = 1; b <= bmax; ++b) {
                    local += r_formula(a, b);
                }
            }
            partial[t] = local;
        });
    }

    for (auto& worker : workers) {
        worker.join();
    }

    long double total = 0.0L;
    for (const long double v : partial) {
        total += v;
    }
    return total;
}

bool nearly_equal(long double lhs,
                  long double rhs,
                  long double rel_tol = 1e-12L,
                  long double abs_tol = 1e-12L) {
    const long double diff = std::fabsl(lhs - rhs);
    if (diff <= abs_tol) {
        return true;
    }
    const long double scale = std::max(std::fabsl(lhs), std::fabsl(rhs));
    return diff <= rel_tol * scale;
}

std::string to_scientific_10_significant(long double value) {
    if (value == 0.0L) {
        return "0.000000000e0";
    }

    const bool negative = (value < 0.0L);
    long double x = negative ? -value : value;

    int exponent = static_cast<int>(std::floor(std::log10(x)));
    long double mantissa = x / std::powl(10.0L, static_cast<long double>(exponent));

    long double scaled = std::floor(mantissa * 1'000'000'000.0L + 0.5L);
    if (scaled >= 10'000'000'000.0L) {
        scaled /= 10.0L;
        ++exponent;
    }

    const u64 digits = static_cast<u64>(scaled);
    const u64 lead = digits / 1'000'000'000ULL;
    const u64 frac = digits % 1'000'000'000ULL;

    std::ostringstream oss;
    if (negative) {
        oss << '-';
    }
    oss << lead << '.' << std::setw(9) << std::setfill('0') << frac << 'e' << exponent;
    return oss.str();
}

struct Options {
    u64 n = kDefaultN;
    bool run_validation = true;
    u64 validation_n = kDefaultValidationN;
    bool allow_multithreading = true;
    unsigned requested_threads = 0;
    bool answer_only = false;
};

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-validation") {
            options.run_validation = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }
        if (arg == "--answer-only") {
            options.answer_only = true;
            continue;
        }

        u64 parsed_u64 = 0;
        if (parse_u64_after_prefix(arg, "--n=", parsed_u64)) {
            options.n = parsed_u64;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--validation-n=", parsed_u64)) {
            options.validation_n = parsed_u64;
            continue;
        }

        unsigned parsed_threads = 0;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_threads)) {
            options.requested_threads = parsed_threads;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return true;
}

bool run_checkpoints(const Options& options) {
    bool ok = true;

    auto check = [&](const bool condition, const std::string& label) {
        std::cout << "[checkpoint] " << label << ": " << (condition ? "PASS" : "FAIL") << '\n';
        if (!condition) {
            ok = false;
        }
    };

    check(nearly_equal(r_formula(3, 1), 0.5L, 1e-15L, 1e-15L), "r(3,1) = 1/2");
    check(nearly_equal(r_formula(6, 2), 1.0L, 1e-15L, 1e-15L), "r(6,2) = 1");
    check(nearly_equal(r_formula(12, 3), 2.0L, 1e-15L, 1e-15L), "r(12,3) = 2");

    const std::string g10 = to_scientific_10_significant(compute_g_fast(10));
    const std::string g100 = to_scientific_10_significant(compute_g_fast(100));

    check(g10 == "2.059722222e1", "G(10) scientific form");
    check(g100 == "1.922360980e4", "G(100) scientific form");

    if (options.validation_n >= 3) {
        const long double fast = compute_g_fast(options.validation_n);
        const long double direct = compute_g_direct_parallel(
            options.validation_n, options.allow_multithreading, options.requested_threads);
        const bool consistent = nearly_equal(fast, direct, 1e-12L, 1e-9L);

        std::ostringstream label;
        label << "fast vs direct at n=" << options.validation_n;
        check(consistent, label.str());

        if (!consistent) {
            std::cout << std::setprecision(20);
            std::cout << "  fast   = " << fast << '\n';
            std::cout << "  direct = " << direct << '\n';
        }
    }

    return ok;
}

} // namespace

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

    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_validation && !options.answer_only) {
        if (!run_checkpoints(options)) {
            std::cerr << "Validation failed.\n";
            return 1;
        }
    }

    const long double answer = compute_g_fast(options.n);
    const std::string answer_scientific = to_scientific_10_significant(answer);

    if (!options.answer_only) {
        std::cout << "[result] G(" << options.n << ") = " << answer_scientific << '\n';
    } else {
        std::cout << answer_scientific << '\n';
    }

    return 0;
}

Python

import sys
import math
from decimal import Decimal, getcontext

# Set sufficient precision for asymptotic expansion
getcontext().prec = 60

def harmonic_asymptotic(n_val):
    if n_val == 0:
        return Decimal(0)
    
    # H_n = log(n) + gamma + 1/(2n) - 1/(12n^2) + 1/(120n^4) - 1/(252n^6) ...
    euler_gamma = Decimal('0.57721566490153286060651209008240243104215933593992')
    
    x = Decimal(n_val)
    inv = Decimal(1) / x
    inv2 = inv * inv
    inv4 = inv2 * inv2
    inv6 = inv4 * inv2
    inv8 = inv4 * inv4
    inv10 = inv8 * inv2
    
    log_x = Decimal(math.log(n_val))
    # Python decimal ln is .ln()
    log_x = x.ln()

    res = log_x + euler_gamma + inv / Decimal(2) \
          - inv2 / Decimal(12) + inv4 / Decimal(120) \
          - inv6 / Decimal(252) + inv8 / Decimal(240) \
          - Decimal(5) * inv10 / Decimal(660)
          
    return res

def harmonic_range(left, right):
    if left > right:
        return Decimal(0)
    
    kDirectHarmonicLimit = 5_000_000
    width = right - left + 1
    
    if right <= kDirectHarmonicLimit or width <= kDirectHarmonicLimit:
        total = Decimal(0)
        for k in range(left, right + 1):
            total += Decimal(1) / Decimal(k)
        return total
    
    return harmonic_asymptotic(right) - harmonic_asymptotic(left - 1)

def compute_g_fast(n):
    if n < 3:
        return Decimal(0)
        
    q = n // 2
    m = (n - 1) // 2
    left = n - m
    right = n - 1
    
    qd = Decimal(q)
    md = Decimal(m)
    nd = Decimal(n)
    ld = Decimal(left)
    rd = Decimal(right)
    
    part_low = qd * (Decimal(2) * qd * qd + Decimal(3) * qd - Decimal(5)) / Decimal(36)
    part_high_poly = md * (md + Decimal(1)) * (md + Decimal(2)) / Decimal(6)
    
    harmonic = harmonic_range(left, right)
    
    acoef = nd * (Decimal(2) * nd + Decimal(1)) * (nd + Decimal(1))
    bcoef = -(Decimal(6) * nd * nd + Decimal(6) * nd + Decimal(1))
    ccoef = Decimal(6) * nd + Decimal(3)
    
    sum_k = (ld + rd) * md / Decimal(2)
    sum_k2 = (rd * (rd + Decimal(1)) * (Decimal(2) * rd + Decimal(1)) -
             (ld - Decimal(1)) * ld * (Decimal(2) * ld - Decimal(1))) / Decimal(6)
             
    rational_tail = acoef * harmonic + bcoef * md + ccoef * sum_k - Decimal(2) * sum_k2
    
    return part_low + part_high_poly - rational_tail / Decimal(6)

def to_scientific_10_significant(value):
    if value == Decimal(0):
        return "0.000000000e0"
        
    negative = value < 0
    x = -value if negative else value
    
    exponent = int(math.floor(math.log10(float(x))))
    mantissa = x / (Decimal(10) ** exponent)
    
    # Scale exactly like C++: floor(mantissa * 10^9 + 0.5)
    scaled = math.floor(float(mantissa) * 1_000_000_000.0 + 0.5)
    if scaled >= 10_000_000_000:
        scaled //= 10
        exponent += 1
        
    lead = scaled // 1_000_000_000
    frac = scaled % 1_000_000_000
    
    sign_str = "-" if negative else ""
    return f"{sign_str}{lead}.{frac:09d}e{exponent}"

def solve():
    n = 100_000_000_000
    ans = compute_g_fast(n)
    return to_scientific_10_significant(ans)

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

Java

import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;

public class Euler471 {

    private static final MathContext MC = new MathContext(60, RoundingMode.HALF_UP);
    private static final BigDecimal EULER_GAMMA = new BigDecimal(
            "0.57721566490153286060651209008240243104215933593992");

    private static double log(BigDecimal val) {
        return Math.log(val.doubleValue());
    }

    private static BigDecimal harmonicAsymptotic(long nVal) {
        if (nVal == 0)
            return BigDecimal.ZERO;

        BigDecimal x = new BigDecimal(nVal);
        BigDecimal inv = BigDecimal.ONE.divide(x, MC);
        BigDecimal inv2 = inv.multiply(inv, MC);
        BigDecimal inv4 = inv2.multiply(inv2, MC);
        BigDecimal inv6 = inv4.multiply(inv2, MC);
        BigDecimal inv8 = inv4.multiply(inv4, MC);
        BigDecimal inv10 = inv8.multiply(inv2, MC);

        BigDecimal logX = new BigDecimal(log(x));

        BigDecimal res = logX.add(EULER_GAMMA, MC)
                .add(inv.divide(new BigDecimal(2), MC), MC)
                .subtract(inv2.divide(new BigDecimal(12), MC), MC)
                .add(inv4.divide(new BigDecimal(120), MC), MC)
                .subtract(inv6.divide(new BigDecimal(252), MC), MC)
                .add(inv8.divide(new BigDecimal(240), MC), MC)
                .subtract(new BigDecimal(5).multiply(inv10, MC).divide(new BigDecimal(660), MC), MC);

        return res;
    }

    private static BigDecimal harmonicRange(long left, long right) {
        if (left > right)
            return BigDecimal.ZERO;

        long kDirectHarmonicLimit = 5_000_000;
        long width = right - left + 1;

        if (right <= kDirectHarmonicLimit || width <= kDirectHarmonicLimit) {
            BigDecimal total = BigDecimal.ZERO;
            for (long k = left; k <= right; ++k) {
                total = total.add(BigDecimal.ONE.divide(new BigDecimal(k), MC), MC);
            }
            return total;
        }

        return harmonicAsymptotic(right).subtract(harmonicAsymptotic(left - 1), MC);
    }

    private static BigDecimal computeGFast(long n) {
        if (n < 3)
            return BigDecimal.ZERO;

        long q = n / 2;
        long m = (n - 1) / 2;
        long left = n - m;
        long right = n - 1;

        BigDecimal qd = new BigDecimal(q);
        BigDecimal md = new BigDecimal(m);
        BigDecimal nd = new BigDecimal(n);
        BigDecimal ld = new BigDecimal(left);
        BigDecimal rd = new BigDecimal(right);

        BigDecimal partLow = qd.multiply(
                new BigDecimal(2).multiply(qd).multiply(qd)
                        .add(new BigDecimal(3).multiply(qd))
                        .subtract(new BigDecimal(5)))
                .divide(new BigDecimal(36), MC);

        BigDecimal partHighPoly = md.multiply(md.add(BigDecimal.ONE))
                .multiply(md.add(new BigDecimal(2)))
                .divide(new BigDecimal(6), MC);

        BigDecimal harmonic = harmonicRange(left, right);

        BigDecimal acoef = nd.multiply(new BigDecimal(2).multiply(nd).add(BigDecimal.ONE))
                .multiply(nd.add(BigDecimal.ONE));
        BigDecimal bcoef = new BigDecimal(6).multiply(nd).multiply(nd)
                .add(new BigDecimal(6).multiply(nd))
                .add(BigDecimal.ONE).negate();
        BigDecimal ccoef = new BigDecimal(6).multiply(nd).add(new BigDecimal(3));

        BigDecimal sumK = ld.add(rd).multiply(md).divide(new BigDecimal(2), MC);

        BigDecimal v1 = rd.multiply(rd.add(BigDecimal.ONE))
                .multiply(new BigDecimal(2).multiply(rd).add(BigDecimal.ONE));
        BigDecimal v2 = ld.subtract(BigDecimal.ONE).multiply(ld)
                .multiply(new BigDecimal(2).multiply(ld).subtract(BigDecimal.ONE));
        BigDecimal sumK2 = v1.subtract(v2).divide(new BigDecimal(6), MC);

        BigDecimal rationalTail = acoef.multiply(harmonic, MC)
                .add(bcoef.multiply(md, MC), MC)
                .add(ccoef.multiply(sumK, MC), MC)
                .subtract(new BigDecimal(2).multiply(sumK2, MC), MC);

        return partLow.add(partHighPoly, MC)
                .subtract(rationalTail.divide(new BigDecimal(6), MC), MC);
    }

    private static String toScientific10Significant(BigDecimal value) {
        if (value.compareTo(BigDecimal.ZERO) == 0)
            return "0.000000000e0";

        boolean negative = value.compareTo(BigDecimal.ZERO) < 0;
        double x = negative ? -value.doubleValue() : value.doubleValue();

        int exponent = (int) Math.floor(Math.log10(x));
        double mantissa = x / Math.pow(10.0, exponent);

        long scaled = (long) Math.floor(mantissa * 1_000_000_000.0 + 0.5);
        if (scaled >= 10_000_000_000L) {
            scaled /= 10L;
            exponent += 1;
        }

        long lead = scaled / 1_000_000_000L;
        long frac = scaled % 1_000_000_000L;

        String signStr = negative ? "-" : "";
        return String.format("%s%d.%09de%d", signStr, lead, frac, exponent);
    }

    public static void main(String[] args) {
        long n = 100_000_000_000L;
        BigDecimal ans = computeGFast(n);
        System.out.println(toScientific10Significant(ans));
    }
}