Se găsesc peste tot, programe care calculează factoriale utilizând direct definiția
$n!=(n-1)!n=n*(n-1)*\cdots*2*1$.
Avem și noi de prin 2007, o poveste [1] – care modelează calculul nu în baza 10, ci pe blocuri de câte 32 de biți; de exemplu, ?100000 produce $10^5!$ în mai puțin de 5 secunde, furnizând rezultatul ca șir de 397177 cifre hexazecimale "120CCAA20..."; dar obținerea reprezentării uzuale, în baza 10 (de exemplu, pentru a constata că $10^5!$ are 456574 cifre zecimale) durează mult mai mult: forma binară a factorialului (de utilizat direct în eventuale calcule ulterioare) este obținută prin limbajul de asamblare, iar transformarea în baza 10 (în scopul banal de a afișa rezultatul în forma uzuală) este lăsată limbajului javaScript.

Dar uneori (de fapt, cel mai adesea) chiar nu ne interesează valoarea exactă a factorialului, ci anumite caracteristici laterale – de exemplu, numărul de cifre (în baza uzuală, 10).
Și nu este necesar să calculăm, de exemplu $2^{1000}$, dacă ne interesează doar numărul de cifre: $\lg(2^{1000})=1000\lg 2\approx 1000*0.30103=301.03$ ne arată că $2^{1000}$ are 302 cifre (acceptând că $\lg 2$ este calculat cu suficientă acuratețe…).
$\lfloor \lg N \rfloor +1$, sau $\lceil \lg(N+1)\rceil$ – adică floor(log10(N))+1, respectiv ceiling(log10(N+1)) – ne dă numărul de cifre ale numărului natural $N$ („oricât” de mare permite sistemul de calcul).

Desigur, există veșnica întrebare, care în cazul de față ar fi: de ce ne-ar interesa numărul de cifre? … Cum-necum, pe OEIS găsim sute de secvențe (asociate câte unui anumit interes cognitiv) legate de "number of digits"; de exemplu, A034886 vizează chiar numărul de cifre din factoriale, A356758 prezintă numărul de cifre nenule din $n!$, iar A035065 – care numere au în factorial un număr prim de cifre.

Adăugăm și noi o întrebare; pentru $n=1, 2, 3, ...$ :
care este cel mai mare $K$ pentru care $(K+1)!$ are cu $n$ cifre mai mult decât are $K!$ ?

$4!$ este primul care are cu o cifră mai mult ca $3!$; dar $6!$ are același număr de cifre ca $5!$… Care este ultimul $n$ pentru care $(n+1)!$ are cu o singură cifră mai mult ca $n!$ ? Putem face o primă investigație astfel (speculând orientarea vectorială specifică limbajului R):

Df <- factorial(1:170) |> log10() |> ceiling() |> diff()

factorial() determină aproximativ valorile factorialelor numerelor din vectorul primit ca argument; prin log10() și ceiling() găsim numărul cifrelor din factorialele respective; aplicând diff() pe vectorul rezultat – aflăm diferențele dintre numărul de cifre din $(n+1)!$ și din $n!$ (precizăm că operatorul "pipe" |> transmite spre dreapta rezultatul obținut în stânga sa).

Frecvența valorilor din vectorul Df este:

table(Df)
  0   1   2   3 
  2  36 122   9 

deci între primele 170 de numere, avem 36 de numere $n$ cu proprietatea că $(n+1)!$ are cu o singură cifră mai mult ca $n!$ (dar… v. Obs. mai jos); aflăm pe cel mai mare dintre acestea:

wh <- which(Df == 1)  # indecșii la care avem valoarea 1
L <- length(wh)  # 36 de indecși (în ordine crescătoare)
print(wh[L])  # 95

Putem și verifica (prin [1], de exemplu) că $95!$ are 149 de cifre și $96!$ are cu una mai mult, 150 de cifre; iar 95 este cel mai mare cu această proprietate, fiindcă este ultimul index la care în Df avem valoarea 1, iar depășind 100=$10^2$ – numărul de cifre ale factorialului crește cu cel puțin 2.

Obs. Pentru Df=1 sunt 35, nu 36. Eroarea provine din faptul că am folosit ceiling(log10(N)), în loc de ceiling(log10(N + 1)) – ceea ce „merge”, dar numai dacă N nu este o putere a lui 10…

