Extending the Verified Range for \(A!B!=C!\) to \(B \leq 10^{1,000,000,000}\)

The factorial Diophantine equation

\[ A!B!=C! \]

asks when a product of two factorials is again a factorial. The numbers involved get very large, but they are not opaque: for every prime \(p\), Legendre's formula gives the exact exponent of \(p\) in \(n!\). For a fixed triple, equality therefore requires matching all prime exponents up to \(C\). This finite but highly coupled system of exact constraints, rather than the size of the factorials alone, drives every advance on the equation.

Because multiplication is commutative, \(A!B!=B!A!\): the triples \((A,B,C)\) and \((B,A,C)\) solve the equation together or not at all. Solutions therefore come in mirror pairs, and it suffices to examine one representative of each pair, so we fix the convention \(A\le B\) in this article. One infinite family of solutions is immediate. If \(C=B+1\), the identity \(C!=(B+1)B!\) turns the equation into \(A!=B+1\), so under the convention \(A\le B\), every \(A\ge3\) yields

\[ B=A!-1, \qquad C=A!, \]

that is,

\[ A!(A!-1)!=(A!)!. \]

The first instances are \(3!5!=6!\), \(4!23!=24!\), and \(5!119!=120!\). Solutions of this shape arise from two consecutive factorials and are called trivial. The equation becomes interesting only when at least two factors separate \(B!\) from \(C!\); accordingly, a solution is called nontrivial when \(A\le B\le C-2\). Exactly one nontrivial solution is known:

\[ 6!\,7!=10!. \]

Cancelling \(7!\) from both sides reduces the identity to \(6!=8\cdot9\cdot10=720\), so it can be checked without expanding either large factorial. The conjecture, attributed in this form to Surányi, is that no other nontrivial solution exists.

Known results progressively confine any hypothetical counterexample. Caldwell verified the conjecture up to \(C\le10^6\). Independently of that computational range, Hajdu, Papp, and Szakács proved that every other nontrivial solution satisfies \(B-A\ge10^6\); since \(A\ge1\), this gives the starting point \(B>10^6\) used in the global estimates below. Habsieger then combined dyadic valuations, explicit Stirling estimates, and a finite computation to show that every further nontrivial solution must satisfy \(B\ge10^{3000}\), and he closed his paper by observing that the computational component of his method was rather basic, so that a smarter algorithm should lead to a much larger bound. The present work takes up that remark. An exact search with a different architecture first reproduces the computation needed at the \(10^{3000}\) boundary. A finite-field method then removes the direct search over \(B\) - exact congruence conditions involving only the two small parameters of the problem - and carries the verified range to \(B\le10^{1{,}000{,}000{,}000}\): any further nontrivial solution must exceed a one followed by one billion zeros.

The question belongs to a broader line of work on products of factorials. Erdős isolated structural consequences of factorial equations; Luca, and later Bhat and Ramachandra, studied when factorials can themselves be products of factorials. Those results explain why the equation is arithmetically rigid, while the explicit ranges needed here come from the later bounds of Hajdu, Papp, Szakács, and Habsieger.

What Must Be Proved

A finite computation can support a theorem only when every hypothetical solution is represented in the search space and every rejection is a logical necessity. The argument therefore has four layers:

1. Rewrite the equation as a product of \(k=C-B\) consecutive integers.
2. Use published analytic inequalities to place every possible \((A,k)\) in a finite rectangle.
3. Apply exact decimal magnitude bounds, then exact congruence conditions modulo several primes, never discarding a pair on approximate evidence.
4. Show that no pair survives all necessary conditions.

For \(10^{3000}\le B\le10^{1{,}000{,}000{,}000}\), conservative forms of Habsieger's bounds leave \(A\le3{,}321{,}928{,}159\) and \(2\le C-B\le33\). Exact decimal magnitude bounds cut this rectangle to \(64{,}660{,}169{,}462\) candidate pairs, and fifty prime-field filters remove every one of them. Subject to Habsieger's published analytic bounds and to the exact computation described below, the resulting computer-assisted statement is

\[ \boxed{ (A,B,C)\ne(6,7,10) \quad\Longrightarrow\quad B>10^{1{,}000{,}000{,}000} } \]

for every nontrivial solution. This does not prove that \((6,7,10)\) is the only nontrivial solution. It extends the explicit verified lower bound on the location of a hypothetical counterexample.

From Three Factorials to Consecutive Integers

Put \(k=C-B\). Dividing \(A!B!=(B+k)!\) by \(B!\) leaves the fundamental equation

\[ \boxed{A!=\prod_{i=1}^{k}(B+i).} \tag{1} \]

Everything that follows works with this form. For fixed \((A,k)\), the product on the right is strictly increasing in \(B\), so at most one integer \(B\) can satisfy (1). Moreover,

\[ B^k<\prod_{i=1}^{k}(B+i)<(B+k)^k, \]

and taking positive \(k\)-th roots,

\[ \boxed{B<(A!)^{1/k}<B+k.} \tag{2} \]

So \(B\) is confined to a window of length \(k\) immediately below \((A!)^{1/k}\), and an unconditional exact fallback needs to test only the \(k+1\) integers from \(\lfloor(A!)^{1/k}\rfloor\) down to \(\lfloor(A!)^{1/k}\rfloor-k\). This removes one dimension of the naive search. It does not remove the real obstacle: constructing \(A!\) is still prohibitive in the terminal range, so a different idea must eliminate \(B\) itself.

Why the Remaining Parameters Are Small

Write \(v_p(N)\) for the exponent of the prime \(p\) in the factorization of \(N\). Legendre's formula evaluates it for factorials:

