// @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

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 conservative integer bounds obtained from Habsieger's
// explicit inequalities for every hypothetical nontrivial solution other than
// (6,7,10):
//
//	A <= log(B+1)/log(2) + 2 log log(B+1)/log(2) + 2.1221
//	k <= log log(B+1)/log(2) + 1.819,  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) {
	L := float64(digits)*math.Ln10 + math.Ln2

	amax = int(math.Floor(
		L/math.Ln2 +
			2*math.Log(L)/math.Ln2 +
			2.1221,
	))

	kmax = int(math.Floor(
		math.Log(L)/math.Ln2 + 1.819,
	))

	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
}

func searchK(
	k, aGlobal int,
	bmin, bmax *big.Int,
	stats *Stats,
	out chan<- Solution,
) {
	// For this fixed k, A! has to lie between the products at bmin/bmax.
	// These exact bounds eliminate most A values before any root work.
	pmin := big.NewInt(1)
	pmax := big.NewInt(1)
	t := new(big.Int)

	for i := 1; i <= k; i++ {
		t.Set(bmin)
		t.Add(t, big.NewInt(int64(i)))
		pmin.Mul(pmin, t)

		t.Set(bmax)
		t.Add(t, big.NewInt(int64(i)))
		pmax.Mul(pmax, t)
	}

	fact := big.NewInt(1)
	mul := new(big.Int)

	for a := 2; a <= aGlobal; a++ {
		mul.SetInt64(int64(a))
		fact.Mul(fact, mul)

		if fact.Cmp(pmin) < 0 {
			continue
		}
		if fact.Cmp(pmax) > 0 {
			break
		}

		stats.pairs.Add(1)

		cands, ok := v2Candidates(fact, a, k, stats)
		if !ok {
			cands = fallbackCandidates(fact, k, stats)
		}

		if b := findAmongCandidates(
			fact, cands, k, bmin, bmax, a, stats,
		); b != nil {
			out <- Solution{A: a, B: new(big.Int).Set(b), K: k}
		}
	}
}

func main() {
	digits := flag.Int(
		"digits", 3000,
		"prove-search hypothetical nontrivial solutions with B <= 10^digits",
	)
	workers := flag.Int(
		"workers", runtime.NumCPU(),
		"maximum number of k values searched concurrently",
	)
	flag.Parse()

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

	start := time.Now()

	amax, kmax := theoremBounds(*digits)
	bmin := new(big.Int).Exp(big.NewInt(10), big.NewInt(6), nil)
	bmax := new(big.Int).Exp(
		big.NewInt(10),
		big.NewInt(int64(*digits)),
		nil,
	)

	fmt.Printf("searching hypothetical solutions other than (6,7,10)\n")
	fmt.Printf("B range : 10^6 .. 10^%d\n", *digits)
	fmt.Printf("bounds  : A <= %d, 2 <= k=C-B <= %d\n", amax, kmax)
	fmt.Printf("workers : %d\n\n", *workers)

	stats := new(Stats)
	out := make(chan Solution, 16)
	sem := make(chan struct{}, *workers)

	var wg sync.WaitGroup
	for k := 2; k <= kmax; k++ {
		wg.Add(1)
		sem <- struct{}{}

		go func(k int) {
			defer wg.Done()
			defer func() { <-sem }()
			searchK(k, amax, bmin, bmax, stats, out)
		}(k)
	}

	go func() {
		wg.Wait()
		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,
		)
	}

	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))
}