Subliniem că urmând schema de mai sus, n-am putea investiga mai departe de 170:

> factorial(170)
[1] 7.257416e+306  # deci 170! are 307 cifre
> factorial(171)
[1] Inf  # atenționare că 171! nu mai poate fi reprezentat

Putem constata cu [1], că $171!$ are 310 cifre – încât nu a mai putut fi reprezentat: IEEE_754 asigură reprezentarea în virgulă mobilă pentru numere care au cel mult 308 cifre.

Ar fi două posibilități: 1) folosim o bibliotecă pentru lucrul cu „numere mari” (de exemplu, pachetul R gmp oferă – via GMP – funcția factorialZ(), prin care putem obține valorile exacte ale factorialelor, urmând să înscriem lungimile acestora într-un vector, căruia să-i aplicăm diff(), ca în secvența de mai sus); 2) determinăm direct numărul de cifre, evitând calculul factorialelor. Optăm deocamdată pentru varianta 2); vom constata eventual, că astfel putem investiga în siguranță „doar” până undeva spre $10^8$ (și dacă va fi cazul să investigăm mai sus de atât, atunci vom putea extinde mai departe folosind prima variantă).

Avem $\lg n!=\sum_{k=1}^n \lg k$, dar s-ar cuveni ca suma să fie calculată nu vectorial, cu sum(log10(1:n)) – ci din aproape în aproape (pentru a evita recalcularea de sume):

N <- 10^7
D <- vector("numeric", N)
sl <- 0
for(n in 1:N) {
    sl <- sl + log10(n)
    D[n] <- floor(sl) + 1
}

În D[n] avem acum numărul de cifre ale lui $n!$, pentru $n\le 10^7$ (ignorăm suspiciunile specifice calculului aproximativ – după care, prin cumularea erorilor de aproximare s-ar ajunge la un D[n] care să nu fie corect, sau reprezentabil); salvând pe disk (prin saveRDS()), obținem un fișier de 23.8 MiB (cam mare – încât am ezitat să trecem de $10^7$).

Obținem diferențele dintre D[n+1] și D[n], pentru n de la $1$ la $10^7-1$:

Df <- diff(D)

table(Df) ne dă:

         0       1       2       3       4       5       6       7 
         3      35     351    3518   35178  351779 3517784 6091351 

Deci avem 35 de numere pentru care Df este 1, apoi 351 pentru Df=2, ș.a.m.d.; dar pentru a vedea câte numere au Df=7 ar trebui să căutăm mai departe, dincolo de $10^7$.
Obs. Ar fi de bănuit (și poate, de conjecturat), că numărul de valori pentru care Df=d are ordinul de mărime $3.5*10^d$ – ceea ce se vede mai sus pentru $d\le 6$.

Următoarea secvență determină cel mai mic și cel mai mare număr, pentru care Df are valoarea $n\le 6$ (și pe cel mai mic pentru care Df=7):

S <- L <- vector("numeric", 7)
for(n in 1:7) {
    wh <- which(Df == n)
    cnt <- length(wh)
    S[n] <- wh[1]
    L[n] <- wh[cnt]
}
print(S)  # Smallest 3      14     103    1042   10158  100502 1000617
print(L)  # Largest 95     957    9841   99497  999382 9993491 9999999

Căutând pe OEIS prima secvență 3,14,103,1042 — aflăm A084420, "Least number k such that between k! and (k+1)! there are n powers of ten". Să verificăm puțin: $3!=6$ și $4!=24$ au între ele o singură putere a lui 10; $14!$ are 11 cifre, deci este cuprins între $10^{10}$ și $10^{11}$, iar $15!$ are 13 cifre – deci între $14!$ și $15!$ avem exact două puteri ale lui 10, $10^{11}$ și $10^{12}$ (dar este mai greu de verificat că 14 este cel mai mic cu această proprietate).

Inspirându-ne de aici, probabil că putem reformula întrebarea noastră astfel: cel mai mare $k$ pentru care între $k!$ și $(k+1)!$ avem exact $n$ puteri ale bazei 10.
Așa este! $k$ este cel mai mare pentru care Df este $n$ (adică $(k+1)!$ are $n$ cifre mai mult decât are $k!$) numai dacă deasupra lui $k$ găsim undeva $10^{n+1}$ (mai încolo de această valoare, numărul de cifre crește cu cel puțin $(n+1)$, deci $k$ este ultimul cu proprietatea cerută).