\[ v_p(n!)=\sum_{j\ge1}\left\lfloor\frac{n}{p^j}\right\rfloor =\frac{n-s_p(n)}{p-1}, \tag{3} \]

where \(s_p(n)\) is the sum of the base-\(p\) digits of \(n\). The sum counts the multiples of \(p\) among \(1,\ldots,n\), then one additional factor from every multiple of \(p^2\), another from every multiple of \(p^3\), and so on. For the digit-sum form, expand \(n=a_0+a_1p+\cdots+a_rp^r\); then

\[ \sum_{j\ge1}\left\lfloor\frac{n}{p^j}\right\rfloor =\sum_{i=1}^{r}a_i(1+p+\cdots+p^{i-1}) =\frac{n-(a_0+\cdots+a_r)}{p-1}, \]

which proves both forms of (3). For \(p=2\) the denominator is one and the digit sum is simply the number of ones in the binary representation, so \(v_2(n!)=n-s_2(n)\). Applying this to \(A!B!=C!\) gives

\[ C=A+B-s_2(A)-s_2(B)+s_2(C). \tag{4} \]

Since \(s_2(C)\ge1\) and \(s_2(n)\le\log_2(n+1)\),

\[ C\ge A+B+1-\frac{\log(A+1)+\log(B+1)}{\log2}. \tag{5} \]

Binary digit sums therefore keep \(C\) within a few bits of \(A+B\). Substituting \(C=B+k\), inequality (5) reads \(A\le k-1+\log_2\bigl((A+1)(B+1)\bigr)\). This dyadic constraint relates \(A\), \(k\), and \(\log B\), but it does not by itself bound \(A\) logarithmically because \(k\) is not yet controlled. The required control comes from combining it with the rapid convex growth of \(\log(C!)=\log(A!)+\log(B!)\). Habsieger made this combination quantitative through an error-controlled form of Stirling's approximation. For \(x>0\), the estimate used in his argument may be written

\[ 0< \log\Gamma(x)-x(\log x-1)-\frac12\log\frac{2\pi}{x} \le\frac{1}{12x}. \]

The central point is not merely that \(n!\) is approximately \(\sqrt{2\pi n}(n/e)^n\), but that the displayed remainder has a known sign and an explicit size, so it can participate in an inequality proof without introducing an uncontrolled approximation. Combining that estimate with (5), Habsieger's Theorem 4 proves, for every hypothetical additional nontrivial solution beyond his verified range,

\[ A\le \frac{\log(B+1)}{\log2} +\frac{2\log\log(B+1)}{\log2} -1.3479, \tag{6} \]

\[ k=C-B\le \frac{\log\log(B+1)}{\log2} -0.8803. \tag{7} \]

These bounds are the analytic front end of the new computation. No stronger asymptotic result is claimed here.

Verification of the Published Record

Reproducing the \(10^{3000}\) result must not assume that result. The lower endpoint follows from the theorem of Hajdu, Papp, and Szakács: \(B-A\ge10^6\) for every nontrivial solution other than \((6,7,10)\), hence \(B>10^6\). Habsieger's Theorems 1 and 2 then provide the globally valid estimates

\[ A\le \frac{\log(B+1)}{\log2} +\frac{2\log\log(B+1)}{\log2} +2.1221, \tag{8} \]

\[ k\le \frac{\log\log(B+1)}{\log2} +1.819. \tag{9} \]

Both right-hand sides increase with \(B\), so it suffices to evaluate them at the conservative upper value \(\log(B+1)<3000\log10+\log2\). To certify the numerical step, use \(\log x=2\sum_{j=0}^{N}z^{2j+1}/(2j+1)+R_N\), where \(z=(x-1)/(x+1)\) and \(0<R_N<2z^{2N+3}/((2N+3)(1-z^2))\). Use 121 terms, \(z=1/3\) for \(\log2\), and \(\log10=3\log2+\log(5/4)\), where the latter uses \(z=1/9\). Before evaluating \(\log L_D\), divide \(L_D\) by \(2^{12}\) for \(D=3000\), or by \(2^{31}\) for \(D=10^9\), so that the series argument lies in \([1,2)\), then add the corresponding multiple of the interval for \(\log2\). Rational interval operations with the displayed tail bound give

\[ 9994.4146766933823 < \frac{L_{3000}}{\log2} +\frac{2\log L_{3000}}{\log2}+2.1221 <9994.4146766933824, \qquad 14.5731460156476 < \frac{\log L_{3000}}{\log2}+1.819 <14.5731460156477, \quad L_{3000}=3000\log10+\log2, \]

and taking integer floors yields the search rectangle

\[ A\le9994, \qquad 2\le k\le14 \qquad (10^6\le B\le10^{3000}). \tag{10} \]

For each \(k\), monotonicity of the product in (1) bounds the possible factorial between exact endpoint products:

\[ \prod_{i=1}^{k}(10^6+i) \le A!\le \prod_{i=1}^{k}(10^{3000}+i). \tag{11} \]

Factorials are generated incrementally only within this interval. Every surviving pair is then localized through an exact integer \(k\)-th root. Floating point supplies nothing but a starting value: Newton's iteration for \(x^k=N\), with an integer floor taken at every step, becomes

\[ x_{n+1}=\left\lfloor \frac{(k-1)x_n+\left\lfloor N/x_n^{k-1}\right\rfloor}{k} \right\rfloor. \]

An approximate logarithm seeds the iteration but is never trusted. After the iteration settles, the code increments or decrements the result until the two integer power comparisons

\[ x^k\le A!<(x+1)^k \tag{12} \]

hold; these comparisons uniquely characterize \(\lfloor N^{1/k}\rfloor\), so the returned root is exact regardless of the quality of the seed.

