// @license Factorial Solver 8/20/2026
// https://raw.org/research/extending-the-verified-range-for-the-factorial-diophantine-equation/
//
// Copyright (c) 2026, Robert Eisele (https://raw.org/)
// Licensed under the MIT license.

//go:build ignore

package main

// Exact hybrid finite-field solver for the factorial Diophantine equation
//
//	A! B! = C!
//
// in the range 10^minDigits <= B <= 10^digits, assuming Habsieger's explicit
// bounds (Fibonacci Quarterly 57(1), 2019) for every hypothetical nontrivial
// solution other than (6,7,10).
//
// Revision 2 keeps every mathematical predicate of revision 1 unchanged and
// re-engineers only how the predicates are evaluated:
//
//   - dense finite-field images are shared between all workers and filled
//     with atomic bit-set operations, instead of one private q-bit image per
//     worker (which used to serialize the dense stage under a memory budget);
//   - the survivor test parallelizes over 64-aligned A segments using a
//     segmented prefix product for A! mod q, and scans only surviving bits;
//   - the exact decimal magnitude filter fills each per-k candidate interval
//     with whole-word operations (the candidate set of every k is a single
//     contiguous A interval, because L(A) and U(A) are nondecreasing);
//   - candidate bitsets are allocated only up to the largest possible A of
//     each k, and survivor counts are maintained incrementally;
//   - when memory is short, active k values are processed in groups sized by
//     the -filter-memory-mb budget instead of throttling workers.
//
// None of these changes can alter the survivor set: dense and sparse mode
// still decide the identical necessary condition A! mod q in image(f_k).

import (
	"flag"
	"fmt"
	"math"
	"math/big"
	"math/bits"
	"runtime"
	"sort"
	"sync"
	"sync/atomic"
	"time"
)

var one = big.NewInt(1)

type Solution struct {
	A int
	B *big.Int
	K int // C-B
}

type Stats struct {
	pairs        atomic.Uint64
	reducedRoots atomic.Uint64
	fallbacks    atomic.Uint64
	comparisons  atomic.Uint64
}

// theoremBounds returns padded integer bounds obtained from Habsieger's tighter
// inequalities for every hypothetical nontrivial solution other than (6,7,10)
// in the range B >= 10^3000:
//
//	A <= log(B+1)/log(2) + 2 log log(B+1)/log(2) - 1.3479
//	k <= log log(B+1)/log(2) - 0.8803,  k=C-B.
//
// For B <= 10^digits we safely use
//
//	log(B+1) < digits*log(10) + log(2).
func theoremBounds(digits int) (amax, kmax int) {
	// We round UP and add two spare integer units.  This intentionally makes
	// the computational search slightly larger; floating point is therefore
	// never used to exclude a boundary candidate.
	L := float64(digits)*math.Ln10 + math.Ln2
	aExpr := L/math.Ln2 + 2*math.Log(L)/math.Ln2 - 1.3479
	kExpr := math.Log(L)/math.Ln2 - 0.8803
	amax = int(math.Ceil(aExpr)) + 2
	kmax = int(math.Ceil(kExpr)) + 2
	if kmax < 2 {
		kmax = 2
	}
	return
}

// approxRoot gives a ~53-bit relative approximation to n^(1/k), solely as a
// Newton seed. Correctness never depends on floating point.
func approxRoot(n *big.Int, k int) *big.Int {
	L := n.BitLen()
	take := 53
	if L < take {
		take = L
	}

	shift := L - take
	topInt := new(big.Int).Rsh(new(big.Int).Set(n), uint(shift))
	top := float64(topInt.Uint64())

	log2n := math.Log2(top) + float64(shift)
	rlog := log2n / float64(k)
	e := int(math.Floor(rlog))
	frac := rlog - float64(e)

	const precision = 52
	mant := uint64(math.Exp2(frac) * float64(uint64(1)<<precision))
	if mant == 0 {
		mant = 1
	}

	x := new(big.Int).SetUint64(mant)
	if e >= precision {
		x.Lsh(x, uint(e-precision))
	} else {
		x.Rsh(x, uint(precision-e))
	}
	if x.Sign() == 0 {
		x.SetInt64(1)
	}
	return x
}

// nthRootFloor returns floor(n^(1/k)) exactly.
func nthRootFloor(n *big.Int, k int) *big.Int {
	if n.Sign() == 0 {
		return new(big.Int)
	}
	if k == 1 {
		return new(big.Int).Set(n)
	}
	if k == 2 {
		return new(big.Int).Sqrt(n)
	}

	x := approxRoot(n, k)
	y := new(big.Int)
	power := new(big.Int)
	quot := new(big.Int)
	tmp := new(big.Int)

	km1 := big.NewInt(int64(k - 1))
	kk := big.NewInt(int64(k))
	expKm1 := big.NewInt(int64(k - 1))
	expK := big.NewInt(int64(k))

	// Newton for x^k=n:
	// y=((k-1)x+n/x^(k-1))/k.
	// One step from either side is on/above the positive real root.
	power.Exp(x, expKm1, nil)
	quot.Quo(n, power)
	y.Mul(x, km1)
	y.Add(y, quot)
	y.Quo(y, kk)
	x.Set(y)

	for {
		power.Exp(x, expKm1, nil)
		quot.Quo(n, power)
		y.Mul(x, km1)
		y.Add(y, quot)
		y.Quo(y, kk)

		if y.Cmp(x) >= 0 {
			break
		}
		x.Set(y)
	}

	// Defensive exact correction. Normally zero steps.
	tmp.Add(x, one)
	power.Exp(tmp, expK, nil)
	if power.Cmp(n) <= 0 {
		x.Set(tmp)
	}

	power.Exp(x, expK, nil)
	for power.Cmp(n) > 0 {
		x.Sub(x, one)
		power.Exp(x, expK, nil)
	}
	return x
}

// v2Factorial is Legendre's formula specialized to p=2:
// v2(n!) = n - popcount(n).
func v2Factorial(n int) int {
	return n - bits.OnesCount(uint(n))
}