Acum problema noastră este următoarea: cum putem scăpa de calcularea în prealabil a tabelului (prea voluminos) D (cu numărul de cifre din factoriale)? Dacă putem evita D, vom putea trece de $10^7$ și eventual, vom putea obține și pentru $n>6$, valorile întrebate.

De unde să vină răspunsul? Întâi de la consola R:

> help(factorial)  # aflăm de lfactorial(), etc. 
> factorial
function (x) 
gamma(x + 1)  # deci factorial(x) = gamma(x+1)
<bytecode: ... >
<environment: namespace:base>

lfactorial(x) aproximează (cu acuratețe) logaritmul natural al lui x (împărțim cu $\ln 10$, pentru a avea logaritmul zecimal) – ceea ce ne scutește de a mai modela însumarea logaritmilor, încât putem proceda foarte direct (și mergem acum până la $10^8$ și chiar mai departe):

S <- diff(floor(lfactorial(1:(10^8))/log(10))+1)
for(n in 1:7)
    print(tail(which(S == n), 1))
[1] 95
[1] 957
[1] 9841
[1] 99497
[1] 999382
[1] 9993491
[1] 99980911

Pe lângă cele 6 valori obținute anterior, avem acum răspunsul și pentru $n=7$.

Un mare avantaj față de versiunea anterioară constă în faptul că acum putem investiga ușor „pe bucăți”; găsim răspunsul pentru $n=8$ în intervalul $[9*10^8+1,\, 11*10^8]$:

S <- diff(floor(lfactorial((9*10^8+1):(11*10^8))/log(10))+1)
print(tail(which(S == 8), 1))  # 99995621

Am observat mai sus că valoarea căutată $k$ este ultima dinaintea lui $10^{n+1}$, care mai are proprietatea cerută – $(k+1)!$ are cu $n$ cifre mai mult decât are $k!$ – iar aceasta conduce la ideea de a pleca de la $(10^{n+1})!$ = S și a coborî cât timp diferența dintre numărul de cifre din S și cel din factorialul curent se păstrează la un multiplu de $(n+1)$.
De exemplu, pentru $n=1$ avem:

(10^2)! are 158 cifre. Coborâm:
-1 (99!)    156 (cu 1×2 mai puțin ca 100!)
-2          154 (cu 2×2 mai puține ca 100!)
-3          152 (cu 3×2 mai puține ca 100!)
-4          150 (cu 4×2 mai puține ca 100!)
-5 (95!)    149  o singură cifră mai puțin ca 96!

și trebuie să ne oprim după cei 5 pași, fiindcă $95!$ are 149 de cifre (și nu cu 5×2 mai puțin ca $100!$); rezultă că 95 este cel mai mare (socotind desigur de la 1 în sus) al cărui factorial are cu o singură cifră (și nu cu două) mai puțin ca factorialul următorului ($96!$).

Ajungem deci la această modelare, prin care se reduce drastic consumul de memorie specific versiunilor anterioare și ne găsește rapid valorile căutate, până la $n=10$:

A <- vector("numeric", 10)
for(n in 1:10) {
    S <- 10^(n+1)
    fs <- floor(lfactorial(S)/log(10)) + 1
    i <- 1
    while(fs - floor(lfactorial(S - i)/log(10)) - 1 == i*(n+1)) 
        i <- i + 1
    A[n] <- S - i
}
print(A)
 [1]          95         957        9841       99497      999382     9993491
 [7]    99980911   999995621  9999829206 99999557080

Putem extinde și mai sus de $n=10$, dar obținem numai valori aproximative – de exemplu, pentru $n=11$ obținem 9.99998e+11 (avem numai primele cifre 999998, dintre cele 12); totuși intern, valoarea respectivă (pentru $n=11$) este reprezentată exact – cum vedem prin format(A[11], digits=16), care afișează 999998019267.