A dyadic congruence for the candidate value of B

The root interval (2) leaves at most \(k+1\) candidates. The dyadic valuation of \(A!\) usually does much better. Let \(E=v_2(A!)=A-s_2(A)\), and suppose the factor \(B+i\) attains the maximum dyadic valuation \(m=v_2(B+i)\) among the \(k\) consecutive factors. If \(2^m>k-1\), no second factor can be divisible by \(2^m\), because two multiples of \(2^m\) differ by at least \(2^m\); the maximum is then unique. For \(j\ne i\), the unequal-valuation identity gives \(v_2(B+j)=v_2(j-i)\), and summing all valuations in (1) yields

\[ m=E-v_2((i-1)!)-v_2((k-i)!). \tag{13} \]

The exact valuation of the distinguished factor forces

\[ \boxed{B\equiv2^m-i\pmod{2^{m+1}}.} \tag{14} \]

An integer with possibly thousands of digits is thereby pinned to a single residue class modulo a large power of two. Because the root interval in (2) has length \(k<2^{m+1}\), it contains at most one representative of this residue class for each possible position \(i\). The implementation nevertheless constructs representatives from two adjacent residue blocks: its reduced root determines only the quotient \(\lfloor(A!)^{1/k}/2^{m+1}\rfloor\), so the valid representative can lie in that block or the immediately preceding one. The final exact product comparison selects at most one of them. A further exact identity makes the reduced root available without constructing all low bits:

\[ \left\lfloor\frac{N^{1/k}}{2^s}\right\rfloor =\left\lfloor \left(\left\lfloor\frac{N}{2^{ks}}\right\rfloor\right)^{1/k} \right\rfloor. \tag{15} \]

To prove (15), put \(t=\lfloor N^{1/k}/2^s\rfloor\). Its definition says

\[ t2^s\le N^{1/k}<(t+1)2^s. \]

All quantities are nonnegative, so raising to the \(k\)-th power preserves the inequalities. Dividing by \(2^{ks}\) gives

\[ t^k\le\frac{N}{2^{ks}}<(t+1)^k, \]

and since \(t^k\) is an integer, taking the floor of the middle term changes neither bound:

\[ t^k\le\left\lfloor\frac{N}{2^{ks}}\right\rfloor<(t+1)^k. \]

The exact floor \(k\)-th root is therefore \(t\), proving (15). Instead of extracting a root of all of \(A!\), the implementation right-shifts by \(ks\) bits first and works on the smaller quotient, preserving the high part of the answer exactly.

If the unique-maximum hypothesis is unavailable, the implementation falls back to (2); the shortcut never decides correctness. Invoked as ./factorial_solver -digits 3000 -workers 10, the exact lower-range solver examined \(87{,}408\) admissible \((A,k)\) pairs, performed \(47{,}248\) reduced exact roots and \(310{,}161\) exact product comparisons, and found no new solution. The reported ten-worker run took \(8.58108925\) seconds; a fresh run under Go 1.26 while preparing this page reproduced all four workload counters exactly in \(8.9863595\) seconds. The empty result reproduces the computational conclusion needed at the published boundary through a search architecture different from Habsieger's reported Maple computation, which covered the earlier rectangle in about 28 hours.

Bounds for One Billion Decimal Exponents

With the published boundary reproduced, Habsieger's stronger inequalities (6) and (7) apply on \(10^{3000}\le B\le10^D\). Since \(B+1<2\cdot10^D\), set \(L_D=D\log10+\log2\). The same outward-rounded logarithm series used above gives, for \(D=10^9\),

\[ 3{,}321{,}928{,}156.7406770021 <\frac{L_D}{\log2} +\frac{2\log L_D}{\log2}-1.3479 <3{,}321{,}928{,}156.7406770022, \]

\[ 30.22030732712 <\frac{\log L_D}{\log2}-0.8803 <30.22030732713. \]

The tight integer consequences are \(A\le3{,}321{,}928{,}156\) and \(k\le30\). The implementation deliberately pads both limits and searches

\[ \boxed{A\le3{,}321{,}928{,}159,\qquad2\le k\le33.} \tag{16} \]

The displayed interval bounds, not the solver's floating-point evaluation, prove that every admissible integer pair lies inside this rectangle. The runtime calculation reproduces the padded values in (16); its two extra units are operational slack, while the certified intervals supply the formal inclusion.

Even so reduced, the problem cannot be finished by enumeration. The interval \(10^{3000}\le B\le10^{1{,}000{,}000{,}000}\) contains on the order of \(10^{1{,}000{,}000{,}000}\) integers, and no constant-factor speedup, vectorization, or distribution changes that. The variable \(B\) itself has to be eliminated. Both sieves that follow do exactly this: each tests a necessary condition that depends on \((A,k)\) alone.

An Exact Decimal Magnitude Sieve

Define integer exponent bounds

\[ L(A)=\sum_{n=1}^{A}\lfloor\log_{10}n\rfloor, \qquad U(A)=\sum_{n=1}^{A}\lceil\log_{10}n\rceil. \]

Their summands change only at powers of ten, so both sums are computed without floating-point logarithms. They give the exact enclosure

\[ 10^{L(A)}\le A!\le10^{U(A)}. \tag{17} \]

If \(B\ge10^d\), equation (1) implies \(A!>B^k\ge10^{dk}\), hence every candidate must satisfy

\[ \boxed{U(A)>dk.} \tag{18} \]

If \(B\le10^D\), then \(B+i<2\cdot10^D\) throughout the relevant range. Since \(2^k<10^k\), equation (1) implies \(A!<10^{(D+1)k}\), and therefore

\[ \boxed{L(A)<(D+1)k.} \tag{19} \]

