În [1] am folosit R într-o chestiune care angaja numere mari și am constatat că rezultatele devin incerte pe măsură ce se depășește $2^{31}-1$; am căuta motive și corecții, pe vreo trei direcții (care de care mai instructive): de elucidat „precizia” de calcul asigurată de R și să zicem, de Python (că poate n-am întrebat pe cine trebuie); de angajat și alte formule de calcul decât cele două vizate în [1] și a confrunta rezultatele; a adapta programul pentru a putea confrunta (prin analogie) rezultatele respective, cu unele cunoscute.

Deocamdată, vedem că (fiind instalat automat – fără ajustări suplimentare – dintr-o distribuție Ubuntu) R ne asigură (privind deocamdată din afară) o precizie de 8-10 cifre semnificative; în schimb sub Python putem avea până la 15-17 cifre exacte:

vb@Home:~$ R  # version 4.3.2; Platform: x86_64-pc-linux-gnu (64-bit)
> 2^31-1  #[1] [1] 2147483647
> 2^63-1  #[1] 9.223372e+18
> gamma(171)  # 170!   (n! = Γ(n+1))  #[1] 7.257416e+306

vb@Home:~$ python  # Python 3.10.12; [GCC 11.4.0] on linux
>>> 2**63-1  # 9223372036854775807
>>> from math import gamma
>>> gamma(171)  # 7.257415615307998e+306

Pentru „calculul științific” se folosește de obicei Python, împreună cu diverse biblioteci de calcul numeric (alături de C și limbaje de asamblare).
Abordăm acum ceva mai riguros chestiunea din [1] și o investigăm folosind sciPy.

Care este cel mai mare $H$ (natural) pentru care $(H+1)!$ are cu $n$ cifre (una, două, trei, etc.) mai mult decât are $H!$ ?

Prin logaritmare în baza 10 (și vizând apoi partea întreagă), condiția pusă revine la $\lfloor\lg(H+1)!\rfloor + 1 \boldsymbol{=} \lfloor\lg(H!)\rfloor + 1 + n \iff \lfloor\lg(H!) + \lg(H+1)\rfloor = \lfloor\lg(H!)\rfloor+ n$; rezultă că $\lfloor\lg(H+1)\rfloor \boldsymbol{\le} n$. Dar $\lfloor\lg(H+1)\rfloor > \lg(H+1)-1$, deci: $\lg(H+1) \boldsymbol{<} n+1$.

Prin urmare, $H+1 \boldsymbol{< 10^{n+1}}$; deci pentru fiecare $n$ dat, putem găsi $H$ coborând pas cu pas de la $10^{n+1}$, până când valoarea curentă are factorialul cu $n$ cifre mai lung decât factorialul valorii aflate dedesubtul ei:

from math import floor, log
from scipy.special import gammaln

def num_dig_fac(n):    
    return floor(gammaln(n+1)/log(10)) + 1

def go_down(n):
    h = 10**(n+1)
    n1 = num_dig_fac(h)
    while(True):
        h = h - 1
        n2 = num_dig_fac(h)
        if(n1 - n2 == n):  # n1! are cu n cifre mai mult ca n2!
            return h
        n1 = n2  # evită un al doilea apel num_dig_fac()

A = [go_down(i) for i in range(1, 13)]
print(A)
# [95, 957, 9841, 99497, 999382, 9993491, 99980911, 999995621, 
   9999829206, 99999557080, 999998019267, 9999994765612]

Primele 11 rezultate coincid respectiv, cu cele obținute, folosind R, în [1] (unde le-am și verificat direct, până la $n$=7).

Adaptăm acum funcția redată mai sus, încât să putem genera într-un același context, niște rezultate cunoscute; dacă în loc să coborâm de la $10^{n+1}$, urcăm de la $10^{n-1}$ (care este primul cu $n$ cifre) – găsim cel mai mic număr al cărui factorial este cu $n$ cifre mai scurt decât factorialul numărului următor (deasupra) lui:

def go_up(n):
    h = 10**(n-1)
    n1 = num_dig_fac(h)
    while(True):
        h = h + 1
        n2 = num_dig_fac(h)
        if(n2 - n1 == n):  # n2! are cu n cifre mai mult ca n1!
            return h - 1
        n1 = n2

A = [go_up(i) for i in range(1, 13)]
print(A)
# [3, 14, 103, 1042, 10158, 100502, 1000617, 10006509, 
   100019088, 1000004377, 10000170792, 100000442933]