// v2Candidates returns a tiny, complete set of possible B values using the
// 2-adic valuation of A!.
//
// Let E=v2(A!). In a solution
//
//	A! = product_{j=1}^k (B+j),
//
// choose B+i with maximal 2-adic valuation m. If 2^m>k-1, it is unique and
// for j!=i
//
//	v2(B+j)=v2(j-i).
//
// Therefore
//
//	m = E - v2((i-1)!) - v2((k-i)!)
//
// and exact valuation m means
//
//	B == 2^m-i  (mod 2^(m+1)).
//
// Also, with x=(A!)^(1/k),
//
//	B < x < B+k.
//
// Since 2^(m+1)>k, only two adjacent residue blocks can contain B.
//
// The key speed trick is
//
//	floor(x/2^s) = floor(( floor(A!/2^(ks)) )^(1/k)).
//
// So we root a heavily right-shifted integer instead of all of A!.
// The bool is false only when the simple unique-max argument cannot be proved;
// the caller then falls back to a full exact kth root.
func v2Candidates(fact *big.Int, a, k int, stats *Stats) ([]*big.Int, bool) {
	E := v2Factorial(a)
	threshold := bits.Len(uint(k - 1)) // 2^threshold > k-1

	// Ensure the actual maximum valuation must reach threshold.
	if (E+k-1)/k < threshold {
		return nil, false
	}

	ms := make([]int, k)
	mMin := int(^uint(0) >> 1)

	for i := 1; i <= k; i++ {
		m := E - v2Factorial(i-1) - v2Factorial(k-i)
		if m < threshold {
			return nil, false
		}
		ms[i-1] = m
		if m < mMin {
			mMin = m
		}
	}

	// Compute qBase=floor(x/2^(mMin+1)) exactly using a much smaller root.
	baseShift := k * (mMin + 1)
	var qBase *big.Int
	if baseShift >= fact.BitLen() {
		qBase = new(big.Int)
	} else {
		scaled := new(big.Int).Rsh(new(big.Int).Set(fact), uint(baseShift))
		qBase = nthRootFloor(scaled, k)
		stats.reducedRoots.Add(1)
	}

	out := make([]*big.Int, 0, 2*k)

	for i := 1; i <= k; i++ {
		m := ms[i-1]

		// floor(x/2^(m+1)) follows by an exact additional right shift.
		q := new(big.Int).Rsh(
			new(big.Int).Set(qBase),
			uint(m-mMin),
		)

		modulus := new(big.Int).Lsh(big.NewInt(1), uint(m+1))
		residue := new(big.Int).Lsh(big.NewInt(1), uint(m))
		residue.Sub(residue, big.NewInt(int64(i)))

		// q*modulus + residue
		b := new(big.Int).Mul(new(big.Int).Set(q), modulus)
		b.Add(b, residue)
		out = append(out, b)

		// Previous residue block. Because B is less than x by <k and
		// modulus>k, no other block can contain a solution.
		if q.Sign() > 0 {
			out = append(out,
				new(big.Int).Sub(new(big.Int).Set(b), modulus),
			)
		}
	}

	return out, true
}

// fallbackCandidates is completely general. Since
//
//	B^k < product(B+i) = A! < (B+k)^k,
//
// B must lie among at most k+1 integers immediately below floor((A!)^(1/k)).
func fallbackCandidates(fact *big.Int, k int, stats *Stats) []*big.Int {
	stats.fallbacks.Add(1)
	q := nthRootFloor(fact, k)

	out := make([]*big.Int, 0, k+1)
	for d := 0; d <= k; d++ {
		b := new(big.Int).Sub(new(big.Int).Set(q), big.NewInt(int64(d)))
		if b.Sign() >= 0 {
			out = append(out, b)
		}
	}
	return out
}

// productCmp compares product_{i=1}^k(B+i) with A! exactly.
func productCmp(fact, b *big.Int, k int, stats *Stats) int {
	stats.comparisons.Add(1)

	p := big.NewInt(1)
	t := new(big.Int)
	for i := 1; i <= k; i++ {
		t.Set(b)
		t.Add(t, big.NewInt(int64(i)))
		p.Mul(p, t)

		// All remaining factors are positive, so we may stop immediately.
		if p.Cmp(fact) > 0 {
			return 1
		}
	}
	return p.Cmp(fact)
}

// Candidates are already a necessary superset. Because the product is strictly
// increasing in B, sort them and binary-search rather than evaluating all of
// them.
func findAmongCandidates(
	fact *big.Int,
	cands []*big.Int,
	k int,
	bmin, bmax *big.Int,
	a int,
	stats *Stats,
) *big.Int {
	aa := big.NewInt(int64(a))
	filtered := cands[:0]

	for _, b := range cands {
		if b.Sign() >= 0 &&
			b.Cmp(bmin) >= 0 &&
			b.Cmp(bmax) <= 0 &&
			b.Cmp(aa) >= 0 {
			filtered = append(filtered, b)
		}
	}
	if len(filtered) == 0 {
		return nil
	}

	sort.Slice(filtered, func(i, j int) bool {
		return filtered[i].Cmp(filtered[j]) < 0
	})

	// Deduplicate in place.
	u := 1
	for i := 1; i < len(filtered); i++ {
		if filtered[i].Cmp(filtered[u-1]) != 0 {
			filtered[u] = filtered[i]
			u++
		}
	}
	filtered = filtered[:u]

	lo, hi := 0, len(filtered)-1
	for lo <= hi {
		mid := (lo + hi) >> 1
		cmp := productCmp(fact, filtered[mid], k, stats)

		switch {
		case cmp < 0:
			lo = mid + 1
		case cmp > 0:
			hi = mid - 1
		default:
			return filtered[mid]
		}
	}
	return nil
}

// modularFilter stores, for every k, the A values that survive several exact
// necessary conditions.  If q > A, then an integer solution must satisfy
//
//	A! mod q = product_{i=1}^k (B+i) mod q.
//
// For each q and k we precompute the image of the rising-factorial polynomial
// over F_q as a bitset.  A miss is a rigorous rejection: no approximation is
// involved.
type modularFilter struct {
	pass  [][]uint64 // pass[k] is a bitset indexed by A, allocated up to hiA[k]
	hiA   []int      // largest A allowed for k by the exact decimal bounds
	count []uint64   // current number of surviving A values per k
}

func bitsetMake(n int, fill bool) []uint64 {
	w := make([]uint64, (n+64)>>6)
	if fill {
		for i := range w {
			w[i] = ^uint64(0)
		}
	}
	return w
}

func bitsetSet(b []uint64, i int)      { b[i>>6] |= uint64(1) << uint(i&63) }
func bitsetClear(b []uint64, i int)    { b[i>>6] &^= uint64(1) << uint(i&63) }
func bitsetHas(b []uint64, i int) bool { return (b[i>>6]>>(uint(i&63)))&1 != 0 }

// bitsetSetAtomic sets bit i under concurrent writers. The load-first check
// makes the common already-set case a plain read.
func bitsetSetAtomic(b []uint64, i uint64) {
	w := &b[i>>6]
	mask := uint64(1) << uint(i&63)
	if atomic.LoadUint64(w)&mask == 0 {
		atomic.OrUint64(w, mask)
	}
}

// bitsetFillRange sets bits lo..hi inclusive using whole-word writes.
func bitsetFillRange(b []uint64, lo, hi int) {
	if hi < lo {
		return
	}
	wLo, wHi := lo>>6, hi>>6
	maskLo := ^uint64(0) << uint(lo&63)
	maskHi := ^uint64(0) >> uint(63-(hi&63))
	if wLo == wHi {
		b[wLo] |= maskLo & maskHi
		return
	}
	b[wLo] |= maskLo
	for i := wLo + 1; i < wHi; i++ {
		b[i] = ^uint64(0)
	}
	b[wHi] |= maskHi
}