Conditions (18) and (19) are necessary integer inequalities: they contain no approximation and cannot reject a genuine solution. Because \(L\) and \(U\) are nondecreasing in \(A\), the surviving \(A\) values of each \(k\) form one contiguous interval, so a single sweep suffices to record every interval endpoint. For \(d=3000\) and \(D=10^9\), the two conditions leave exactly \(64{,}660{,}169{,}462\) pairs \((A,k)\).

The Prime-Field Image Sieve

The main reduction comes from reading equation (1) modulo a prime. A congruence modulo \(q\) retains only the remainder upon division by \(q\): the notation \(x\equiv y\pmod q\) means that \(q\) divides \(x-y\). When \(q\) is prime, the remainders \(0,1,\ldots,q-1\), with addition and multiplication modulo \(q\), form the finite field \(\mathbb F_q\), in which every nonzero element has a multiplicative inverse. These elementary facts are all the modular part of the computation needs.

Let \(q>A_{\max}\) be prime and define the rising-factorial polynomial over \(\mathbb F_q\),

\[ f_k(X)=\prod_{i=1}^{k}(X+i). \tag{20} \]

A polynomial maps each field element \(x\) to one field element \(f_k(x)\). Its image \(f_k(\mathbb F_q)\) is the set of remainders that occur as \(x\) runs through \(0,1,\ldots,q-1\); it need not contain every remainder, and membership in it is precisely the information the sieve uses. If (1) holds and \(x=B\bmod q\), reduction modulo \(q\) gives \(A!\equiv f_k(x)\pmod q\), so every genuine solution satisfies the exact necessary condition

\[ \boxed{A!\bmod q\in f_k(\mathbb F_q).} \tag{21} \]

Proof of the finite-field condition

Assume that integers \(A,B,k\) satisfy (1), and let \(x\) be the remainder of \(B\) upon division by \(q\). Then \(B\equiv x\pmod q\), so adding the same integer \(i\) to both sides gives \(B+i\equiv x+i\pmod q\). Congruences may be multiplied, hence

\[ A!=\prod_{i=1}^{k}(B+i) \equiv\prod_{i=1}^{k}(x+i) =f_k(x) \pmod q. \]

The value \(f_k(x)\) belongs to the image by definition, proving (21). Its contrapositive is the filter actually implemented: if \(A!\bmod q\) is absent from the image, no integer \(B\) whatsoever can rescue the pair \((A,k)\). One modular failure excludes every integer in every residue class simultaneously, including integers far too large to construct.

The predicate depends on \((A,k)\), not on the size of \(B\); this is the asymmetry that lets the method reach extreme ranges. Choosing \(q>A_{\max}\) ensures \(A!\not\equiv0\pmod q\), avoiding the weak condition produced when \(q\le A\), where the residue \(0\) lies in almost every image. For \(D=10^9\), the first filter prime is \(3{,}321{,}928{,}171\).

Dense image construction

For fixed \(x\), all active values of \(k\) are generated in one pass through the recurrence \(p_0=1\) and \(p_k=p_{k-1}(x+k)\bmod q\), so the value for \(k+1\) costs one multiplication beyond the value for \(k\). A bitset records every value attained by each \(f_k\): bit number \(y\) is set exactly when some \(x\) has produced \(f_k(x)=y\), so a later membership query costs one array lookup and one mask. All workers share one bitset per \(k\) and publish bits atomically - a bit is only ever set, never cleared, during construction, so concurrent workers cannot lose a value. The factorial residues are then swept with \(r_A=r_{A-1}A\bmod q\), parallelized through a segmented prefix product over 64-bit-aligned \(A\) ranges, and \((A,k)\) survives precisely when the bit for \(r_A\) is set.

for x := lo; x < hi; x++ { // each worker owns one range of 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)
        }
    }
}

One image requires \(q\) bits rather than \(q\) machine words. At \(q\approx3.32\cdot10^9\), that is about \(415\) MB per active image, so the implementation processes the active \(k\) values in groups whose images fit into a configurable memory budget. Grouping changes scheduling only; the intersection of necessary conditions is the same in any order.

How Fast a Prime Filter Eliminates

How many filters will be needed is a question about the size of the images \(f_k(\mathbb F_q)\), and the answer has a clean structure. Substituting \(x\mapsto-(k+1)-x\) sends each factor \(x+i\) to \(-(x+k+1-i)\), so the factors are permuted up to sign and

\[ f_k(-(k+1)-x)=(-1)^k f_k(x). \]

For odd \(k\) this only says the image is symmetric under negation, and the image behaves like that of a generic polynomial map. If every value's preimage count behaved like a Poisson variable with mean one - the standard heuristic for a map of a large finite set to itself - a fraction \(1-e^{-1}\approx0.632\) of the field would be covered. For even \(k\), however, the symmetry has a fixed shape: completing the square with \(u=x+\tfrac{k+1}{2}\) (the division by two is exact in \(\mathbb F_q\) for odd \(q\)) turns the product into

\[ f_k(x)=\prod_{j=1}^{k/2}\Bigl(u^2-\bigl(j-\tfrac12\bigr)^2\Bigr), \]

a polynomial in \(u^2\). Its image is therefore the image of a degree-\(k/2\) polynomial evaluated only on the \((q+1)/2\) squares of the field, and applying the same Poisson model to half a domain suggests density \(1-e^{-1/2}\approx0.393\). The smallest cases are exact: \(k=2\) gives a quadratic with image size \((q+1)/2\). For \(k=3\), the map is a translate of \(u\mapsto u^3-u\); its image density approaches \(2/3\) with finite-field corrections rather than equalling \(2/3\) for every prime (for example, it is \(3/5\) over \(\mathbb F_5\)). Direct enumeration at \(q=100{,}003\) illustrates the predicted pattern:

k       2      3      4      5      6      7      8      9     10
|image|/q   0.5000 0.6667 0.3750 0.6337 0.3950 0.6324 0.3936 0.6293 0.3930

Under this model, a filter applied to a mixture of parities retains roughly half of all candidate pairs, the even-\(k\) classes tend to disappear first at rate \(\approx0.393\) per prime, and later filters approach the odd-\(k\) rate \(\approx0.632\). The observed \(D=10^8\) survivor ratios are consistent with that picture:

0.5167  0.5444  0.5688  0.5884  0.6031
0.6134  0.6203  0.6249  0.6278  0.6298

they climb from the mixed rate toward \(1-e^{-1}=0.6321\ldots\) as the even classes vanish. This does not establish a distribution law for factorial residues or independence between filters. If those additional assumptions are used as a runtime model, an odd-\(k\) population of size \(N\) is expected to require about \(\ln N/(-\ln(1-e^{-1}))\approx\ln N/0.459\) filters, followed by fluctuations in the small tail. Every actual rejection remains an exact congruence failure, so the heuristic affects performance expectations only.

Sparse Filtering by a Polynomial GCD

Dense image construction is wasteful once few candidates remain: it answers \(q\) membership questions when a handful are needed. For a single candidate, set \(y=A!\bmod q\) and \(F(X)=f_k(X)-y\); the question is whether \(F\) has a root in \(\mathbb F_q\). Fermat's theorem gives \(a^q=a\) for every field element, and a monic degree-\(q\) polynomial cannot have more than \(q\) roots, so

\[ X^q-X=\prod_{a\in\mathbb F_q}(X-a), \]

a factorization containing every linear polynomial over \(\mathbb F_q\) exactly once. The greatest common divisor of two polynomials is their highest-degree common factor, normalized to leading coefficient one; a linear factor \(X-a\) divides \(F(X)\) exactly when \(F(a)=0\). Thus \(F\) has a root in \(\mathbb F_q\) exactly when it shares a nonconstant factor with \(X^q-X\):

\[ \boxed{ y\in f_k(\mathbb F_q) \iff \deg\gcd\bigl(f_k(X)-y,X^q-X\bigr)>0. } \tag{22} \]

The degree-\(q\) polynomial is never constructed. Binary exponentiation computes \(R(X)=X^q\bmod F(X)\), after which the test is \(\gcd(F,R-X)\). All intermediate polynomials have degree below \(k\le33\), giving roughly \(O(k^2\log q)\) field operations with the elementary arithmetic used here, in place of \(O(qk)\) work for a dense image. At \(q\approx3.3\cdot10^9\) one query costs on the order of \(50\) microseconds, so the sparse test overtakes a dense pass as soon as fewer than roughly \(10^8\) pairs remain.

To justify the reduction, polynomial division writes \(X^q-X=Q(X)F(X)+H(X)\), where \(H=(X^q-X)\bmod F\). A polynomial divides both \(F\) and \(X^q-X\) exactly when it divides both \(F\) and \(H\), because subtracting \(QF\) does not change common divisors. Thus \(\gcd(F,X^q-X)=\gcd(F,H)\). Combining this with the linear-factor argument proves (22), rather than merely motivating it.

func hasRisingFactorialRoot(base []uint64, y, q uint64) bool {
    f := append([]uint64(nil), base...)
    f[0] = subMod(f[0], y, q)
    f = polyTrim(f)

    xq := polyPowMod([]uint64{0, 1}, q, f, q)
    h := polySub(xq, []uint64{0, 1}, q)
    return polyGCDDegree(f, h, q) > 0
}

Dense bitset membership and sparse GCD membership implement the same predicate (21). Changing the crossover threshold changes performance only.

Factorial Residues in Sparse Mode

The sparse test still needs \(A!\bmod q\). The implemented segmented forward path first computes every segment product and then traverses the segments containing requested values again; its cost is linear in the traversed range, with two modular-product passes over the relevant segments. No isolated benchmark is used here to assign it a wall-clock time. When the requested \(A\) values lie close to \(q\), Wilson's theorem supplies a shorter reverse route. Define

\[ T_A=\prod_{n=A+1}^{q-1}n. \]

Since \((q-1)!\equiv-1\pmod q\),

\[ \boxed{A!\equiv-T_A^{-1}\pmod q.} \tag{23} \]

Wilson's theorem itself follows by pairing every nonzero residue with its multiplicative inverse: all pairs contribute \(1\), and the only residues equal to their own inverses solve \(x^2=1\), namely \(1\) and \(-1\), so the full product is \(-1\). Splitting \((q-1)!=A!\,T_A\) and multiplying by \(T_A^{-1}\) gives (23). The implementation compares a forward cost of \(A\) multiplications with a reverse cost of \(q-1-A\), uses the cheaper sweep, and batch-inverts residues when several values of \(A\) are needed.

Batch inversion replaces many expensive inversions by one. For nonzero residues \(z_1,\ldots,z_n\), form prefix products \(P_i=z_1\cdots z_i\), invert only \(P_n\), and walk backward: at step \(i\), multiplying the current inverse of \(P_i\) by \(P_{i-1}\) gives \(z_i^{-1}\), and multiplying it by \(z_i\) gives \(P_{i-1}^{-1}\). Here \(q>A_{\max}\ge A\), so every factor from \(A+1\) through \(q-1\), every tail \(T_A\), and every prefix product is nonzero in \(\mathbb F_q\); all displayed inverses therefore exist. The procedure uses one field inversion and \(O(n)\) multiplications, meaning that the operation count grows at most proportionally to \(n\).

