Oarecum din întâmplare — scriind din greșală "-" în loc de "+", undeva în textul funcției estimate_H1() — am observat la un moment dat că putem trata „unitar” (și nu doar ca formulare, ceea ce chiar observasem) cele două chestiuni puse în [1]:
Cum găsim cel mai mic și cel mai mare număr, al cărui factorial este cu $n$ cifre mai scurt decât factorialul succesorului său?

În [1] separasem cazurile și găsisem câte o formulă estimativă pentru „cel mai mic” și respectiv, pentru „cel mai mare” (pe care le notasem pentru un $n$ dat, cu $H_1$ și $H_2$); cele două programe $\mathbf{R}$ pe care le asociam acestor formule nimereau termenii respectivi, dar nu neapărat exact (și-i puteam preciza doar pentru $n \le 10$).

Bineînțeles că ne-am gândit să mărim „precizia” de calcul (permițând și $n > 10$) și în plus, să instrumentăm cumva acreditarea și (când ar fi cazul) corectarea estimărilor produse prin cele două formule. După mai multe încercări (ba în $\mathbf{R}$, ba în $\mathbf{C}$ sau în Python) de a adopta una sau alta dintre bibliotecile de calcul în „precizie arbitrară” — am găsit drept cel mai convenabil, să folosim Python împreună cu gmpy2.

Dar mai întâi să redăm prima încercare, în care doar am translatat codul formulei pentru $H_1$, din $\mathbf{R}$ (unde ne lovisem de limitarea $2^{31}-1$ a întregilor) în Python:

from math import log, log10, sqrt, lgamma, floor, fmod
def estimate_H1 (n):  # estimează H1
    S1 = 10**(n-1)
    delta = fmod(lgamma(S1 + 1) / log(10), 1.0)
    D = sqrt(1 + 8*S1*(1-delta)*log(10)) / 2  # v. formula în [1]
    return floor(S1 - D - 0.5)  # corect era "S1 + D"... (plus, nu minus!)
H1 = [estimate_H1(n+1) for n in range(1, 11)]
print(H1)
# [4, 95, 956, 9840, 99496, 999381, 9993490, 99980910, 999995621, 9999829204]

Fiindcă am greșit la scriere un semn, cum se mai întâmplă — trebuia "S1 + D", nu "S1 - D" — nu am obținut rezultatele $H_1$ așteptate (v. [1]), dar… recunoaștem imediat că ne-au rezultat pe neașteptate valorile $H_2$ (cu unele scăpări la ultima cifră, v. [1]).
Așa ne-am dat seama, că nu sunt necesare două formule (una pentru $H_1$ și alta pentru $H_2$); ajunsesem la $H_1$ urcând din capătul stâng al intervalului $\mathcal{I}_n$, apoi la $H_2, coborând din capătul drept (v. [1]) — dar în principiu, a coborâ dintr-un capăt este echivalent cu a urca din celălalt (desigur… ar fi de meditat: cum de n-am observat aceasta, decât experimentând și având noroc să mai și greșim câte un semn? E drept că nici în [2] nu s-a observat, că formula prin care se estimează „cel mai mic” poate estima — schimbând un semn — și pe „cel mai mare”).

Rescriem funcția de mai sus, vizând acum ambele cazuri (și angajând gmpy2):

import gmpy2 as gy
gy.set_context(gy.ieee(128))  # Quadruple-precision floating-point format

def estim12(n, SL = True):  # Smallest (H1), or Largest (H2)
    E0 = gy.exp10(n-1) if SL else gy.exp10(n) 
    LG = gy.log(10)
    delta = gy.frac(gy.lngamma(E0 + 1) / LG)
    D = gy.sqrt(1 + 8*E0*(1-delta)*LG) / 2
    res = gy.round2(E0 + D - 0.5) if SL else gy.round2(E0 - D - 0.5)
    return gy.mpz(res)

H1 = [estim12(n) for n in range(1, 24)]  # [mpz(3), mpz(14), mpz(103), ...]
H2 = [estim12(n, False) for n in range(2, 21)]  # [mpz(96), mpz(957), ...]

Pentru a afișa rezultatele în forma obișnuită, când va fi cazul, vom folosi:

def eval_frm(Lst):
    return [eval(format(gy.mpz(i))) for i in Lst]

print("H1 =", eval_frm(H1));  print("H2 =", eval_frm(H2))
H1 = [3, 14, 103, 1042, 10158, 100502, 1000617, 10006508, 100019088, 1000004377, 10000170793, 100000442970, 1000001981666, 10000005339905, 100000018997256, 1000000065392525, 10000000201014294, 100000000631520633, 1000000001077349719, 10000000005309690910, 100000000018130580501, 1000000000032533084356, 10000000000202461315940]
H2 = [96, 957, 9841, 99497, 999382, 9993491, 99980911, 999995622, 9999829206, 99999557029, 999998018333, 9999994660094, 99999981002743, 999999934607474, 9999999798985705, 99999999368479366, 999999998922650280, 9999999994690309089, 99999999981869419498]

Rezultatele pentru $H_1$ coincid (cu doar două excepții, pe cifrele boldate mai sus) cu cele din [2]; iar rezultatele $H_2$ obținute în [1] pentru $n \le 10$, se regăsesc în lista $H_2$ de mai sus (cu vreo două diferențe, pe ultima cifră).

În scopul de a găsi valoarea $H_1$ corectă, plecând de la estimarea dată de formula de calcul respectivă, în programul $\mathbf{C}$ din [2] este aplicată cu finețe o tehnică (specifică pentru biblioteca $MPFR$) de rotunjire în sus și în jos, îmbinată cu „metoda înjumătățirii” — și ne-am tot documentat o bună bucată de timp, căutând să imităm (pentru valorile $H_2$) calea arătată; dar la un moment dat… ne-am dat seama că totuși, nu este nevoie să procedăm astfel.