// mulMod64 computes a*b mod m exactly for arbitrary uint64 values with m>0.
func mulMod64(a, b, m uint64) uint64 {
	hi, lo := bits.Mul64(a, b)
	_, rem := bits.Div64(hi, lo, m)
	return rem
}

func powMod64(a, e, m uint64) uint64 {
	r := uint64(1)
	for e != 0 {
		if e&1 != 0 {
			r = mulMod64(r, a, m)
		}
		e >>= 1
		if e != 0 {
			a = mulMod64(a, a, m)
		}
	}
	return r
}

// Deterministic Miller-Rabin for uint64.
func isPrime64(n uint64) bool {
	if n < 2 {
		return false
	}
	for _, p := range []uint64{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37} {
		if n == p {
			return true
		}
		if n%p == 0 {
			return false
		}
	}
	d, s := n-1, 0
	for d&1 == 0 {
		d >>= 1
		s++
	}
	// This base set is deterministic for all 64-bit integers.
	for _, a := range []uint64{2, 325, 9375, 28178, 450775, 9780504, 1795265022} {
		if a%n == 0 {
			continue
		}
		x := powMod64(a%n, d, n)
		if x == 1 || x == n-1 {
			continue
		}
		witness := true
		for r := 1; r < s; r++ {
			x = mulMod64(x, x, n)
			if x == n-1 {
				witness = false
				break
			}
		}
		if witness {
			return false
		}
	}
	return true
}

func nextPrime64(n uint64) uint64 {
	if n <= 2 {
		return 2
	}
	if n&1 == 0 {
		n++
	}
	for !isPrime64(n) {
		n += 2
	}
	return n
}

// initializeSizePass applies only exact, integer-valued magnitude bounds.
//
// Write
//
//	L(A) = sum_{n=1}^A floor(log10 n),
//	U(A) = sum_{n=1}^A ceil (log10 n).
//
// Then, without any floating point,
//
//	10^L(A) <= A! <= 10^U(A).
//
// If B >= 10^minDigits, then
//
//	product(B+i) > 10^(minDigits*k).
//
// Hence U(A) <= minDigits*k is impossible.
//
// If B <= 10^maxDigits, then (for our huge ranges and tiny k)
//
//	product(B+i) < (2*10^maxDigits)^k < 10^((maxDigits+1)k).
//
// Hence L(A) >= (maxDigits+1)k is impossible.
//
// Because L and U are nondecreasing in A, the surviving A values of each k
// form one contiguous interval.  A single sweep records both endpoints and
// the intervals are then filled with whole-word bitset writes.
func initializeSizePass(amax, kmax, minDigits, maxDigits int) (*modularFilter, uint64) {
	aLo := make([]int, kmax+1) // first a with U(a) > minDigits*k; 0 = never
	aHi := make([]int, kmax+1) // last a with L(a) < (maxDigits+1)*k

	var lowerExp, upperExp int64
	floorD := int64(0) // floor(log10 a)
	ceilD := int64(0)  // ceil(log10 a)
	nextFloor := int64(10)
	nextCeil := int64(2) // ceil(log10 a)=1 starts at a=2
	ceilUpper := int64(10)

	kReady := 1 // largest k whose lower endpoint has been recorded
	kDead := 1  // largest k whose upper condition has already failed

	for a := 2; a <= amax; a++ {
		aa := int64(a)
		if aa == nextFloor {
			floorD++
			if nextFloor <= int64(amax)/10 {
				nextFloor *= 10
			} else {
				nextFloor = int64(amax) + 1
			}
		}
		if aa == nextCeil {
			if aa == 2 {
				ceilD = 1
			} else {
				ceilD++
			}
			if ceilUpper < int64(amax) {
				nextCeil = ceilUpper + 1
				if ceilUpper <= int64(amax)/10 {
					ceilUpper *= 10
				} else {
					ceilUpper = int64(amax)
				}
			} else {
				nextCeil = int64(amax) + 1
			}
		}

		lowerExp += floorD
		upperExp += ceilD

		// U > minDigits*k  <=>  k <= (U-1)/minDigits.
		maxK := int((upperExp - 1) / int64(minDigits))
		if maxK > kmax {
			maxK = kmax
		}
		for kReady < maxK {
			kReady++
			aLo[kReady] = a
		}

		// L < (maxDigits+1)*k fails exactly for k <= L/(maxDigits+1).
		failK := int(lowerExp / int64(maxDigits+1))
		if failK > kmax {
			failK = kmax
		}
		for kDead < failK {
			kDead++
			aHi[kDead] = a - 1
		}
	}
	for k := kDead + 1; k <= kmax; k++ {
		aHi[k] = amax
	}

	mf := &modularFilter{
		pass:  make([][]uint64, kmax+1),
		hiA:   make([]int, kmax+1),
		count: make([]uint64, kmax+1),
	}
	var pairs uint64
	for k := 2; k <= kmax; k++ {
		lo, hi := aLo[k], aHi[k]
		if lo < 2 {
			continue // never became possible
		}
		if hi > amax {
			hi = amax
		}
		if hi < lo {
			continue
		}
		mf.hiA[k] = hi
		mf.pass[k] = bitsetMake(hi+1, false)
		bitsetFillRange(mf.pass[k], lo, hi)
		mf.count[k] = uint64(hi - lo + 1)
		pairs += mf.count[k]
	}
	return mf, pairs
}

// totalPass returns the incrementally maintained survivor count.
func totalPass(mf *modularFilter, kmax int) uint64 {
	var n uint64
	for k := 2; k <= kmax; k++ {
		n += mf.count[k]
	}
	return n
}

// countPassExact recounts survivors from the bitsets; used by the self-test to
// validate the incremental counters.
func countPassExact(mf *modularFilter, kmax int) uint64 {
	var n uint64
	for k := 2; k <= kmax; k++ {
		for _, w := range mf.pass[k] {
			n += uint64(bits.OnesCount64(w))
		}
	}
	return n
}

type akPair struct {
	a int
	k int
}

func addMod64(a, b, q uint64) uint64 {
	a += b
	if a >= q {
		a -= q
	}
	return a
}

func subMod64(a, b, q uint64) uint64 {
	if a >= b {
		return a - b
	}
	return q - (b - a)
}

func mulModFast(a, b, q uint64) uint64 {
	if q <= uint64(^uint32(0)) {
		return (a * b) % q
	}
	return mulMod64(a, b, q)
}

func polyTrim(p []uint64) []uint64 {
	if len(p) == 0 {
		return []uint64{0}
	}
	i := len(p) - 1
	for i > 0 && p[i] == 0 {
		i--
	}
	return p[:i+1]
}