Dense and Sparse Work Are the Same Test

Dense mode spends \(q\) evaluations once and then answers many membership queries cheaply. Sparse mode spends roughly \(O(k^2\log q)\) operations per queried pair but allocates no \(q\)-bit image. The implementation switches to the sparse test when the survivor count falls below a configurable threshold.

This crossover is not part of the proof. Dense enumeration asks whether some \(x\) satisfies \(f_k(x)=y\) by trying every \(x\); sparse mode asks the identical question through (22). A wrong performance estimate can make the run slower, but cannot alter the mathematically defined survivor set when both implementations are correct.

Correctness and Validation

Let \(\mathcal C_0\) be the candidate set after (16), (18), and (19). For each selected prime \(q\), define

\[ \mathcal C_q= \{(A,k)\in\mathcal C_0:A!\bmod q\in f_k(\mathbb F_q)\}. \]

A real solution belongs to every \(\mathcal C_q\). The primes need not behave independently, because no probability estimate is involved. For any finite set of primes \(q_i>A_{\max}\), if

\[ \mathcal C_0\cap\mathcal C_{q_1}\cap\cdots\cap\mathcal C_{q_r}=\varnothing, \tag{24} \]

then no candidate pair can come from an integer solution in the searched range. Conversely, modular survival would not prove a solution; it would only trigger exact reconstruction of \(B\) and comparison of the consecutive product.

Theorem. Let \(D\ge3000\), and let conservative integers \(A_{\max}\) and \(k_{\max}\) contain all pairs permitted by Habsieger's inequalities for \(10^{3000}\le B\le10^D\). Retain every pair satisfying the magnitude conditions (18) and (19). Let \(Q\) be any finite set of primes greater than \(A_{\max}\). If every retained pair fails condition (21) for at least one \(q\in Q\), then no nontrivial solution other than \((6,7,10)\) exists with \(10^{3000}\le B\le10^D\).

Proof. A hypothetical solution lies inside the analytic rectangle and satisfies both magnitude conditions by (1). Reducing the same integer equality modulo every \(q\in Q\) forces (21), so its pair would survive every filter. An empty intersection contradicts its existence. \(\square\)

The implementation preserves this argument with exact guards: the analytic limits are enclosed by the rational intervals above; the decimal sieve uses integer sums; integer powers certify every Newton root; both finite-field modes compute (21) exactly; atomic image publication and disjoint bitset words prevent lost concurrent updates; and grouping changes only work order. No density model, random-prime assumption, floating-point equality, or timeout contributes to the empty intersection.

The self-test compares the interval decimal sieve with direct per-\(A\) evaluation, segmented factorial residues with a direct sweep, the reverse Wilson path with direct residues, and dense image membership with the sparse GCD test for every field value over primes from \(29\) through \(101\) and \(2\le k\le8\). These checks exercise the posted implementation, but they are not a separately implemented exclusion verifier. The downloadable source and stdout therefore form a reproducible computation record rather than an independently checkable certificate.

The Verified Runs

Two generations of the finite-field solver exist. Revision 1 established the method and originally covered \(D=10^8\). Revision 2, the posted source, evaluates the identical predicates with a re-engineered core: shared atomically-filled images instead of per-worker copies, a parallel survivor sweep over aligned segments, and interval-based construction of the decimal sieve. At \(D=10^6\), the archived logs show identical survivor counts after every one of the 39 primes and the same final empty result. Equal counts do not prove that the intermediate candidate sets were identical, and the revisions share predicates and design, so this is a useful regression check rather than independent verification.

All runs listed on this page were executed on the same ten-core test system (Apple M1 Max, 64 GB, Go 1.26). The full revision-1 source is not part of the archive; its recorded digest identifies it but does not make the cross-check reproducible from the current downloads.

The record computation

For \(D=10^9\), the run used fifty primes greater than \(A_{\max}=3{,}321{,}928{,}159\). The complete validated sequence was

3321928171  3321928307  3321928417  3321928531  3321928661  3321928777
3321928919  3321929029  3321929129  3321929257  3321929357  3321929461
3321929563  3321929683  3321929789  3321929909  3321930011  3321930173
3321930289  3321930397  3321930517  3321930671  3321930793  3321930893
3321931007  3321931111  3321931213  3321931313  3321931429  3321931531
3321931637  3321931763  3321931897  3321932063  3321932179  3321932353
3321932473  3321932591  3321932701  3321932819  3321932947  3321933047
3321933169  3321933271  3321933379  3321933541  3321933661  3321933781
3321933893  3321934003

Every entry was rechecked outside the filtering loop with the deterministic Miller-Rabin test for 64-bit integers, and each equals the value demanded by the prime rule, so the sequence can be regenerated prime for prime by anyone repeating the run. Sixteen filters ran in dense image mode and carried the candidate set from \(64.7\) billion pairs to the sparse crossover below thirty million; thirty-four polynomial-GCD filters completed the elimination.

size-bound candidates   : 64660169462 pairs
filters actually used   : 50
filter build            : 3h38m8.345710959s

modular survivors       : 0 pairs
inside B-range          : 0 pairs
known nontrivial solution: 6! * 7! = 10!
new solutions found     : 0
(A,k) pairs examined    : 0
reduced kth roots       : 0
full-root fallbacks     : 0
exact product compares  : 0
elapsed                 : 3h38m8.345874s

The terminal zero counters matter: after the finite-field stage, no candidate required construction of a large factorial, an exact root, or a final product comparison. The result was determined entirely by necessary conditions. The observed survivor trajectory is consistent with the heuristic model above: it begins at ratio \(0.516\) on the first filter and approaches \(1-e^{-1}\) after the even classes disappear. This agreement is descriptive and is not used in the exclusion argument.