Putem verifica rezultatele redate mai sus, folosind GMP: calculăm factorialul lui $A[n]$ (chit că ajunge la câteva sute de milioane de cifre) și determinăm numărul de cifre, apoi înmulțim factorialul respectiv cu $(A[n]+1)$ (obținând $(A[n]+1)!$ ) și determinăm numărul de cifre ale rezultatului – trebuie să fie cu $n$ mai mare ca în cazul lui $(A[n])!$ ). De exemplu, la $n=7$:

> a7 <- 99980911
> fa7 <- gmp::factorialZ(a7)  # foarte rapid!
> nchar(as.character(fa7))  # socotirea numărului de cifre ia ceva timp...
[1] 756417846  # numărul de cifre din (a7)!
> fa71 <- gmp::mul.bigz(fa7, a7+1)  # factorialul lui (a7+1)
> nchar(as.character(fa71))
[1] 756417853  # cu 7 cifre mai mult ca (a7)!

Însă pentru $n=8$ numărul de cifre ale factorialului devine deja prea mare, încât (pe sistemul meu) n-am mai putut confirma astfel (prin gmp::factorialZ()) valorile A[n] cu n≥8; dar ideea de a chinui astfel calculatorul (pentru a obține cifrele factorialului pentru $n=8$, în număr de ordinul miliardelor) este chiar oribilă! Iar adevărul este că nu prea putem investiga mai departe de $n=8$, fiindcă ajungem să depășim posibilitățile uzuale de reprezentare a numerelor (cel puțin în cazul de aici, când folosim R).

Revăzând (din zi în zi, ca de obicei) cele de mai sus, observăm în final că puteam angaja din capul locului doar formula de bază $\lg n!\approx\sum_1^n \lg k$; nu era necesar să determinăm întâi vectorul uriaș D – de ordinul a 23.6 MiB, conținând numărul de cifre pentru fiecare $n!$ – și să folosim apoi diff(), pentru a găsi în sfârșit, vectorul Df al diferențelor D[n+1]-D[n] – cum avem în primul program de mai sus, până $10^7$.
Este chiar foarte simplu să obținem direct, diferențele respective – și anume, într-un vector de tip integer (ceea ce scurtează la numai vreo 300 Kib volumul de memorie necesar și aduce sub 10 secunde timpul necesar obținerii vectorului Df):

N <- 10^7
Df <- vector("integer", N)
slg <- s2 <- 0
for(n in 2:N) {
    s1 <- s2
    slg <- slg + log10(n)
    s2 <- floor(slg)
    Df[n] <- as.integer(s2 - s1)
}
Df <- Df[2:length(Df)]
saveRDS(Df, "diff10_7.RDS")  # 302.6 KiB
print(table(Df))
#      0       1       2       3       4       5       6       7 
#      3      35     351    3518   35178  351779 3517784 6091351 
S <- L <- vector("numeric", 7)
for(n in 1:7) {
    wh <- which(Df == n)
    cnt <- length(wh)
    S[n] <- wh[1]
    L[n] <- wh[cnt]
}
print(S)  # Smallest:  3      14     103    1042   10158  100502 1000617
print(L)  # Largest:   95     957    9841   99497  999382 9993491

Iar al doilea program de mai sus – în care coboram de la $10^{n+1}$ – se poate formula mai bine folosind una sau alta dintre formulele de calcul redate în A034886:

num_dig_fac <- function(n)  # după formula clasică a lui Stirling
    ceiling((log(2*pi) - 2*n + (1+2*n)*log(n))/(2*log(10)))  

Largest <- function(n) {
    h <- 10^(n+1)
    Zd <- num_dig_fac(h)
    i <- 1
    repeat{
        hi <- h - i
        Zi <- num_dig_fac(hi)
        if(Zd - Zi < i*(n+1))
            break
        i <- i + 1
    }
    hi    
}

L <- sapply(1:8, Largest)
#[1] 95  957  9841  99497  999382  9993491  99980911  999995620

Rezultatul redat se obține într-o fracțiune de secundă; dar să observăm că al 8-lea termen diferă cu 1 față de cel obținut prin programele anterioare – ceea ce înseamnă că nu putem fi siguri de rezultat decât pentru $n\le 7$ (pentru care am găsit aceleași valori utilizând formule de calcul diferite și pentru care dispunem de diverse posibilități de verificare).