// polyMulRem multiplies two polynomials of degree < deg(modPoly) and reduces
// modulo the monic polynomial modPoly. Coefficients are in F_q.
func polyMulRem(a, b, modPoly []uint64, q uint64) []uint64 {
	k := len(modPoly) - 1
	if k <= 0 {
		return []uint64{0}
	}
	tmp := make([]uint64, 2*k-1)
	for i, ai := range a {
		if ai == 0 {
			continue
		}
		for j, bj := range b {
			if bj == 0 {
				continue
			}
			v := mulModFast(ai, bj, q)
			tmp[i+j] = addMod64(tmp[i+j], v, q)
		}
	}
	for d := len(tmp) - 1; d >= k; d-- {
		c := tmp[d]
		if c == 0 {
			continue
		}
		for j := 0; j < k; j++ {
			v := mulModFast(c, modPoly[j], q)
			tmp[d-k+j] = subMod64(tmp[d-k+j], v, q)
		}
	}
	out := make([]uint64, k)
	copy(out, tmp[:k])
	return polyTrim(out)
}

// polyRemGeneral returns a mod b over F_q.
func polyRemGeneral(a, b []uint64, q uint64) []uint64 {
	a = append([]uint64(nil), a...)
	a = polyTrim(a)
	b = polyTrim(b)
	db := len(b) - 1
	if db == 0 && b[0] == 0 {
		panic("polynomial division by zero")
	}
	if len(a)-1 < db {
		return a
	}
	invLead := powMod64(b[db], q-2, q)
	for len(a)-1 >= db && !(len(a) == 1 && a[0] == 0) {
		da := len(a) - 1
		c := mulModFast(a[da], invLead, q)
		shift := da - db
		if c != 0 {
			for j := 0; j <= db; j++ {
				v := mulModFast(c, b[j], q)
				a[shift+j] = subMod64(a[shift+j], v, q)
			}
		}
		a = polyTrim(a)
	}
	return a
}

func polyGCDDegree(a, b []uint64, q uint64) int {
	a = polyTrim(append([]uint64(nil), a...))
	b = polyTrim(append([]uint64(nil), b...))
	for !(len(b) == 1 && b[0] == 0) {
		r := polyRemGeneral(a, b, q)
		a, b = b, r
	}
	return len(polyTrim(a)) - 1
}

// risingPoly returns f_k(X)=prod_{i=1}^k (X+i) over F_q.
func risingPoly(k int, q uint64) []uint64 {
	p := []uint64{1}
	for i := 1; i <= k; i++ {
		n := make([]uint64, len(p)+1)
		ii := uint64(i) % q
		for j, c := range p {
			n[j] = addMod64(n[j], mulModFast(c, ii, q), q)
			n[j+1] = addMod64(n[j+1], c, q)
		}
		p = n
	}
	return p
}

// hasRisingFactorialRoot tests exactly whether
//
//	prod_{i=1}^k (x+i) == y (mod q)
//
// has x in F_q. Let F(X)=prod(X+i)-y. F has an F_q root iff
//
//	gcd(F, X^q-X)
//
// has positive degree.
func hasRisingFactorialRoot(base []uint64, y, q uint64) bool {
	F := append([]uint64(nil), base...)
	F[0] = subMod64(F[0], y, q)
	F = polyTrim(F)
	k := len(F) - 1
	if k <= 0 {
		return F[0] == 0
	}

	res := []uint64{1}
	x := []uint64{0, 1}
	if k == 1 {
		x = polyRemGeneral(x, F, q)
	}
	e := q
	for e != 0 {
		if e&1 != 0 {
			res = polyMulRem(res, x, F, q)
		}
		e >>= 1
		if e != 0 {
			x = polyMulRem(x, x, F, q)
		}
	}

	xred := []uint64{0, 1}
	if k == 1 {
		xred = polyRemGeneral(xred, F, q)
	}
	hLen := len(res)
	if len(xred) > hLen {
		hLen = len(xred)
	}
	h := make([]uint64, hLen)
	copy(h, res)
	for i, v := range xred {
		h[i] = subMod64(h[i], v, q)
	}
	h = polyTrim(h)

	return polyGCDDegree(F, h, q) > 0
}

func collectPairs(mf *modularFilter, kmax int) ([]akPair, []int) {
	pairs := make([]akPair, 0)
	seenA := make(map[int]struct{})
	for k := 2; k <= kmax; k++ {
		for wi, w := range mf.pass[k] {
			for w != 0 {
				bit := bits.TrailingZeros64(w)
				a := (wi << 6) + bit
				if a >= 2 {
					pairs = append(pairs, akPair{a: a, k: k})
					seenA[a] = struct{}{}
				}
				w &= w - 1
			}
		}
	}
	as := make([]int, 0, len(seenA))
	for a := range seenA {
		as = append(as, a)
	}
	sort.Ints(as)
	return pairs, as
}

// alignedSegments splits [2..maxA] into up to parts segments whose interior
// boundaries all have the form 64m-1.  Segment i covers (b[i-1], b[i]], so two
// segments never touch the same bitset word and workers cannot race on clears.
func alignedSegments(maxA, parts int) []int {
	b := []int{1}
	if maxA < 2 {
		return b
	}
	if parts < 1 {
		parts = 1
	}
	totalWords := maxA>>6 + 1
	if parts > totalWords {
		parts = totalWords
	}
	for i := 1; i < parts; i++ {
		w := totalWords * i / parts
		e := w<<6 - 1
		if e <= b[len(b)-1] {
			continue
		}
		if e >= maxA {
			break
		}
		b = append(b, e)
	}
	b = append(b, maxA)
	return b
}

// segmentedFactorialPrefix returns, for the segment boundaries b, the residues
//
//	factAt[i] = b[i]! mod q,
//
// computing the per-segment products in parallel and combining them with one
// short sequential prefix pass. Since b[0]=1, factAt[0]=1.
func segmentedFactorialPrefix(b []int, q uint64, workers int) []uint64 {
	n := len(b) - 1
	segProd := make([]uint64, n+1)
	var wg sync.WaitGroup
	idx := make(chan int)
	if workers < 1 {
		workers = 1
	}
	for w := 0; w < workers; w++ {
		wg.Add(1)
		go func() {
			defer wg.Done()
			for i := range idx {
				p := uint64(1)
				for a := b[i-1] + 1; a <= b[i]; a++ {
					p = mulModFast(p, uint64(a), q)
				}
				segProd[i] = p
			}
		}()
	}
	for i := 1; i <= n; i++ {
		idx <- i
	}
	close(idx)
	wg.Wait()

	factAt := make([]uint64, n+1)
	factAt[0] = 1 % q
	for i := 1; i <= n; i++ {
		factAt[i] = mulModFast(factAt[i-1], segProd[i], q)
	}
	return factAt
}

