Problem 729: Range of Periodic Sequence

Problem 729 source code

C++

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

namespace {

using u32 = std::uint32_t;
using i64 = std::int64_t;

struct Options {
    int p = 25;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

struct KahanSum {
    long double sum = 0.0L;
    long double c = 0.0L;

    void add(const long double x) {
        const long double y = x - c;
        const long double t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }
};

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

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

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

    value = static_cast<u32>(parsed);
    return true;
}

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

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

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u32 parsed_u32 = 0U;
        if (parse_u32_after_prefix(arg, "--p=", parsed_u32)) {
            if (parsed_u32 > static_cast<u32>(std::numeric_limits<int>::max())) {
                std::cerr << "--p is too large.\n";
                return false;
            }
            options.p = static_cast<int>(parsed_u32);
            continue;
        }

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

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

    if (options.p < 2) {
        std::cerr << "--p must be at least 2.\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload) {
    if (!allow_multithreading || workload < 2048ULL) {
        return 1U;
    }

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

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

long double apply_inverse_word(long double y, const u32 bits, const int m) {
    for (int i = 0; i < m; ++i) {
        const long double d = std::sqrt(y * y + 4.0L);
        if (((bits >> i) & 1U) != 0U) {
            y = 0.5L * (y + d);
        } else {
            y = 0.5L * (y - d);
        }
    }
    return y;
}

long double solve_fixed_point(const u32 bits, const int m) {
    long double y = 0.0L;

    for (int it = 0; it < 80; ++it) {
        const long double ny = apply_inverse_word(y, bits, m);
        if (std::fabsl(ny - y) < 1e-28L) {
            return ny;
        }
        y = ny;
    }

    for (int it = 0; it < 16; ++it) {
        long double v = y;
        long double deriv = 1.0L;
        for (int i = 0; i < m; ++i) {
            const long double d = std::sqrt(v * v + 4.0L);
            const int s = (((bits >> i) & 1U) != 0U) ? 1 : -1;
            deriv *= 0.5L * (1.0L + static_cast<long double>(s) * (v / d));
            v = 0.5L * (v + static_cast<long double>(s) * d);
        }

        const long double phi = v - y;
        const long double dphi = deriv - 1.0L;
        if (std::fabsl(dphi) < 1e-30L) {
            break;
        }
        const long double step = phi / dphi;
        y -= step;
        if (std::fabsl(step) < 1e-28L) {
            break;
        }
    }

    return y;
}

long double orbit_range_from_inverse_word(const long double x0, const u32 bits, const int m) {
    long double mn = x0;
    long double mx = x0;
    long double y = x0;

    // The fixed point x0 of g_{s_{m-1}} ... g_{s_0} generates the same orbit
    // by repeatedly applying g_{s_0}, g_{s_1}, ..., g_{s_{m-2}}.
    for (int i = 0; i < m - 1; ++i) {
        const long double d = std::sqrt(y * y + 4.0L);
        if (((bits >> i) & 1U) != 0U) {
            y = 0.5L * (y + d);
        } else {
            y = 0.5L * (y - d);
        }
        mn = std::min(mn, y);
        mx = std::max(mx, y);
    }

    return mx - mn;
}

int mobius(const int n) {
    int x = n;
    int prime_count = 0;
    for (int p = 2; p * p <= x; ++p) {
        if (x % p != 0) {
            continue;
        }
        int exp = 0;
        while (x % p == 0) {
            x /= p;
            ++exp;
        }
        if (exp > 1) {
            return 0;
        }
        ++prime_count;
    }
    if (x > 1) {
        ++prime_count;
    }
    return (prime_count % 2 == 0) ? 1 : -1;
}

i64 expected_primitive_necklaces(const int n) {
    i64 sum = 0;
    for (int d = 1; d <= n; ++d) {
        if (n % d != 0) {
            continue;
        }
        const int mu = mobius(d);
        if (mu == 0) {
            continue;
        }
        sum += static_cast<i64>(mu) * (1LL << (n / d));
    }
    return sum / n;
}

void generate_primitive_necklaces(const int n, std::vector<u32>& out) {
    std::vector<int> a(static_cast<std::size_t>(n) + 1ULL, 0);
    out.clear();
    out.reserve(static_cast<std::size_t>(expected_primitive_necklaces(n)));

    std::function<void(int, int)> gen = [&](const int t, const int p) {
        if (t > n) {
            if (n % p == 0 && p == n) {
                u32 bits = 0U;
                for (int i = 1; i <= n; ++i) {
                    bits |= static_cast<u32>(a[i]) << static_cast<u32>(i - 1);
                }
                out.push_back(bits);
            }
            return;
        }

        a[t] = a[t - p];
        gen(t + 1, p);

        for (int j = a[t - p] + 1; j < 2; ++j) {
            a[t] = j;
            gen(t + 1, t);
        }
    };

    gen(1, 1);
}

long double solve_for_p(const int p,
                        const bool allow_multithreading,
                        const unsigned requested_threads,
                        const bool validate_counts) {
    KahanSum total;

    for (int m = 2; m <= p; ++m) {
        std::vector<u32> necklaces;
        generate_primitive_necklaces(m, necklaces);

        if (validate_counts) {
            const i64 expected = expected_primitive_necklaces(m);
            if (static_cast<i64>(necklaces.size()) != expected) {
                std::cerr << "Necklace count mismatch at m=" << m << ": expected " << expected
                          << ", got " << necklaces.size() << "\n";
                std::exit(1);
            }
        }

        const unsigned threads =
            choose_thread_count(allow_multithreading, requested_threads, necklaces.size());

        if (threads == 1U) {
            for (const u32 bits : necklaces) {
                const long double x0 = solve_fixed_point(bits, m);
                const long double r = orbit_range_from_inverse_word(x0, bits, m);
                total.add(static_cast<long double>(m) * r);
            }
            continue;
        }

        std::vector<KahanSum> partials(threads);
        auto worker = [&](const unsigned tid) {
            const std::size_t begin = (necklaces.size() * tid) / threads;
            const std::size_t end = (necklaces.size() * (tid + 1U)) / threads;
            for (std::size_t idx = begin; idx < end; ++idx) {
                const u32 bits = necklaces[idx];
                const long double x0 = solve_fixed_point(bits, m);
                const long double r = orbit_range_from_inverse_word(x0, bits, m);
                partials[tid].add(static_cast<long double>(m) * r);
            }
        };

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

        for (const auto& ps : partials) {
            total.add(ps.sum);
            total.add(ps.c);
        }
    }

    return total.sum;
}

bool run_checkpoints(const bool allow_multithreading, const unsigned requested_threads) {
    struct Checkpoint {
        int p;
        long double expected;
    };

    const std::vector<Checkpoint> checks = {
        {2, 2.8284271247461900976L},
        {3, 14.646120160892684L},
        {5, 124.10555788745795L},
    };

    for (const auto& cp : checks) {
        const long double got = solve_for_p(cp.p, allow_multithreading, requested_threads, true);
        const long double err = std::fabsl(got - cp.expected);
        if (err > 5e-10L) {
            std::cerr << std::setprecision(18)
                      << "Checkpoint failed for P=" << cp.p << ": expected " << cp.expected
                      << ", got " << got << ", |err|=" << err << "\n";
            return false;
        }
    }

    std::cerr << "Validation checkpoints passed.\n";
    return true;
}

}  // 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_checkpoints &&
        !run_checkpoints(options.allow_multithreading, options.requested_threads)) {
        return 1;
    }

    const long double ans =
        solve_for_p(options.p, options.allow_multithreading, options.requested_threads, false);

    std::cout << std::fixed << std::setprecision(10) << ans << '\n';
    return 0;
}