Scaling

Nothing number-theoretic happens at \(10^{1{,}000{,}000{,}000}\); the limit is computational. If \(B\le10^D\), Habsieger's bound gives

\[ A_{\max}\sim\log_2(10)D\approx3.321928094D, \qquad k_{\max}=O(\log D), \]

so the computation is governed by the decimal exponent \(D\), not by the number \(10^D\). Dense image construction costs roughly \(O(rqk_{\mathrm{active}})\) modular multiplications for \(r\) filters with \(q\approx A_{\max}\); only under the distribution assumptions above does the model suggest \(r\approx\ln N/0.459\) for \(N\) odd-\(k\) candidates. The sparse stage removes the \(O(q)\) image cost for the tail. The single-word multiplication path requires \(q<2^{32}\), which holds up to \(D\approx1.29\cdot10^9\). Above that threshold the current code uses exact 128-by-64-bit remainder reduction, not Montgomery multiplication.

Memory has three distinct components. The persistent candidate sets occupy \(8\sum_k\lceil(A_{\mathrm{hi}}(k)+1)/64\rceil\) bytes, while one dense image occupies \(8\lceil q/64\rceil\) bytes. At \(D=10^9\), these are approximately \(7.53\) GiB and \(0.387\) GiB per image. Extrapolating the same exact decimal intervals to \(D=10^{10}\) gives approximately \(84.3\) GiB of candidate bitsets and \(3.87\) GiB per image. The configured image-memory budget limits only the number of simultaneous images, so grouping alone cannot fit that projected run into the present 64-GB system. Dense mode additionally uses about \(8\cdot2^{16}\) bytes per worker for residue buffers plus small segment arrays; sparse mode allocates pair, residue, and rejection structures proportional to the current survivor count. A further order of magnitude therefore requires a different candidate representation or external/distributed storage in addition to arithmetic and parallelism improvements.

Relation to the Earlier Method

Habsieger's analysis supplies the indispensable bounds on \(A\) and \(k\); the improvement here is computational, not analytic. His computation also develops a more specialized localization of \(B\), commonly denoted \(B_k(A)\), around \(\lceil(A!)^{1/k}-(k+1)/2\rceil\). The published lemma states a range through \(k=20\), while its following proof text reports a Maple polynomial verification for \(2\le k\le12\). The extension reported here does not depend on that finite-range assertion. It uses only the unconditional interval (2) and the independently derived dyadic congruence (14); in the large-range stage it discards \((A,k)\) pairs through the universal image condition (21). This separation makes the computational implication short enough to audit directly.

Habsieger also noted that his argument used no \(p\)-adic information beyond \(p=2\). The image condition (21) responds to that observation obliquely rather than literally: instead of computing valuations at odd primes, it forces the entire residue of \(A!\) into the image of a degree-\(k\) polynomial for many primes \(q\) at once, which constrains the multiplicative structure of all \(k\) factors simultaneously. The current sieve uses the 2-adic valuation only in the lower-range verification and claims no new odd-prime \(p\)-adic theorem; such restrictions could reduce later computations, but they would require their own proofs and are listed as a possible optimization rather than silently assumed.

Scope of the Result

The computation does not settle Surányi's conjecture. It relies on Habsieger's analytic inequalities and leaves the infinite region \(B>10^{1{,}000{,}000{,}000}\) untouched. It also does not improve the asymptotic estimates for \(A\) or \(C-B\). The appropriate conclusion is narrower:

\[ \boxed{ A!B!=C!,\ A\le B\le C-2,\ B\le10^{1{,}000{,}000{,}000} \quad\Longrightarrow\quad (A,B,C)=(6,7,10). } \]

The decimal exponent increases from \(3000\) to \(10^9\), a factor of \(333{,}333.\overline{3}\). The corresponding numerical ranges differ by the factor \(10^{999{,}997{,}000}\). To the best of our knowledge, no later published computation has superseded Habsieger's \(10^{3000}\) bound for this equation; that literature statement is separate from the mathematical theorem.

Reproducibility

Core implementation reference

The following excerpts preserve the remaining mathematical core of the posted Go implementation. Production code additionally needs error handling, concurrency control, memory budgeting, and regression tests. Dense image construction and the sparse polynomial root test already appear in their respective sections above.

A.1 Conservative theorem bounds

func theoremBounds(digits int) (amax, kmax int) {
  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
}

These runtime values reproduce the padded rectangle (16), but they are not its proof. The rational logarithm series and outward-rounded intervals displayed above certify the unpadded integer limits; the extra two units then provide explicit operational slack.

A.2 Exact decimal candidate intervals

for a := 2; a <= amax; a++ {
  updateExactDecimalExponents(a, &lowerExp, &upperExp)

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

  // L < (maxDigits+1)*k fails exactly for k <= L/(maxDigits+1).
  failK := int(lowerExp / int64(maxDigits+1))
  for kDead < min(failK, kmax) {
    kDead++
    aHi[kDead] = a - 1
  }
}
// The candidate set of k is the interval [aLo[k], aHi[k]],
// filled into its bitset with whole-word writes.

The helper updates the sums in (17) only when \(a\) crosses an exact power of ten. The interval endpoints are the strict integer forms of (18)-(19); monotonicity of \(L\) and \(U\) makes each candidate set one interval, and the self-test compares this construction bit for bit against the direct per-\(A\) evaluation.

A.3 Rising-factorial polynomial

func risingPoly(k int, q uint64) []uint64 {
  p := []uint64{1}
  for i := 1; i <= k; i++ {
    next := make([]uint64, len(p)+1)
    ii := uint64(i) % q
    for j, coefficient := range p {
      next[j] = addMod64(
        next[j], mulModFast(coefficient, ii, q), q,
      )
      next[j+1] = addMod64(next[j+1], coefficient, q)
    }
    p = next
  }
  return p
}