// buildGroupImages enumerates the exact image bitsets of f_k over F_q for all
// k in ks. All workers share one bitset per k and publish bits atomically, so
// parallelism costs no extra memory.
func buildGroupImages(q uint64, ks []int, workers int) [][]uint64 {
	kTop := ks[len(ks)-1]
	inGroup := make([]bool, kTop+1)
	for _, k := range ks {
		inGroup[k] = true
	}
	words := int((q + 64) >> 6)
	images := make([][]uint64, kTop+1)
	for _, k := range ks {
		images[k] = make([]uint64, words)
	}
	if workers < 1 {
		workers = 1
	}
	if uint64(workers) > q {
		workers = int(q)
	}

	var wg sync.WaitGroup
	for w := 0; w < workers; w++ {
		lo := q * uint64(w) / uint64(workers)
		hi := q * uint64(w+1) / uint64(workers)
		wg.Add(1)
		go func(lo, hi uint64) {
			defer wg.Done()
			for x := lo; x < hi; x++ {
				f := x + 1
				if f == q {
					f = 0
				}
				prod := uint64(1)
				for k := 1; k <= kTop; k++ {
					prod = mulModFast(prod, f, q)
					f++
					if f == q {
						f = 0
					}
					if k >= 2 && inGroup[k] {
						bitsetSetAtomic(images[k], prod)
					}
				}
			}
		}(lo, hi)
	}
	wg.Wait()
	return images
}

// applyDensePrimeFilter intersects the candidate sets of every k in ks with
// the exact image condition modulo q. The survivor sweep runs in parallel
// over 64-aligned A segments: a segmented prefix product supplies A! mod q at
// each boundary, and inside a segment only surviving bits are inspected.
func applyDensePrimeFilter(mf *modularFilter, ks []int, q uint64, workers int) {
	if len(ks) == 0 {
		return
	}
	images := buildGroupImages(q, ks, workers)

	maxHiA := 0
	for _, k := range ks {
		if mf.hiA[k] > maxHiA {
			maxHiA = mf.hiA[k]
		}
	}
	if maxHiA < 2 {
		return
	}
	if workers < 1 {
		workers = 1
	}
	b := alignedSegments(maxHiA, workers*4)
	factAt := segmentedFactorialPrefix(b, q, workers)

	kTop := ks[len(ks)-1]
	cleared := make([]uint64, kTop+1)
	var mu sync.Mutex
	const blockLen = 1 << 16

	idx := make(chan int)
	var wg sync.WaitGroup
	for w := 0; w < workers; w++ {
		wg.Add(1)
		go func() {
			defer wg.Done()
			buf := make([]uint64, blockLen)
			local := make([]uint64, kTop+1)
			for i := range idx {
				lo, hi := b[i-1], b[i]
				fact := factAt[i-1]
				for bs := lo + 1; bs <= hi; bs += blockLen {
					be := bs + blockLen - 1
					if be > hi {
						be = hi
					}
					for a := bs; a <= be; a++ {
						fact = mulModFast(fact, uint64(a), q)
						buf[a-bs] = fact
					}
					for _, k := range ks {
						hiK := mf.hiA[k]
						if bs > hiK {
							continue
						}
						e := be
						if e > hiK {
							e = hiK
						}
						bm := mf.pass[k]
						img := images[k]
						for wi := bs >> 6; wi <= e>>6; wi++ {
							w := bm[wi]
							if w == 0 {
								continue
							}
							base := wi << 6
							for w != 0 {
								bit := bits.TrailingZeros64(w)
								w &= w - 1
								a := base + bit
								if a < bs || a > e {
									continue
								}
								if !bitsetHas(img, int(buf[a-bs])) {
									bm[wi] &^= uint64(1) << uint(bit)
									local[k]++
								}
							}
						}
					}
				}
			}
			mu.Lock()
			for k := range local {
				cleared[k] += local[k]
			}
			mu.Unlock()
		}()
	}
	for i := 1; i < len(b); i++ {
		idx <- i
	}
	close(idx)
	wg.Wait()

	for _, k := range ks {
		mf.count[k] -= cleared[k]
	}
}

// factorialResiduesAt computes A! mod q only for requested A values.
// It chooses the cheaper of a forward sweep from 1 or a reverse Wilson sweep
// from q-1. The forward sweep parallelizes through a segmented prefix product;
// the reverse path uses one batch inversion for all requested tails.
func factorialResiduesAt(as []int, q uint64, workers int) map[int]uint64 {
	out := make(map[int]uint64, len(as))
	if len(as) == 0 {
		return out
	}
	minA, maxA := as[0], as[len(as)-1]
	forwardCost := uint64(maxA)
	reverseCost := q - 1 - uint64(minA)

	if forwardCost <= reverseCost {
		b := alignedSegments(maxA, workers*4)
		factAt := segmentedFactorialPrefix(b, q, workers)

		res := make([]uint64, len(as))
		idx := make(chan int)
		var wg sync.WaitGroup
		if workers < 1 {
			workers = 1
		}
		for w := 0; w < workers; w++ {
			wg.Add(1)
			go func() {
				defer wg.Done()
				for i := range idx {
					lo, hi := b[i-1], b[i]
					// Requested indices inside (lo, hi].
					first := sort.SearchInts(as, lo+1)
					last := sort.SearchInts(as, hi+1)
					if first == last {
						continue
					}
					fact := factAt[i-1]
					p := first
					for a := lo + 1; a <= hi && p < last; a++ {
						fact = mulModFast(fact, uint64(a), q)
						for p < last && as[p] == a {
							res[p] = fact
							p++
						}
					}
				}
			}()
		}
		for i := 1; i < len(b); i++ {
			idx <- i
		}
		close(idx)
		wg.Wait()

		for i, a := range as {
			out[a] = res[i]
		}
		return out
	}

	tails := make([]uint64, len(as))
	idx := len(as) - 1
	tail := uint64(1)
	for a := int(q - 1); a >= minA; a-- {
		for idx >= 0 && as[idx] == a {
			tails[idx] = tail
			idx--
		}
		if a > minA {
			tail = mulModFast(tail, uint64(a), q)
		}
	}

	prefix := make([]uint64, len(tails)+1)
	prefix[0] = 1
	for i, t := range tails {
		prefix[i+1] = mulModFast(prefix[i], t, q)
	}
	invAll := powMod64(prefix[len(tails)], q-2, q)
	invs := make([]uint64, len(tails))
	for i := len(tails) - 1; i >= 0; i-- {
		invs[i] = mulModFast(invAll, prefix[i], q)
		invAll = mulModFast(invAll, tails[i], q)
	}
	for i, a := range as {
		out[a] = q - invs[i]
		if out[a] == q {
			out[a] = 0
		}
	}
	return out
}