Python

import math
import multiprocessing

class KahanSum:
    def __init__(self):
        self.sum = 0.0
        self.c = 0.0

    def add(self, x):
        y = x - self.c
        t = self.sum + y
        self.c = (t - self.sum) - y
        self.sum = t

def apply_inverse_word(y, bits, m):
    for i in range(m):
        d = math.sqrt(y * y + 4.0)
        if (bits >> i) & 1:
            y = 0.5 * (y + d)
        else:
            y = 0.5 * (y - d)
    return y

def solve_fixed_point(bits, m):
    y = 0.0
    for it in range(80):
        ny = apply_inverse_word(y, bits, m)
        if abs(ny - y) < 1e-15:
            return ny
        y = ny
        
    for it in range(16):
        v = y
        deriv = 1.0
        for i in range(m):
            d = math.sqrt(v * v + 4.0)
            s = 1 if (bits >> i) & 1 else -1
            deriv *= 0.5 * (1.0 + s * (v / d))
            v = 0.5 * (v + s * d)
            
        phi = v - y
        dphi = deriv - 1.0
        if abs(dphi) < 1e-15:
            break
            
        step = phi / dphi
        y -= step
        if abs(step) < 1e-15:
            break
            
    return y

def orbit_range_from_inverse_word(x0, bits, m):
    mn = x0
    mx = x0
    y = x0
    
    for i in range(m - 1):
        d = math.sqrt(y * y + 4.0)
        if (bits >> i) & 1:
            y = 0.5 * (y + d)
        else:
            y = 0.5 * (y - d)
        if y < mn: mn = y
        if y > mx: mx = y
        
    return mx - mn

def mobius(n):
    x = n
    prime_count = 0
    p = 2
    while p * p <= x:
        if x % p == 0:
            exp = 0
            while x % p == 0:
                x //= p
                exp += 1
            if exp > 1: return 0
            prime_count += 1
        p += 1
    if x > 1:
        prime_count += 1
    return 1 if prime_count % 2 == 0 else -1

def expected_primitive_necklaces(n):
    total = 0
    for d in range(1, n + 1):
        if n % d == 0:
            mu = mobius(d)
            if mu == 0: continue
            total += mu * (1 << (n // d))
    return total // n

def generate_primitive_necklaces(n):
    out = []
    a = [0] * (n + 1)
    
    def gen(t, p):
        if t > n:
            if n % p == 0 and p == n:
                bits = 0
                for i in range(1, n + 1):
                    bits |= (a[i] << (i - 1))
                out.append(bits)
            return
            
        a[t] = a[t - p]
        gen(t + 1, p)
        
        for j in range(a[t - p] + 1, 2):
            a[t] = j
            gen(t + 1, t)
            
    gen(1, 1)
    return out

def worker(m, necklaces):
    total = KahanSum()
    for bits in necklaces:
        x0 = solve_fixed_point(bits, m)
        r = orbit_range_from_inverse_word(x0, bits, m)
        total.add(m * r)
    return total.sum

def solve():
    p = 25
    total = KahanSum()
    
    threads = max(1, multiprocessing.cpu_count())
    pool = multiprocessing.Pool(threads)
    
    for m in range(2, p + 1):
        necklaces = generate_primitive_necklaces(m)
        
        chunk_size = max(1, len(necklaces) // threads)
        chunks = [necklaces[i:i + chunk_size] for i in range(0, len(necklaces), chunk_size)]
        
        results = pool.starmap(worker, [(m, chunk) for chunk in chunks])
        for val in results:
            total.add(val)
            
    pool.close()
    
    return f"{total.sum:.10f}"

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

Java

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

public class Euler729 {
    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("Euler729.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(".euler729_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 Euler729 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("Euler729 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("Euler729 C++ bridge produced empty output.");
        }
        return parsed;
    }

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