Rezultatele obținute acum sunt (în principiu) cunoscute: căutând pe OEIS secvența "3 14 103 1042", găsim A084420 – cu o formulare care este echivalentă cu a noastră: dacă $H!$ are $\boldsymbol{q}$ cifre și $(H+1)!$ are cu $n$ cifre mai mult (formularea noastră), atunci $H!\boldsymbol{<10^q}$ și $(H+1)!\boldsymbol{>10^{n+q-1}}$ deci (formularea din A084420) între $H!$ și $(H+1)!$ există exact $n$ puteri ale bazei 10.
Deosebirea ține de analogie: A084420 – și funcția de mai sus go_up() – indică pe cel mai mic, iar noi – în formularea inițială go_down() – indicam pe cel mai mare, $H$.

Cele două funcții constituite mai sus folosesc o aceeași formulă de calcul și diferă numai prin sensul în care o aplică (de sus în jos, sau invers); prin urmare, dacă rezultatele din go_up() sunt corecte, atunci este firesc să considerăm că și acelea din go_down() sunt deasemenea, corecte; iar acum putem valida rezultatele din go_up(), prin confruntare cu A084420 – unde găsim aceste valori:

[3, 14, 103, 1042, 10158, 100502, 1000617, 10006509,
 100019088, 1000004377, 10000170793, 100000442970,
 1000001981666, 10000005339905]

Constatăm că rezultatele coincid până la $n$=10 (și le-am boldat mai sus, în ambele cazuri); dar pentru $n$=11, rezultatul din go_up() diferă la ultima cifră ("2" și respectiv "3") față de A084420, iar pentru $n$=12 – diferă la ultimele două cifre ("33", respectiv "70").
Se pare că problema noastră de validare se reduce la: cine a greșit?…

În pagina secvenței A084420 este redat acest program, pentru Wolfram Mathematica (aici îl formatăm cumva pe linii):

LogBase10Stirling[n_] := Floor[ Log[10, 2 Pi n]/2 + n * Log[10, n/E] + 
    Log[10, 1 + 1/(12n) + 1/(288n^2) - 139/(51840n^3) - 571/(2488320n^4) + 
        163879/(209018880n^5)]];
Do[
    k = 10^(n - 1); 
    While[ LogBase10Stirling[k + 1] - LogBase10Stirling[k] < n, k++ ]; 
    a[n] = k; 
    Print[ a[n]], 
{n, 1, 14}]

Ca și în go_up(), se urcă pas cu pas de la $10^{n-1}$, până când diferența de cifre dintre factorialele a două valori consecutive este $n$ (doar că acum… se calculează de câte două ori un același LogBase10Stirling[], o dată ca scăzător și a doua oară ca descăzut, ceea ce – exceptând cazul improbabil că ar fi vorba de vreo subtilitate de limbaj sau de calcul numeric – este totuși, un exemplu clasic de „rele practici”).

În fond (ca funcționare), cele două programe sunt identice și până la $n$=10, produc același rezultat – și deducem prin analogie, că rezultatele din go_down() (obținute și în [1] prin R) sunt corecte pentru $n\boldsymbol{< 10}$ (valoarea pentru $n=10$ rămâne incertă: în go_down() se coboară de la o valoare foarte mare și „analogia” cu go_up() devine discutabilă).
Numai începând de la $n=11$, ele dau răspunsuri diferite – și de fapt, acestea (valori de ordinul $10^{12}$ și $10^{13}$) diferă numai la ultima sau la ultimele două cifre.

Diferența arătată nu decurge – cum am fi tentați să credem – din calculul „diferit” al logaritmului factorialului: în go_up() folosim direct scipy.special.gammaln(), iar programul din A084420 angajează primii 6 termeni din seria lui Stirling (v. de exemplu Stirling's_approximation); se poate constata consultând codul expus la github.com/scipy, că gammaln() implică deasemenea, o aproximare polinomială a seriei lui Stirling.

Deci întrebarea justă ar fi nu „cine a greșit”, ci unde (pe ce sistem de calcul) rulează fiecare dintre cele două programe comparate mai sus; cel mai probabil, pe calculatorul pe care s-a rulat LogBase10Stirling[] este asigurată o precizie de calcul extinsă ("128-bit floating-point values"; v. Quadruple-precision), în timp ce pe sistemul propriu dispunem de precizia obișnuită (double, sau "binary64").

Ne-am făcut astfel „de lucru”… Ce este aceea „precizie” și ce sunt numerele "floating point" (dincolo de ce știm desigur, pentru un examen teoretic sumar)? Care sunt valorile $H$ corecte, pentru $n$=10, 11 și 12 (și cum ajungem la ele pe sistemul propriu)? Nu mai zicem și de teoria funcțiilor speciale (pentru $\Gamma(z)$ și $\ln\Gamma(z)$)…