// applySparsePrimeFilter applies the same exact finite-field necessary
// condition as the image-bitset filter, but candidate-by-candidate via a
// polynomial root test. Its cost depends on survivor count and k*log(q), not q.
func applySparsePrimeFilter(mf *modularFilter, kmax int, q uint64, workers int) uint64 {
	pairs, as := collectPairs(mf, kmax)
	if len(pairs) == 0 {
		return 0
	}
	residues := factorialResiduesAt(as, q, workers)

	base := make(map[int][]uint64)
	for _, p := range pairs {
		if _, ok := base[p.k]; !ok {
			base[p.k] = risingPoly(p.k, q)
		}
	}

	if workers < 1 {
		workers = 1
	}
	if workers > len(pairs) {
		workers = len(pairs)
	}
	reject := make([]bool, len(pairs))
	var wg sync.WaitGroup
	for w := 0; w < workers; w++ {
		lo := len(pairs) * w / workers
		hi := len(pairs) * (w + 1) / workers
		wg.Add(1)
		go func(lo, hi int) {
			defer wg.Done()
			for idx := lo; idx < hi; idx++ {
				p := pairs[idx]
				if !hasRisingFactorialRoot(base[p.k], residues[p.a], q) {
					reject[idx] = true
				}
			}
		}(lo, hi)
	}
	wg.Wait()

	// Clear sequentially: different pairs can share the same uint64 word in a
	// bitset, so concurrent read-modify-write clears would otherwise race.
	for i, bad := range reject {
		if bad {
			p := pairs[i]
			bitsetClear(mf.pass[p.k], p.a)
			mf.count[p.k]--
		}
	}
	return totalPass(mf, kmax)
}

// buildModularFilter starts from the exact decimal-magnitude filter above and
// intersects it with exact finite-field image conditions.  It stops as soon as
// no (A,k) pair remains, so increasing -filter-primes is only a maximum.
// Active k values are processed in groups whose image bitsets fit into the
// configured memory budget; grouping changes scheduling only, never the
// resulting survivor set.
func buildModularFilter(
	amax, kmax, minDigits, maxDigits, nprimes int,
	sparseThreshold uint64,
	filterWorkers int,
	filterMemoryMB int,
	progress bool,
) (*modularFilter, []uint64, uint64) {
	mf, initialPairs := initializeSizePass(amax, kmax, minDigits, maxDigits)
	if initialPairs == 0 || nprimes == 0 {
		return mf, nil, initialPairs
	}

	primes := make([]uint64, 0, nprimes)
	q := nextPrime64(uint64(amax + 1))

	for round := 0; round < nprimes; round++ {
		leftBefore := totalPass(mf, kmax)
		if leftBefore == 0 {
			break
		}
		primes = append(primes, q)

		if sparseThreshold > 0 && leftBefore <= sparseThreshold {
			left := applySparsePrimeFilter(mf, kmax, q, filterWorkers)
			if progress {
				fmt.Printf(
					"filter %2d q=%d [sparse gcd] -> %d surviving (A,k) pairs\n",
					round+1, q, left,
				)
			}
			if left == 0 {
				break
			}
			q = nextPrime64(q + 100)
			continue
		}

		activeKs := make([]int, 0, kmax-1)
		for k := 2; k <= kmax; k++ {
			if mf.count[k] > 0 {
				activeKs = append(activeKs, k)
			}
		}
		if len(activeKs) == 0 {
			break
		}

		groupCap := len(activeKs)
		if filterMemoryMB > 0 {
			imageBytes := (uint64((q + 64) >> 6)) * 8
			byMem := int((uint64(filterMemoryMB) << 20) / imageBytes)
			if byMem < 1 {
				byMem = 1
			}
			if groupCap > byMem {
				groupCap = byMem
			}
		}
		groups := 0
		for g := 0; g < len(activeKs); g += groupCap {
			end := g + groupCap
			if end > len(activeKs) {
				end = len(activeKs)
			}
			applyDensePrimeFilter(mf, activeKs[g:end], q, filterWorkers)
			groups++
		}

		left := totalPass(mf, kmax)
		if progress {
			fmt.Printf(
				"filter %2d q=%d [dense image x%d g%d] -> %d surviving (A,k) pairs\n",
				round+1, q, filterWorkers, groups, left,
			)
		}
		if left == 0 {
			break
		}
		q = nextPrime64(q + 100)
	}
	return mf, primes, initialPairs
}

type searchTask struct {
	a    int
	fact *big.Int
	ks   []int
}

func buildProductBounds(kmax int, bmin, bmax *big.Int) ([]*big.Int, []*big.Int) {
	pmin := make([]*big.Int, kmax+1)
	pmax := make([]*big.Int, kmax+1)
	lo := big.NewInt(1)
	hi := big.NewInt(1)
	t := new(big.Int)
	for k := 1; k <= kmax; k++ {
		t.Add(bmin, big.NewInt(int64(k)))
		lo.Mul(lo, t)
		t.Add(bmax, big.NewInt(int64(k)))
		hi.Mul(hi, t)
		if k >= 2 {
			pmin[k] = new(big.Int).Set(lo)
			pmax[k] = new(big.Int).Set(hi)
		}
	}
	return pmin, pmax
}

func workerSearch(
	tasks <-chan searchTask,
	bmin, bmax *big.Int,
	stats *Stats,
	out chan<- Solution,
	wg *sync.WaitGroup,
) {
	defer wg.Done()
	for task := range tasks {
		for _, k := range task.ks {
			cands, ok := v2Candidates(task.fact, task.a, k, stats)
			if !ok {
				cands = fallbackCandidates(task.fact, k, stats)
			}
			if b := findAmongCandidates(task.fact, cands, k, bmin, bmax, task.a, stats); b != nil {
				out <- Solution{A: task.a, B: new(big.Int).Set(b), K: k}
			}
		}
	}
}

// searchFiltered performs only one exact factorial sweep for all k values.
func searchFiltered(
	amax, kmax, workers int,
	bmin, bmax *big.Int,
	filter *modularFilter,
	stats *Stats,
	out chan<- Solution,
) (filterPairs, filterAs, boundedPairs uint64) {
	pmin, pmax := buildProductBounds(kmax, bmin, bmax)

	// Collect only (A,k) pairs surviving every modular filter.
	pairsList, as := collectPairs(filter, kmax)
	byA := make(map[int][]int, len(as))
	for _, p := range pairsList {
		byA[p.a] = append(byA[p.a], p.k)
		filterPairs++
	}
	filterAs = uint64(len(as))

	tasks := make(chan searchTask, workers*2)
	var wg sync.WaitGroup
	for i := 0; i < workers; i++ {
		wg.Add(1)
		go workerSearch(tasks, bmin, bmax, stats, out, &wg)
	}

	// Compute A! only once.  Gaps between surviving A values are multiplied by
	// a balanced product tree via MulRange, which is much faster than doing all
	// small multiplications separately for each k.
	fact := big.NewInt(1)
	chunk := new(big.Int)
	prev := 1
	for _, a := range as {
		if a > prev {
			chunk.MulRange(int64(prev+1), int64(a))
			fact.Mul(fact, chunk)
			prev = a
		}

		ks := make([]int, 0, len(byA[a]))
		for _, k := range byA[a] {
			if fact.Cmp(pmin[k]) < 0 || fact.Cmp(pmax[k]) > 0 {
				continue
			}
			boundedPairs++
			stats.pairs.Add(1)
			ks = append(ks, k)
		}
		if len(ks) == 0 {
			continue
		}

		// The worker may outlive the next factorial update, so give it one
		// immutable copy per A (shared by all surviving k for that A).
		tasks <- searchTask{
			a:    a,
			fact: new(big.Int).Set(fact),
			ks:   ks,
		}
	}
	close(tasks)
	wg.Wait()
	return
}