Formula de estimare pe care am stabilit-o în [1] (parcă mai bună decât cea din [2], fiindcă ia în seamă și termenii de gradul doi din dezvoltarea lui $\ln(1 + x)$) ne dă fie chiar valoarea corectă ($H_1$, respectiv $H_2$), fie (în unele cazuri) pe cea vecină, deasupra sau dedesubt, acesteia (anume, în acele cazuri în care "$\sum\delta_h$" — cum însemnam în [1] — ajunge foarte aproape de 1, încât rotunjirea specifică reprezentării „floating point” poate rata cu $\pm 1$ valoarea corectă). Nu-i cazul atunci (sau nu prea ar fi cazul), să invocăm „metoda înjumătățirii”…

Ne putem baza pe formula lui Stirling (într-o formă sau alta), pentru a determina când vom avea acuși nevoie, numărul de cifre ale factorialului:

e = gy.exp(1)
pi = gy.const_pi()
def numdig_Stirling(n):
    return gy.mpz(gy.ceil(n*gy.log10(n/e) + (gy.log10(2*n*pi))/2))  # +
                                          # gy.log10(1 + 1/(12*n)))

Inițial am folosit direct return gy.floor(gy.lngamma(n+1)/gy.log(10))+1 — dar așa (și nu mai căutăm acum, explicații) în dicționarul pe care-l definim mai jos am obținut doar "error" pentru $n=23$ (în timp ce, adoptând formula lui Stirling, am găsit $H_1$ corect și pentru $n \ge 23$).

În următoarea funcție se testează valorile $h$ din lista indicată (H1, sau H2): pentru fiecare $h$ se determină numărul de cifre din $h!$ și din $(h+1)!$ și se compară diferența cu $n$ (unde $n=i+1$, $i\ge 0$ fiind indexul valorii $h$ în lista respectivă); dacă diferența respectivă nu este egală cu $n$, atunci se calculează numdig_Stirling(h-1) și respectiv (în celălalt caz) numdig_Stirling(h+2) — valoarea corectă este fie $h-1$, fie $h+1$ (și afară poate, de cazul „teoretic” când $n$ chiar nu ar mai fi abordabil, implicând valori prea mari pentru a mai putea fi verificate — nu avem a ne teme că estimarea $h$ este mai departe decât cu $\pm 1$, față de valoarea corectă).
Și se returnează un dicționar care asociază fiecărui $n$ la care estimarea din lista indicată trebuie corectată (cu $\pm 1$), câte o listă conținând: valoarea corectată, numărul de cifre ale factorialului acesteia și numărul de cifre ale factorialului succesorului ei (pentru a constata mintal că diferența este $n$); dar în cazul de excepție, când pentru un anumit $n$ calculele eșuează dintr-un motiv sau altul, pe cheia $n$ înscriem doar "error":

def to_correct(H12):
    R = {}
    for i in range(0, len(H12)):
        h = H12[i]
        nd1 = numdig_Stirling(h)
        nd2 = numdig_Stirling(h + 1)
        sig = gy.cmp(nd2 - nd1, i + 1)
        if(sig != 0):
            Prev = numdig_Stirling(h-1)
            Next = numdig_Stirling(h+2)
            if(nd1 - Prev == i+1):
                R[i+1] = eval_frm([h-1, Prev, nd1])
            elif(Next - nd2 == i+1):
                R[i+1] = eval_frm([h+1, nd2, Next])
            else: 
                R[i+1] = "error"
    return R

dH1 = to_correct(H1);  print("corecturi în H1:", dH1, sep="\n")
dH2 = to_correct(H2);  print("\ncorecturi în H2:", dH2, sep="\n")

corecturi în H1:
{8: [10006509, 65702623, 65702631], 
23: [10000000000202461315939, 215657055185421630674159, 215657055185421630674182]}

corecturi în H2:
{1: [95, 149, 150], 
8: [999995621, 8565666113, 8565666121], 
14: [999999934607475, 14565704537208883, 14565704537208897]}

Corecturile trebuie aplicate prin constructorul mpz(); de exemplu pentru H1:

H1[7] = gy.mpz(dH1[8][0]);  H1[22] = gy.mpz(dH1[23][0]) 

Lista $H_1$ conține acum exact primele 23 de valori date în [1]; prin analogie, putem considera că și cale 19 valori $H_2$ sunt corecte (fiind produse printr-un același program).

… nu ne dă pace întrebarea filozofică iscată mai sus: cum de n-am observat „la timp”, că $H_1$ și $H_2$ se pot găsi (sau estima) printr-o aceeași formulă? Și parcă ar fi nevoie și de o justificare „teoretică” (chiar dacă este simplă) a faptului că $H_2$ se poate estima „urcând” de la $10^n$ (în loc de a coborâ de la $10^{n+1}-1$, ca în [1]).

… altfel (dar tot filozofic), trebuie să recunoaștem că în sine — ignorând (să zicem că se poate) aparatul investigativ folosit (fel de fel de mici raționamente matematice, experimente și cunoștințe de "floating-point", limbaje, biblioteci și programare) — chestiunea tratată în [1] și mai sus, privitoare la numărul de cifre dintre factoriale, nu are în sine, vreo însemnătate (…veșnica întrebare: de ce ne-ar interesa „numărul de cifre”?); mai concis, elegant și memorabil, în formularea unuia dintre editorii de la OEIS: "This is quite contrived...".