Care este cel mai mare $H$ pentru care $(H+1)!$ are cu $n$ cifre mai mult ca $H!$ ?
Ca și în [1], coborâm prin go_down() de la $10^{n+1}$ – dar de data aceasta folosim C (v. GCC 11.4.0), întâi cu biblioteca math uzuală (pe tipul de numere long double):
#include <stdlib.h> #include <stdio.h> #include <math.h> long double num_dig_fac(long double n) { return floorl(lgammal(n+1) / M_LN10) + 1; } long double go_down(int n) { long double h = powl(10, n+1); while(num_dig_fac(h) - num_dig_fac(h-1) > n) h--; /* coboară până când diferența lungimilor este n */ return h - 1; } int main() { long double A[12], B[12], C[12]; for(int i=0; i < 12; i++) { A[i] = go_down(i+1); /* H */ B[i] = num_dig_fac(A[i]); /* lungimea lui H! */ C[i] = num_dig_fac(A[i]+1); /* lungimea lui (H+1)! */ } printf("%-18s %-16s %-16s %-4s\n", "H", "len H!", "len (H+1)!", "n"); for(int i=0; i < 12; i++) { printf("%-18.0Lf %-16.0Lf %-16.0Lf %d\n", A[i], B[i], C[i], i+1); } return 0; } vb@Home:~/24mar$ gcc -o a1 1.c -lm vb@Home:~/24mar$ ./a1 H len H! len (H+1)! n 95 149 150 1 957 2440 2442 2 9841 35025 35028 3 99497 454060 454064 4 999382 5562002 5562007 5 9993491 65611498 65611504 6 99980911 756417846 756417853 7 999995621 8565666113 8565666121 8 9999829205 95655347238 95655347247 9 99999556977 1056565678564 1056565678574 10 999998017067 11565681722909 11565681722920 11 9999994609320 125656985102136 125656985102148 12
Rezultatele $H$ sunt corecte numai până la $n=8$ – apoi sunt nesigure, fiindcă diferă la ultimele cifre, de cele obținute anterior în R și în scipy (v. [1]).
Sistemul standard de reprezentare în virgulă mobilă (și în particular, a valorilor de tip long double) se dovedește insuficient pentru a obține valorile $H$ corecte pentru $n\ge9$. Iar faptul că (nici acum) nu s-a furnizat vreo atenționare de depășire, ne arată că este necesară prevederea explicită în program a unui mecanism de supraveghere a calculelor (pentru a depista inexactitatea și măcar, a furniza ceva informații contextuale).
Probabil, calculul a angajat coprocesorul matematic (rezultând un timp de execuție foarte scurt, sub 3 secunde); prin biblioteca libquadmath putem folosi modelul (soft) de calcul pe 128 biți (în loc de regiștrii de 80 biți și instrucțiunile coprocesorului), acceptând desigur o creștere de peste zece ori, a timpului de execuție:
#include <quadmath.h> #include <stdlib.h> #include <stdio.h> __float128 num_dig_fac(unsigned long int n) { return floorq(lgammaq(n+1)/M_LN10q) + 1; } __float128 go_down(int n) { unsigned long int h = powq(10, n+1); for(h; num_dig_fac(h) - num_dig_fac(h-1) > n; --h); return --h; } int main() { __float128 A[13]; for(int i=0; i < 12; i++) A[i] = go_down(i+1); A[12] = num_dig_fac(9999994609321); /* v. H anterior, la n=12 */ char buf[128]; for(int i=0; i < 13; i++) { int h = quadmath_snprintf(buf, sizeof buf, "%.30Qg", A[i]); printf("%s\n", buf); } } vb@Home:~/24mar$ gcc -o a2 2.c -lquadmath vb@Home:~/24mar$ time ./a2 /* 0m33,268s (×11 față de "1.c") */ 95 957 9841 99497 999382 9993491 99980911 999995621 9999829206 /* pentru n=9 */ 99999557029 999998018333 9999994660094 125656985102149 /* len (H+1)! pentru H=9999994609320 */
Rezultatele păstrează caracteristicile evidențiate anterior; dar să observăm că acum (lucrând pe 128 de biți), pentru $n\ge 10$ se coboară mai puțin (uneori, mai mult) până ce se decide valoarea $H$. De exemplu, pentru $n=12$ cu "1.c" se găsea $H=99999946\boldsymbol{09320}$, cu semnificația că acesta este cel mai mare pentru care $(H+1)!$ are cu 12 cifre mai mult ca $H!$ – ori verificarea pentru acest $H$ făcută în al doilea program ("2.c") arată că de fapt, $(H+1)!$ are 13 cifre și nu 12, mai mult ca $H!$; al doilea program găsește în loc $H=99999946\boldsymbol{60094}$, adică mai sus cu vreo 5000 de locuri decât găsise "1.c".
Rezultatele din "2.c" sunt mai credibile decât cele din "1.c" (fiindcă în "2.c" se vizează mai multe cifre exacte, decât în "1.c"), dar pentru ambele programe – nu avem de ce să fim siguri de rezultatele produse (pentru $n\ge 9$), în absența unui control al valorilor logaritmice care intervin în calculul num_dig_fac(); când acestea ajung foarte apropiate (fie mai sus, fie mai jos) de o valoare întreagă, nu putem fi siguri de care parte a acesteia este rezultatul corect (este „13 cifre”, sau este „12 cifre”?).
O idee mai bună pare a fi angajarea bibliotecii matematice MPFR, care introduce o reprezentare în „virgulă mobilă” analogă celei standard, dar pe zone de memorie cât se poate de mari (încât calculul angajează un număr suficient, controlabil, de cifre exacte); iar pentru comoditate (putând viza noul tip de numere la fel ca tipurile obișnuite), putem folosi modelul asociat pentru C++, "mpreal.h" (v. MPFR C++):
#include <iomanip> #include <iostream> #include "mpreal.h" using mpfr::mpreal; using std::cout, std::setw, std::setprecision; mpreal num_dig_fac(mpreal n) { return floor(lgamma(n+1) / M_LN10) + 1; } mpreal go_down(int n) { mpreal h = pow(10, n+1); while(num_dig_fac(h) - num_dig_fac(h-1) > n) h--; return h - 1; } int main() { mpreal::set_default_prec(256); mpreal A[12], B[12], C[12]; for(int i=0; i < 12; i++) { A[i] = go_down(i+1); B[i] = num_dig_fac(A[i]); C[i] = num_dig_fac(A[i]+1); } for(int i=0; i < 12; i++) { cout << setw(3) << (i+1) << setw(20) << setprecision(19) << A[i] << setw(20) << setprecision(19) << B[i] << setw(20) << setprecision(19) << C[i] << '\n'; } return 0; } vb@Home:~/24mar$ g++ -o 3cc 3.cc -lmpfr vb@Home:~/24mar$ time ./3cc /* 5m20,709s */ 1 95 149 150 2 957 2440 2442 3 9841 35025 35028 4 99497 454060 454064 5 999382 5562002 5562007 6 9993491 65611498 65611504 7 99980911 756417846 756417853 8 999995621 8565666113 8565666121 9 9999829205 95655347238 95655347247 10 99999556977 1056565678564 1056565678574 11 999998017067 11565681722909 11565681722920 12 9999994609254 125656985101278 125656985101290
Pentru $n=10$ și $n=11$ s-au găsit aceleași valori ca în "1.c" (contrazicând însă "2.c"), încât acestea devin mai credibile.
Timpul de lucru este însă enorm, peste 5 minute pentru $n\le 12$, iar rezultatele, câte s-au putut obține, sunt totuși incerte… Ar trebui să evităm coborârea pas cu pas, în care verificăm fiecare număr în parte – încercând în loc, să estimăm cumva numărul de pași cu care ar trebui să coborâm de la $10^{n+1}$ pentru a ajunge într-o vecinătate, acceptabil de mică, a numărului $H$ pe care îl căutăm (urmând să explorăm cam ca în "3.cc" – în orice caz, folosind MPFR – numai numerele din acea vecinătate).
Estimarea numărului de pași cu care să coborâm necesită o analiză matematică prealabilă a lucrurilor (probabil mai simplă decât pare la prima vedere); apoi, probabil că trebuie să vedem mai precis, cam câte cifre exacte ar fi necesare (sau suficiente) pentru valorile logaritmilor implicați, pentru a ne încredința că valoarea calculată prin MPFR pentru $H$ este nici mai sus, nici mai jos decât cea corectă.
vezi Cărţile mele (de programare)