Multiplying the current coefficient array by \(X+i\) at each iteration constructs \((X+1)(X+2)\cdots(X+k)\) over \(\mathbb F_q\) without symbolic algebra software.

A.4 Wilson-based sparse residues

forwardCost := uint64(maxA)
reverseCost := q - 1 - uint64(minA)

if forwardCost <= reverseCost {
  // Segmented parallel prefix product over 1, 2, ..., maxA.
} else {
  // Sweep backward from q-1 and use
  // A! * tail(A) = (q-1)! = -1 mod q.
}

A.5 Exact dyadic valuation and candidates

func v2Factorial(n int) int {
  return n - bits.OnesCount(uint(n))
}

E := v2Factorial(a)
for i := 1; i <= k; i++ {
  m := E - v2Factorial(i-1) - v2Factorial(k-i)
  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)))
  // B is congruent to residue modulo modulus.
}

A.6 Reduced exact root

baseShift := k * (mMin + 1)
scaled := new(big.Int).Rsh(
  new(big.Int).Set(fact), uint(baseShift),
)
qBase := nthRootFloor(scaled, k)

This is the direct implementation of (15). The function nthRootFloor performs the post-correction certified by (12).

A.7 Universal fallback

func fallbackCandidates(fact *big.Int, k int) []*big.Int {
  root := nthRootFloor(fact, k)
  out := make([]*big.Int, 0, k+1)
  for delta := 0; delta <= k; delta++ {
    b := new(big.Int).Sub(
      new(big.Int).Set(root), big.NewInt(int64(delta)),
    )
    if b.Sign() >= 0 {
      out = append(out, b)
    }
  }
  return out
}

The list may contain values outside the strict interval, but (2) proves that it cannot omit a solution.

A.8 Final exact comparison

func productCmp(fact, b *big.Int, k int) int {
  product := big.NewInt(1)
  term := new(big.Int)
  for i := 1; i <= k; i++ {
    term.Set(b)
    term.Add(term, big.NewInt(int64(i)))
    product.Mul(product, term)
    if product.Cmp(fact) > 0 {
      return 1
    }
  }
  return product.Cmp(fact)
}

Early exit is exact because all remaining factors are positive. If this stage is reached, equality of arbitrary- precision integers, not a numerical tolerance, makes the final decision.

The excerpts above isolate the mathematical steps, while the complete downloadable solvers provide the executable verification. Their exact source files are identified by the SHA-256 digests 80a6abf56ec92c1b4c3b8d06ea171840b28a880fb1a9a8ce597dac7e62f71168 for the lower-range solver and 5a110d4547ce8bffde944b6d0fe5e3987f282bb09c1c65ae12645850dba0a557 for the hybrid finite-field solver (revision 2). The recorded runs used byte-identical source bodies beginning at //go:build ignore, with respective SHA-256 digests a8814c1d9b146550871349ea1b6fa555537c194d4ea8d5b75b07d6aaabb4d216 and 3f196533a42bacf36093e955090ac1fb34a2391b98213a0ed183891913b1efcd. The current downloads differ only by six prepended comment lines containing publication-license metadata. The following commands run the hybrid solver's internal validation, build the executable, and start the standalone extension at the published boundary:

go run factorial_solver_hybrid.go -self-test

go build -o factorial_solver_hybrid factorial_solver_hybrid.go

./factorial_solver_hybrid \
    -min-digits 3000 \
    -digits 1000000000 \
    -workers 10 \
    -filter-primes 100 \
    -sparse-threshold 30000000 \
    -filter-memory-mb 16384 \
    -progress

The sparse crossover and dense-memory limits affect running time and memory consumption, but not the mathematical survivor set.

Computational materials

The excerpts on this page explain the algorithms; the linked complete sources are the reviewable research artifacts. The complete stdout logs of every run reported on this page - the lower-range reproduction, both engine revisions at \(D=10^6\) with their full survivor trajectories, the posted source at \(D=10^7\) and \(D=10^8\), and the \(D=10^9\) record run - are archived alongside the sources. The revision-2 and lower-range source bodies are tied to those logs by the hashes above; the original binaries and the revision-1 source are not archived. The logs can be reproduced but cannot by themselves establish that the recorded execution was correct, and no separate independent verifier is currently supplied.

References

• Bhat2010K. G. Bhat and K. Ramachandra, A Remark on Factorials That Are Products of Factorials, Mathematical Notes 88 (2010), 317-320.
• Caldwell1994C. Caldwell, The Diophantine Equation A!B!=C!, Journal of Recreational Mathematics 26 (1994), 128-133.
• Erdos1976P. Erdős, Problems and Results on Number Theoretic Properties of Consecutive Integers and Related Questions, Congressus Numerantium XVI (1976), 25-44.
• Erdos1993P. Erdős, A Consequence of a Factorial Equation, American Mathematical Monthly 100 (1993), 407-408.
• Guy2004R. K. Guy, Unsolved Problems in Number Theory, 3rd ed., Springer, 2004, Section B23.
• Hajdu2018L. Hajdu, Á. Papp, and T. Szakács, On the Equation A!B!=C!, Journal of Number Theory 187 (2018), 160-165.
• Habsieger2019L. Habsieger, Explicit Bounds for the Diophantine Equation A!B!=C!, Fibonacci Quarterly 57(1) (2019), 21-28.
• Luca2007F. Luca, On Factorials Which Are Products of Factorials, Mathematical Proceedings of the Cambridge Philosophical Society 143 (2007), 533-542.