func selfTest() error {
	// The zero polynomial has a canonical representation throughout the GCD
	// implementation, including when a caller supplies an empty coefficient list.
	if degree := polyGCDDegree([]uint64{1}, nil, 29); degree != 0 {
		return fmt.Errorf("zero-polynomial GCD check failed: degree=%d", degree)
	}

	// Known nontrivial solution.
	f6 := new(big.Int).MulRange(1, 6)
	f7 := new(big.Int).MulRange(1, 7)
	f10 := new(big.Int).MulRange(1, 10)
	lhs := new(big.Int).Mul(new(big.Int).Set(f6), f7)
	if lhs.Cmp(f10) != 0 {
		return fmt.Errorf("known solution check failed")
	}

	// Verify the exact decimal envelopes on small factorials.
	fact := big.NewInt(1)
	pow10L := new(big.Int)
	pow10U := new(big.Int)
	for a := 2; a <= 500; a++ {
		fact.Mul(fact, big.NewInt(int64(a)))
		L := decFloorLogFactSlow(a)
		U := decCeilLogFactSlow(a)
		pow10L.Exp(big.NewInt(10), big.NewInt(L), nil)
		pow10U.Exp(big.NewInt(10), big.NewInt(U), nil)
		if fact.Cmp(pow10L) < 0 || fact.Cmp(pow10U) > 0 {
			return fmt.Errorf("decimal envelope failed at A=%d", a)
		}
	}

	// The interval form of the decimal size filter must agree bit-for-bit
	// with the direct per-A evaluation of U(A) > d*k and L(A) < (D+1)*k.
	amaxT, kmaxT, dT, DT := 4000, 9, 5, 40
	mfT, pairsT := initializeSizePass(amaxT, kmaxT, dT, DT)
	var refPairs uint64
	for a := 2; a <= amaxT; a++ {
		L := decFloorLogFactSlow(a)
		U := decCeilLogFactSlow(a)
		for k := 2; k <= kmaxT; k++ {
			want := U > int64(dT*k) && L < int64((DT+1)*k)
			got := mfT.pass[k] != nil && a <= mfT.hiA[k] && bitsetHas(mfT.pass[k], a)
			if got != want {
				return fmt.Errorf("size filter mismatch at A=%d k=%d: got=%v want=%v", a, k, got, want)
			}
			if want {
				refPairs++
			}
		}
	}
	if pairsT != refPairs || countPassExact(mfT, kmaxT) != refPairs || totalPass(mfT, kmaxT) != refPairs {
		return fmt.Errorf("size filter pair count mismatch: %d vs %d", pairsT, refPairs)
	}

	// The parallel shared-image dense filter, applied in two groups, must
	// agree with a direct enumeration of every image and every residue.
	qT := nextPrime64(uint64(amaxT + 1))
	activeT := make([]int, 0)
	for k := 2; k <= kmaxT; k++ {
		if mfT.count[k] > 0 {
			activeT = append(activeT, k)
		}
	}
	mid := len(activeT) / 2
	if mid < 1 {
		mid = 1
	}
	applyDensePrimeFilter(mfT, activeT[:mid], qT, 4)
	applyDensePrimeFilter(mfT, activeT[mid:], qT, 4)
	for _, k := range activeT {
		image := make([]bool, qT)
		for x := uint64(0); x < qT; x++ {
			prod := uint64(1)
			for i := 1; i <= k; i++ {
				prod = mulModFast(prod, (x+uint64(i))%qT, qT)
			}
			image[prod] = true
		}
		r := uint64(1)
		for a := 2; a <= amaxT; a++ {
			r = mulModFast(r, uint64(a), qT)
			L := decFloorLogFactSlow(a)
			U := decCeilLogFactSlow(a)
			inSize := U > int64(dT*k) && L < int64((DT+1)*k)
			want := inSize && image[r]
			got := mfT.pass[k] != nil && a <= mfT.hiA[k] && bitsetHas(mfT.pass[k], a)
			if got != want {
				return fmt.Errorf("dense filter mismatch at A=%d k=%d q=%d: got=%v want=%v", a, k, qT, got, want)
			}
		}
	}
	if countPassExact(mfT, kmaxT) != totalPass(mfT, kmaxT) {
		return fmt.Errorf("dense filter survivor counters diverged")
	}

	// The parallel segmented factorial residues must match a direct sweep.
	{
		q := uint64(10007)
		asT := []int{2, 3, 5, 100, 4093, 7001, 9990}
		got := factorialResiduesAt(asT, q, 3)
		r := uint64(1)
		want := make(map[int]uint64)
		for a := 2; a <= 9990; a++ {
			r = mulModFast(r, uint64(a), q)
			want[a] = r
		}
		for _, a := range asT {
			if got[a] != want[a] {
				return fmt.Errorf("factorial residue mismatch at A=%d: got=%d want=%d", a, got[a], want[a])
			}
		}
		// Force the reverse Wilson path and compare as well.
		asW := []int{9980, 9990, 10001}
		gotW := factorialResiduesAt(asW, q, 3)
		r = uint64(1)
		for a := 2; a <= 10001; a++ {
			r = mulModFast(r, uint64(a), q)
			if a == 9980 || a == 9990 || a == 10001 {
				if gotW[a] != r {
					return fmt.Errorf("wilson residue mismatch at A=%d: got=%d want=%d", a, gotW[a], r)
				}
			}
		}
	}

	// Verify the finite-field necessary condition on 6! = 8*9*10.
	const q uint64 = 101
	r := uint64(1)
	for a := uint64(2); a <= 6; a++ {
		r = (r * a) % q
	}
	p := uint64(1)
	for i := uint64(1); i <= 3; i++ {
		p = (p * ((7 + i) % q)) % q
	}
	if r != p {
		return fmt.Errorf("finite-field check failed")
	}

	// Exhaustively validate the sparse polynomial-GCD image test against
	// direct enumeration over several small prime fields.
	for _, qq := range []uint64{29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101} {
		for k := 2; k <= 8; k++ {
			base := risingPoly(k, qq)
			image := make([]bool, qq)
			for x := uint64(0); x < qq; x++ {
				prod := uint64(1)
				for i := 1; i <= k; i++ {
					prod = mulModFast(prod, (x+uint64(i))%qq, qq)
				}
				image[prod] = true
			}
			for y := uint64(0); y < qq; y++ {
				got := hasRisingFactorialRoot(base, y, qq)
				if got != image[y] {
					return fmt.Errorf("sparse GCD test failed: q=%d k=%d y=%d got=%v want=%v", qq, k, y, got, image[y])
				}
			}
		}
	}

	// Exact integer root sanity checks.
	for n := int64(1); n <= 1000; n++ {
		N := big.NewInt(n)
		for k := 2; k <= 7; k++ {
			rt := nthRootFloor(N, k)
			pk := new(big.Int).Exp(new(big.Int).Set(rt), big.NewInt(int64(k)), nil)
			next := new(big.Int).Add(new(big.Int).Set(rt), one)
			pkNext := new(big.Int).Exp(next, big.NewInt(int64(k)), nil)
			if pk.Cmp(N) > 0 || pkNext.Cmp(N) <= 0 {
				return fmt.Errorf("nthRootFloor failed for n=%d k=%d", n, k)
			}
		}
	}
	return nil
}

