euler/967.py

#!/usr/bin/env pypy3

# Generated by GPT-5 on 2025-11-05. Tested on Ubuntu 25.04, pypy3 3.11.11, runtime: 0m1.695s

# Python implementation of the meet-in-the-middle algorithm for Project Euler 967 (B-trivisible)
# This will print the example values and the final answer F(10**18, 120).
# It is written to be reasonably efficient in CPython.
import bisect, sys
sys.setrecursionlimit(10000)

def primes_upto(B):
    sieve = bytearray(b'\x01')*(B+1)
    sieve[0:2] = b'\x00\x00'
    for i in range(2, int(B**0.5)+1):
        if sieve[i]:
            step = i
            start = i*i
            sieve[start:B+1:step] = b'\x00'*(((B - start)//step) + 1)
    return [p for p in range(2, B+1) if sieve[p]]

def extend_vec_mod3(cls, s0, s1, s2):
    # cls is prime % 3: 1 or 2
    a,b,c = s0, s1, s2
    if cls == 1:
        return a - c, b - a, c - b
    else:
        return a - b, b - c, c - a

def gen_half(primes, N):
    # Enumerate squarefree products <= N from this list of primes.
    # Return list of tuples: (prod, t, s0, s1, s2)
    out = []
    # iterative DFS stack: index, prod, t, s0,s1,s2
    stack = [(0, 1, 1, 1, 0, 0)]
    while stack:
        i, prod, t, s0, s1, s2 = stack.pop()
        out.append((prod, t, s0, s1, s2))
        for j in range(i, len(primes)):
            p = primes[j]
            np = prod * p
            if np > N:
                # since primes is ascending, further j will only be larger
                break
            ns0, ns1, ns2 = extend_vec_mod3(p % 3, s0, s1, s2)
            stack.append((j+1, np, -t, ns0, ns1, ns2))
    return out

def solve(N, B):
    allp = primes_upto(B)
    # remove 3 (contributes 0 mod 3, and can be ignored per analysis)
    allp = [p for p in allp if p != 3]
    mid = len(allp) // 2
    L = allp[:mid]
    R = allp[mid:]
    A = gen_half(L, N)
    Bv = gen_half(R, N)
    # sort Bv by prod
    Bv.sort(key=lambda x: x[0])
    bp = [x[0] for x in Bv]
    m = len(Bv)
    pref0 = [0] * (m+1)
    pref1 = [0] * (m+1)
    pref2 = [0] * (m+1)
    for i,(prod,t,s0,s1,s2) in enumerate(Bv, start=1):
        pref0[i] = pref0[i-1] + t * s0
        pref1[i] = pref1[i-1] + t * s1
        pref2[i] = pref2[i-1] + t * s2

    def sum_range(l, r):
        return (pref0[r] - pref0[l], pref1[r] - pref1[l], pref2[r] - pref2[l])

    ans = 0
    for a_prod, a_t, a_s0, a_s1, a_s2 in A:
        if a_prod > N: 
            continue
        limit = N // a_prod
        # find number of right items with b <= limit
        r = bisect.bisect_right(bp, limit)
        i = 0
        # iterate in blocks where floor(N/(a_prod*b)) is constant
        while i < r:
            v = limit // bp[i]
            if v == 0:
                break
            maxb = limit // v
            j = bisect.bisect_right(bp, maxb, i, r)
            s0, s1, s2 = sum_range(i, j)
            comb = a_s0 * s0 + a_s1 * s2 + a_s2 * s1
            ans += a_t * v * comb
            i = j
    return ans

if __name__ == "__main__":
    # examples
    print(solve(10, 4))    # expect 5
    print(solve(10, 10))   # expect 3
    print(solve(100, 10))  # expect 41

    # final requested value
    N = 10**18
    B = 120
    print(solve(N, B))