func decFloorLogFactSlow(a int) int64 {
	var s int64
	for n := 1; n <= a; n++ {
		x := n
		for x >= 10 {
			s++
			x /= 10
		}
	}
	return s
}

func decCeilLogFactSlow(a int) int64 {
	var s int64
	for n := 2; n <= a; n++ {
		d := int64(1)
		p := 10
		for n > p {
			d++
			p *= 10
		}
		s += d
	}
	return s
}

func main() {
	digits := flag.Int(
		"digits", 1000000,
		"upper exponent D in B <= 10^D",
	)
	minDigits := flag.Int(
		"min-digits", 3000,
		"lower exponent d in 10^d <= B; keep 3000 for a standalone extension of Habsieger",
	)
	workers := flag.Int(
		"workers", runtime.NumCPU(),
		"parallel workers for filters and the exact search",
	)
	filterPrimes := flag.Int(
		"filter-primes", 64,
		"maximum number of exact finite-field filters (stops early at zero survivors)",
	)
	sparseThreshold := flag.Uint64(
		"sparse-threshold", 2000000,
		"switch from q-sized image bitsets to exact polynomial-GCD filters at or below this survivor count",
	)
	filterMemoryMB := flag.Int(
		"filter-memory-mb", 8192,
		"memory budget for dense image bitsets; active k values are filtered in groups that fit",
	)
	progress := flag.Bool(
		"progress", false,
		"print survivor count after each exact modular filter",
	)
	selfTestFlag := flag.Bool(
		"self-test", false,
		"run internal exactness/regression checks and exit",
	)
	flag.Parse()

	if *selfTestFlag {
		if err := selfTest(); err != nil {
			panic(err)
		}
		fmt.Println("self-test: OK")
		return
	}

	if *minDigits < 3000 {
		panic("-min-digits must be >= 3000")
	}
	if *digits < *minDigits {
		panic("-digits must be >= -min-digits")
	}
	if *workers < 1 {
		panic("-workers must be >= 1")
	}

	start := time.Now()

	amax, kmax := theoremBounds(*digits)

	fmt.Printf("searching hypothetical solutions other than (6,7,10)\n")
	fmt.Printf("B range : 10^%d .. 10^%d\n", *minDigits, *digits)
	fmt.Printf("bounds  : A <= %d, 2 <= k=C-B <= %d\n", amax, kmax)
	fmt.Printf("workers : %d\n", *workers)
	fmt.Printf("filters : up to %d exact prime-field filters\n", *filterPrimes)
	fmt.Printf("hybrid  : sparse polynomial-GCD mode at <= %d survivors\n", *sparseThreshold)
	fmt.Printf("memory  : dense image bitsets grouped within %d MiB\n\n", *filterMemoryMB)

	filterStart := time.Now()
	filter, filterQs, initialPairs := buildModularFilter(
		amax, kmax, *minDigits, *digits, *filterPrimes, *sparseThreshold, *workers, *filterMemoryMB, *progress,
	)
	fmt.Printf("size-bound candidates   : %d pairs\n", initialPairs)
	fmt.Printf("filters actually used  : %d\n", len(filterQs))
	if len(filterQs) <= 16 {
		fmt.Printf("filter primes          : %v\n", filterQs)
	} else {
		fmt.Printf("filter primes          : %v ... %v\n", filterQs[:8], filterQs[len(filterQs)-8:])
	}
	fmt.Printf("filter build           : %v\n\n", time.Since(filterStart))

	remaining := totalPass(filter, kmax)
	if remaining == 0 {
		fmt.Printf("modular survivors      : 0 pairs\n")
		fmt.Printf("inside B-range         : 0 pairs\n")
		fmt.Printf("known nontrivial solution: 6! * 7! = 10!\n")
		fmt.Printf("new solutions found    : 0\n")
		fmt.Printf("(A,k) pairs examined   : 0\n")
		fmt.Printf("reduced kth roots      : 0\n")
		fmt.Printf("full-root fallbacks    : 0\n")
		fmt.Printf("exact product compares : 0\n")
		fmt.Printf("elapsed                : %v\n", time.Since(start))
		return
	}

	// The endpoint powers of ten are only needed once modular survivors exist;
	// building 10^digits eagerly would cost minutes of big-integer time for a
	// number that the zero-survivor case never uses.
	bmin := new(big.Int).Exp(big.NewInt(10), big.NewInt(int64(*minDigits)), nil)
	bmax := new(big.Int).Exp(
		big.NewInt(10),
		big.NewInt(int64(*digits)),
		nil,
	)

	stats := new(Stats)
	out := make(chan Solution, 16)
	type searchSummary struct{ filterPairs, filterAs, boundedPairs uint64 }
	done := make(chan searchSummary, 1)
	go func() {
		fp, fa, bp := searchFiltered(amax, kmax, *workers, bmin, bmax, filter, stats, out)
		done <- searchSummary{fp, fa, bp}
		close(out)
	}()

	found := 0
	for s := range out {
		found++
		c := new(big.Int).Add(s.B, big.NewInt(int64(s.K)))
		fmt.Printf(
			"NEW SOLUTION: %d! * %s! = %s!  (k=%d)\n",
			s.A, s.B.String(), c.String(), s.K,
		)
	}
	ss := <-done

	fmt.Printf("modular survivors      : %d pairs on %d A values\n", ss.filterPairs, ss.filterAs)
	fmt.Printf("inside B-range         : %d pairs\n", ss.boundedPairs)
	fmt.Printf("known nontrivial solution: 6! * 7! = 10!\n")
	fmt.Printf("new solutions found    : %d\n", found)
	fmt.Printf("(A,k) pairs examined   : %d\n", stats.pairs.Load())
	fmt.Printf("reduced kth roots      : %d\n", stats.reducedRoots.Load())
	fmt.Printf("full-root fallbacks    : %d\n", stats.fallbacks.Load())
	fmt.Printf("exact product compares : %d\n", stats.comparisons.Load())
	fmt.Printf("elapsed                : %v\n", time.Since(start))
}