Difference between revisions of "Psihologia concursurilor de informatică/2 Lucrul cu numere mari"

From Algopedia
Jump to: navigation, search
(Înmulțirea a doi vectori)
(Înmulțirea a doi vectori)
Line 277: Line 277:
  
 
care duce prin eliminarea recurenței la <math>T(N) \in O(N^2)</math>. Cu alte cuvinte, încă n-am câștigat nimic. Trebuie să reușim cumva să reducem numărul de înmulțiri de la 4 la 3, chiar dacă prin aceasta vom mări numărul de adunări necesare. Să definim produsul
 
care duce prin eliminarea recurenței la <math>T(N) \in O(N^2)</math>. Cu alte cuvinte, încă n-am câștigat nimic. Trebuie să reușim cumva să reducem numărul de înmulțiri de la 4 la 3, chiar dacă prin aceasta vom mări numărul de adunări necesare. Să definim produsul
 
  
Atunci putem scrie:
+
<math>G = (C + D) \times (E + F) = CE + CF + DE + DF = CE + DF + (CF + DE)</math>
 +
 
 +
Atunci putem scrie:
  
 +
<math>A \times B = CE \times 10^{N} + (G - CE - DF) \times 10^{N / 2} + DF</math>
 
 
 +
Pentru această variantă, folosim doar trei înmulțiri, și chiar dacă avem nevoie de șase adunări și scăderi și două înmulțiri cu puteri ale lui 10, complexitatea se va reduce la <math>O(N^{\log_2 3})</math>. În cazul în care '''''N''''' este o putere a lui 2, împărțirea în două a numerelor se poate face fără probleme la fiecare pas, până se ajunge la numere de o singură cifră, care se înmulțesc direct. În cazul în care '''''N''''' nu este o putere a lui 2, este comod să se completeze numerele cu zerouri până la o putere a lui 2. În funcțiile descrise mai jos, <tt>MultRec</tt> nu face decât înmulțirea recursivă, pe când <tt>MultVect2</tt> se ocupă și de corectarea numărului de cifre (incrementarea până la o putere a lui 2). Pentru calculul produselor <math>C \times E</math> și <math>D \times F</math>, procedura <tt>MultRec</tt> se autoapelează; pentru calcularea produsului <math>(C+D) \times (E+F)</math>, însă, este nevoie să fie apelată procedura <tt>MultVect2</tt>, deoarece prin cele două adunări poate să apară o creștere a numărului de cifre al factorilor, care în acest caz trebuie readuși la un număr de cifre egal cu o putere a lui 2.
  
Pentru această variantă, folosim doar trei înmulțiri, și chiar dacă avem nevoie de șase adunări și scăderi și două înmulțiri cu puteri ale lui 10, complexitatea se va reduce la . În cazul în care N este o putere a lui 2, împărțirea în două a numerelor se poate face fără probleme la fiecare pas, până se ajunge la numere de o singură cifră, care se înmulțesc direct. În cazul în care N nu este o putere a lui 2, este comod să se completeze numerele cu zerouri până la o putere a lui 2. În funcțiile descrise mai jos, MultRec nu face decât înmulțirea recursivă, pe când MultVect2 se ocupă și de corectarea numărului de cifre (incrementarea până la o putere a lui 2). Pentru calculul produselor C x E și D x F, procedura MultRec se autoapelează; pentru calcularea produsului (C+D) x (E+F), însă, este nevoie să fie apelată procedura MultVect2, deoarece prin cele două adunări poate să apară o creștere a numărului de cifre al factorilor, care în acest caz trebuie readuși la un număr de cifre egal cu o putere a lui 2.
+
<syntaxhighlight lang=c>
 
 
 
 
 
void MultHuge2(Huge A, Huge B, Huge P);
 
void MultHuge2(Huge A, Huge B, Huge P);
  
 
void MultRec(Huge A, Huge B, Huge P)
 
void MultRec(Huge A, Huge B, Huge P)
Î Huge C,D,E,F,CE,DF;
+
{ Huge C,D,E,F,CE,DF;
  
   if (Aî0ș==1)
+
   if (A[0]==1)
     Î Pî1ș=Aî1ș*Bî1ș;
+
     { P[1]=A[1]*B[1];
       Pî0ș=(Pî2ș=Pî1ș/10)>0 ? 2 : 1;
+
       P[0]=(P[2]=P[1]/10)>0 ? 2 : 1;
       Pî1ș%=10;
+
       P[1]%=10;
     Ș
+
     }
     else Î Pî0ș=0;
+
     else { P[0]=0;
           AtribHuge(C,A);Shr(C,Aî0ș/2);
+
           AtribHuge(C,A);Shr(C,A[0]/2);
           AtribHuge(D,A);Dî0ș=Aî0ș/2;
+
           AtribHuge(D,A);D[0]=A[0]/2;
           AtribHuge(E,B);Shr(E,Bî0ș/2);
+
           AtribHuge(E,B);Shr(E,B[0]/2);
           AtribHuge(F,B);Fî0ș=Bî0ș/2;
+
           AtribHuge(F,B);F[0]=B[0]/2;
 
           MultRec(C,E,CE);MultRec(D,F,DF);
 
           MultRec(C,E,CE);MultRec(D,F,DF);
 
           Add(C,D);Add(E,F);
 
           Add(C,D);Add(E,F);
 
           MultHuge2(C,E,P);
 
           MultHuge2(C,E,P);
 
           Subtract(P,CE);Subtract(P,DF);
 
           Subtract(P,CE);Subtract(P,DF);
           Shl(P,Aî0ș/2);
+
           Shl(P,A[0]/2);
           Shl(CE,Aî0ș);Add(P,CE);
+
           Shl(CE,A[0]);Add(P,CE);
 
           Add(P,DF);
 
           Add(P,DF);
         Ș
+
         }
Ș
+
}
  
 
void MultHuge2(Huge A, Huge B, Huge P)
 
void MultHuge2(Huge A, Huge B, Huge P)
/* P <- A x B, varianta (lg 3) */
+
/* P <- A x B, varianta N^(lg 3) */
Î int i,j,NDig=Aî0ș>Bî0ș ? Aî0ș : Bî0ș,Needed=1,SaveA,SaveB;
+
{ int i,j,NDig=A[0]>B[0] ? A[0] : B[0],Needed=1,SaveA,SaveB;
  
 
   /* Corecteaza numarul de cifre */
 
   /* Corecteaza numarul de cifre */
 
   while (Needed<NDig) Needed<<=1;
 
   while (Needed<NDig) Needed<<=1;
   SaveA=Aî0ș;SaveB=Bî0ș;Aî0ș=Bî0ș=Needed;
+
   SaveA=A[0];SaveB=B[0];A[0]=B[0]=Needed;
   for (i=SaveA+1;i<=Needed;) Aîi++ș=0;
+
   for (i=SaveA+1;i<=Needed;) A[i++]=0;
   for (i=SaveB+1;i<=Needed;) Bîi++ș=0;
+
   for (i=SaveB+1;i<=Needed;) B[i++]=0;
 
   MultRec(A,B,P);
 
   MultRec(A,B,P);
  
 
   /* Restaureaza numarul de cifre */
 
   /* Restaureaza numarul de cifre */
   Aî0ș=SaveA;Bî0ș=SaveB;
+
   A[0]=SaveA;B[0]=SaveB;
   while (!PîPî0șș && Pî0ș>1) Pî0ș--;
+
   while (!P[P[0]] && P[0]>1) P[0]--;
Ș
+
}
 +
</syntaxhighlight>
  
 
=== Împărțirea unui vector la un scalar ===
 
=== Împărțirea unui vector la un scalar ===

Revision as of 19:27, 4 March 2018

Psihologia concursurilor de informatică

Capitolul II: Lucrul cu numere mari

De multe ori, în probleme, apar situații când este nevoie să memorăm numere întregi foarte mari (de ordinul sutelor de cifre), iar uneori trebuie să efectuăm și operații aritmetice cu aceste numere. Iată un asemenea exemplu:

ENUNȚ: Se dă un număr natural cu N\leq 1000 cifre. Se cere să se extragă rădăcina cubică a numărului. Se garantează că numărul citit este cub perfect.

Intrarea: Fișierul INPUT.TXT conține un singur rând, terminat cu EOF, pe care se află numărul, cifrele fiind nedespărțite.

Ieșirea: În fișierul OUTPUT.TXT se va tipări rădăcina cubică a numărului, pe o singură linie terminată cu EOF.

Exemplu:

INPUT.TXT OUTPUT.TXT
2097152 128

Timp de implementare: 1h 30’.

Timp de rulare: 10 secunde.

Complexitate cerută: O(N^{2}).

REZOLVARE: Problema este cât se poate de simplă din punct de vedere matematic; dificultatea apare la implementare, atât datorită structurilor de date necesare, cât mai ales datorită complexității cerute. Despre lucrul cu numere întregi (chiar și Longint) nici nu poate fi vorba, iar la lucrul cu numere reale apar erori de calcul.

Structura de date propusă pentru abordarea acestui gen de probleme este următoarea: numerele vor fi reprezentate printr-un vector de cifre zecimale. Prima poziție (poziția 0) din vector este rezervată pentru a memora numărul de cifre. Definiția C a tipului de date este:

typedef int Huge[1001];

Iată cum s-ar memora numărul „1997” într-un asemenea vector:

Psycho-fig01.png

Se observă că vectorul este oarecum „întors pe dos”. Totuși, această formă este cea mai convenabilă, pentru că ea permite implementarea cu o mai mare ușurință a operațiilor aritmetice.

Mai trebuie remarcat că pe fiecare poziție K în vector se află cifra care îl reprezintă pe 10^{{K-1}} în numărul reprezentat: în H[1] se află cifra unităților (10^{0}), în H[2] se află cifra zecilor (10^{1}), în H[3] se află cifra sutelor (10^{2}) ș.a.m.d. Formatul zecimal nu este cea mai fericită (a se citi „eficientă”) alegere. El ocupă doar patru biți din cei opt ai unui octet, deci face risipă de memorie. Dacă am alege baza de numerație 256, am folosi la maximum memoria, și în plus operațiile ar fi cu mult mai rapide (deoarece 256 este o putere a lui 2, înmulțirile și împărțirile la 256 sunt simple deplasări la stânga și la dreapta ale reprezentărilor binare). Să facem următorul calcul ca să vedem câte cifre are în baza 256 un număr care are 1000 de cifre în baza 10:

10^{{1000}}=10\cdot (10^{3})^{{333}}\approx 10\cdot (2^{{10}})^{{333}}\approx 10\cdot (2^{8})^{{335}}=10\cdot 256^{{335}}[1]

Așadar, numărul de cifre s-a redus cam de trei ori. Un algoritm liniar ar funcționa de trei ori mai repede pe reprezentări în baza 256, iar unul pătratic ar funcționa de nouă ori mai repede. Inconvenientul major este dificultatea depanării unui program care operează într-o bază aritmetică atât de străină nouă. Vom rămâne deci la baza 10, cu mențiunea că acei mai temerari dintre voi pot încerca folosirea bazei 256.

Numerele reprezentate pe vectori le vom numi pur și simplu „vectori”, iar numerele reprezentate printr-un tip ordinal de date le vom numi, printr-o analogie ușor forțată cu matematica, „scalari”. Să vedem acum cum se efectuează operațiile elementare pe aceste numere.

Inițializarea

Un vector poate fi inițializat în trei feluri: cu 0, cu un scalar sau cu un alt vector.

La inițializarea cu 0, singurul lucru pe care îl avem de făcut este să setăm numărul de cifre pe 0. De aceea, este practic inutil să implementăm această funcție ca atare; putem folosi în loc singura instrucțiune pe care ea o conține.

void Atrib0(Huge H)
{ H[0]=0; }

La inițializarea cu un scalar nenul, trebuie să așezăm fiecare cifră pe poziția corespunzătoare, aflând în paralel și numărul de cifre. Se începe cu cifra unităților, și la fiecare pas se pune în vector cifra cea mai puțin semnificativă, după care numărul de reprezentat se împarte la 10 (neglijându-se restul), iar numărul de cifre se incrementează.

void AtribValue(Huge H, unsigned long X)
{ H[0]=0;
  while (X)
    { H[++H[0]]=X%10;
      X/=10;
    }
}

Iată, de exemplu, cum se pune pe vector numărul 195:

Psycho-fig02.png

În sfârșit, inițializarea unui vector cu altul se face printr-o simplă copiere (se pot folosi cu succes rutine de lucru cu memoria, cum ar fi FillChar în Pascal sau memmove în C). Pentru eleganță, poate fi folosită și atribuirea cifră cu cifră:

void AtribHuge(Huge H, Huge X)
{ int i;

  for (i=0;i<=X[0];i++) H[i]=X[i];
}

Compararea

Pentru a compara două numere „uriașe”, începem prin a afla numărul de cifre semnificative (deoarece, în urma anumitor operații pot rezulta zerouri nesemnificative care „atârnă” totuși la numărul de cifre). Aceasta se face decrementând numărul de cifre al fiecărui număr atâta timp cât cifra cea mai semnificativă este 0. După ce ne-am asigurat asupra acestui punct, comparăm numărul de cifre al celor două cifre. Numărul cu mai multe cifre este cel mai mare. Dacă ambele numere au același număr de cifre, pornim de la cifra cea mai semnificativă și comparăm cifrele celor două numere până la prima diferență întâlnită. În acest moment, numărul a cărui cifră este mai mare este el însuși mai mare. Dacă toate cifrele numerelor sunt egale două câte două, atunci în mod evident numerele sunt egale.

După cum se vede, algoritmul seamănă foarte bine cu ceea ce s-a învățat la matematică prin clasa a II-a (doar că atunci nu ni s-a spus că este vorba de un „algoritm”). Rutina de mai jos compară două numere uriașe H_{1} și H_{2} și întoarce -1, 0 sau 1, după H_{1} este mai mic, egal sau mai mare decât H_{2}.

int Sgn(Huge H1,Huge H2)
{ int i;

  while (!H1[H1[0]] && H1[0]) H1[0]--;
  while (!H2[H2[0]] && H2[0]) H2[0]--;
  if (H1[0]!=H2[0])
    return H1[0]<H2[0] ? -1 : 1;
  i=H1[0];
  while (H1[i]==H2[i] && i) i--;
  return H1[i]<H2[i] ? -1 : H1[i]==H2[i] ? 0 : 1;
}

Adunarea a doi vectori

Fiind dați doi vectori, A cu M cifre și B cu N cifre, adunarea lor se face în mod obișnuit, ca la aritmetică. Pentru a evita testele de depășire, se recomandă să se completeze mai întâi vectorul cu mai puține cifre cu zerouri până la dimensiunea vectorului mai mare. La sfârșit, vectorul sumă va avea dimensiunea vectorului mai mare dintre A și B, sau cu 1 mai mult dacă apare transport de la cifra cea mai semnificativă. Procedura de mai jos adaugă numărul B la numărul A.

void Add(Huge A, Huge B)
/* A <- A+B */
{ int i,T=0;

  if (B[0]>A[0])
    { for (i=A[0]+1;i<=B[0];) A[i++]=0;
      A[0]=B[0];
    }
    else for (i=B[0]+1;i<=A[0];) B[i++]=0;
  for (i=1;i<=A[0];i++)
    { A[i]+=B[i]+T;
      T=A[i]/10;
      A[i]%=10;
    }
  if (T) A[++A[0]]=T;
}

Scăderea a doi vectori

Se dau doi vectori A și B și se cere să se calculeze diferența A - B. Se presupune B\leq A (acest lucru se poate testa cu funcția Sgn). Pentru aceasta, se pornește de la cifra unităților și se memorează la fiecare pas „împrumutul” care trebuie efectuat de la cifra de ordin imediat superior (împrumutul poate fi doar 0 sau 1). Deoarece B\leq A, avem garanția că pentru a scădea cifra cea mai semnificativă a lui B din cifra cea mai semnificativă a lui A nu e nevoie de împrumut.

void Subtract(Huge A, Huge B)
/* A <- A-B */
{ int i, T=0;

  for (i=B[0]+1;i<=A[0];) B[i++]=0;
  for (i=1;i<=A[0];i++)
    A[i]+= (T=(A[i]-=B[i]+T)<0) ? 10 : 0;
    /* Adica A[i]=A[i]-(B[i]+T);
       if (A[i]<0) T=1; else T=0;
       if (T) A[i]+=10; */
  while (!A[A[0]]) A[0]--;
}

Înmulțirea și împărțirea cu puteri ale lui 10

Aceste funcții sunt uneori utile. Ele pot folosi și funcțiile de înmulțire a unui vector cu un scalar, care vor fi prezentate mai jos, dar se pot face și prin deplasarea întregului număr spre stânga sau spre dreapta. De exemplu, înmulțirea unui număr cu 100 presupune deplasarea lui cu două poziții înspre cifra cea mai semnificativă și adăugarea a două zerouri la coadă. Principalul avantaj al scrierii unor funcții separate pentru înmulțirea cu 10, 100, ..., este că se pot folosi rutinele de acces direct al memoriei (FillChar, respectiv memmove). Iată funcțiile care realizează deplasarea vectorilor, atât prin mutarea blocurilor de memorie, cât și prin atribuiri succesive.

void Shl(Huge H, int Count)
/* H <- H*10^Count */
{ 
  /* Shifteaza vectorul cu Count pozitii */
  memmove(&H[Count+1],&H[1],sizeof(int)*H[0]);
  /* Umple primele Count pozitii cu 0 */
  memset(&H[1],0,sizeof(int)*Count);
  /* Incrementeaza numarul de cifre */
  H[0]+=Count;
}

void Shl2(Huge H, int Count)
/* H <- H*10^Count */
{ int i;

  /* Shifteaza vectorul cu Count pozitii */
  for (i=H[0];i;i--) H[i+Count]=H[i];
  /* Umple primele Count pozitii cu 0 */
  for (i=1;i<=Count;) H[i++]=0;
  /* Incrementeaza numarul de cifre */
  H[0]+=Count;
}

void Shr(Huge H, int Count)
/* H <- H/10^Count */
{ 
  /* Shifteaza vectorul cu Count pozitii */
  memmove(&H[1],&H[Count+1],sizeof(int)*(H[0]-Count));
  /* Decrementeaza numarul de cifre */
  H[0]-=Count;
}

void Shr2(Huge H, int Count)
/* H <- H/10^Count */
{ int i;

  /* Shifteaza vectorul cu Count pozitii */
  for (i=Count+1;i<=H[0];i++) H[i-Count]=H[i];
  /* Decrementeaza numarul de cifre */
  H[0]-=Count;
}

Înmulțirea unui vector cu un scalar

Și această operație este o simplă implementare a modului manual de efectuare a calculului. La înmulțirea „de mână” a unui număr mare cu unul de o singură cifră, noi parcurgem deînmulțitul de la sfârșit la început, și pentru fiecare cifră efectuăm următoarele operații:

  • Înmulțim cifra respectivă cu înmulțitorul;
  • Adăugăm „transportul” de la înmulțirea precedentă;
  • Separăm ultima cifră a rezultatului și o trecem la produs;
  • Celelalte cifre are rezultatului constituie transportul pentru următoarea înmulțire;
  • La sfârșitul înmulțirii, dacă există transport, acesta are o singură cifră, care se scrie înaintea rezultatului.

Exact același procedeu se poate aplica și dacă înmulțitorul are mai mult de o cifră. Singura deosebire este că transportul poate avea mai multe cifre (poate fi mai mare ca 9). Din această cauză, la sfârșitul înmulțirii, poate rămâne un transport de mai multe cifre, care se vor scrie înaintea rezultatului. Iată de exemplu cum se calculează produsul 312 x 87:

Psycho-fig03.png

Procedura este descrisă mai jos:

void Mult(Huge H, unsigned long X)
/* H <- H*X */
{ int i;
  unsigned long T=0;

  for (i=1;i<=H[0];i++)
    { H[i]=H[i]*X+T;
      T=H[i]/10;
      H[i]=H[i]%10;
    }
  while (T) /* Cat timp exista transport */
    { H[++H[0]]=T%10;
      T/=10;
    }
}

Înmulțirea a doi vectori

Dacă ambele numere au dimensiuni mari și se reprezintă pe tipul de date Huge, produsul lor se calculează înmulțind fiecare cifră a deînmulțitului cu fiecare cifră a înmulțitorului și trecând rezultatul la ordinul de mărime (exponentul lui 10) cuvenit. De fapt, același lucru îl facem și noi pe hârtie. Considerând același exemplu, în care ambele numere sunt „uriașe”, produsul lor se calculează de mână astfel:

Psycho-fig04.png

S-a luat deci fiecare cifră a înmulțitorului și s-a efectuat produsul parțial corespunzător, corectând la fiecare pas rezultatul prin calculul transportului. Rezultatul pentru fiecare produs parțial s-a scris din ce în ce mai în stânga, pentru a se alinia corect ordinele de mărime. Acest procedeu este oarecum incomod de implementat. Se pot face însă unele observații care ușurează mult scrierea codului:

  • Prin înmulțirea cifrei cu ordinul de mărime 10^{i} din primul număr cu cifra cu ordinul de mărime 10^{j} din al doilea număr se obține o cifră corespunzătoare ordinului de mărime 10^{{i+j}} în rezultat (sau se obține un număr cu mai mult de o singură cifră, caz în care transportul merge la cifra corespunzătoare ordinului de mărime 10^{{i+j+1}}).
  • Dacă numerele au M și respectiv N cifre, atunci produsul lor va avea fie M + N fie M + N - 1 cifre. Într-adevăr, dacă numărul A are M cifre, atunci 10^{{M-1}}\leq A<10^{M} și 10^{{N-1}}\leq B<10^{N}, de unde rezultă 10^{{M+N-2}}\leq A\times B<10^{{M+N}}.
  • La calculul produselor parțiale se poate omite calculul transportului, acesta urmând a se face la sfârșit. Cu alte cuvinte, într-o primă fază putem pur și simplu să înmulțim cifră cu cifră și să adunăm toate produsele de aceeași putere, obținând un număr cu „cifre” mai mari ca 9, pe care îl aducem la forma normală printr-o singură parcurgere. Să reluăm același exemplu:

Psycho-fig05.png

Această operație efectuează M\times N produse de cifre și M + N (sau M + N - 1, după caz) „transporturi” pentru aflarea rezultatului, deci are complexitatea O(M\times N). Iată și implementarea:

void MultHuge(Huge A, Huge B, Huge C)
/* C <- A x B */
{ int i,j,T=0;

  C[0]=A[0]+B[0]-1;
  for (i=1;i<=A[0]+B[0];) C[i++]=0;
  for (i=1;i<=A[0];i++)
    for (j=1;j<=B[0];j++)
      C[i+j-1]+=A[i]*B[j];
  for (i=1;i<=C[0];i++)
    { T=(C[i]+=T)/10;
      C[i]%=10;
    }
  if (T) C[++C[0]]=T;
}

Mai există o altă modalitate de a înmulți două numere de câte N cifre fiecare, care are complexitatea O(N^{{\log _{2}3}})\approx O(N^{{1,58}})\approx O(N{\sqrt  {N}}). Ea derivă de un algoritm propus de Strassen în 1969 pentru înmulțirea matricelor. Diferența se face simțită, ce-i drept pentru valori mari ale lui N, dar constanta multiplicativă crește destul de mult și, în plus, soluția e mai greu de implementat; de aceea nu recomandăm implementarea ei în timpul concursului. Ideea de bază este să se micșoreze numărul de înmulțiri și să se mărească numărul de adunări, deoarece adunarea a doi vectori se face în O(N), pe când înmulțirea se face în O(N^{2}). Să considerăm întregii A și B, fiecare de câte N cifre. Trebuie să-i înmulțim într-un timp T(N) mai bun decât O(N^{2}). Să împărțim numărul A în două „bucăți” C și D, fiecare de câte N / 2 cifre, iar întregul B în două bucăți E și F, tot de câte N / 2 cifre (presupunem că N este par):

Psycho-fig06.png

Atunci se poate scrie relația

A\times B=(C\times 10^{{N/2}}+D)\times (E\times 10^{{N/2}}+F)=CE\times 10^{{N}}+(CF+DE)\times 10^{{N/2}}+DF

Pentru a putea calcula produsul A\times B avem, prin urmare, nevoie de patru produse parțiale, de trei adunări și de două înmulțiri cu puteri ale lui 10. Adunările și înmulțirile cu puteri ale lui 10 se fac în timp liniar. Dacă efectuăm cele patru produse parțiale prin patru înmulțiri, rezultă formula recurentă de calcul

T(N)=4T(N/2)+O(N)

care duce prin eliminarea recurenței la T(N)\in O(N^{2}). Cu alte cuvinte, încă n-am câștigat nimic. Trebuie să reușim cumva să reducem numărul de înmulțiri de la 4 la 3, chiar dacă prin aceasta vom mări numărul de adunări necesare. Să definim produsul

G=(C+D)\times (E+F)=CE+CF+DE+DF=CE+DF+(CF+DE)

Atunci putem scrie:

A\times B=CE\times 10^{{N}}+(G-CE-DF)\times 10^{{N/2}}+DF

Pentru această variantă, folosim doar trei înmulțiri, și chiar dacă avem nevoie de șase adunări și scăderi și două înmulțiri cu puteri ale lui 10, complexitatea se va reduce la O(N^{{\log _{2}3}}). În cazul în care N este o putere a lui 2, împărțirea în două a numerelor se poate face fără probleme la fiecare pas, până se ajunge la numere de o singură cifră, care se înmulțesc direct. În cazul în care N nu este o putere a lui 2, este comod să se completeze numerele cu zerouri până la o putere a lui 2. În funcțiile descrise mai jos, MultRec nu face decât înmulțirea recursivă, pe când MultVect2 se ocupă și de corectarea numărului de cifre (incrementarea până la o putere a lui 2). Pentru calculul produselor C\times E și D\times F, procedura MultRec se autoapelează; pentru calcularea produsului (C+D)\times (E+F), însă, este nevoie să fie apelată procedura MultVect2, deoarece prin cele două adunări poate să apară o creștere a numărului de cifre al factorilor, care în acest caz trebuie readuși la un număr de cifre egal cu o putere a lui 2.

void MultHuge2(Huge A, Huge B, Huge P);

void MultRec(Huge A, Huge B, Huge P)
{ Huge C,D,E,F,CE,DF;

  if (A[0]==1)
    { P[1]=A[1]*B[1];
      P[0]=(P[2]=P[1]/10)>0 ? 2 : 1;
      P[1]%=10;
    }
    else { P[0]=0;
           AtribHuge(C,A);Shr(C,A[0]/2);
           AtribHuge(D,A);D[0]=A[0]/2;
           AtribHuge(E,B);Shr(E,B[0]/2);
           AtribHuge(F,B);F[0]=B[0]/2;
           MultRec(C,E,CE);MultRec(D,F,DF);
           Add(C,D);Add(E,F);
           MultHuge2(C,E,P);
           Subtract(P,CE);Subtract(P,DF);
           Shl(P,A[0]/2);
           Shl(CE,A[0]);Add(P,CE);
           Add(P,DF);
         }
}

void MultHuge2(Huge A, Huge B, Huge P)
/* P <- A x B, varianta N^(lg 3) */
{ int i,j,NDig=A[0]>B[0] ? A[0] : B[0],Needed=1,SaveA,SaveB;

  /* Corecteaza numarul de cifre */
  while (Needed<NDig) Needed<<=1;
  SaveA=A[0];SaveB=B[0];A[0]=B[0]=Needed;
  for (i=SaveA+1;i<=Needed;) A[i++]=0;
  for (i=SaveB+1;i<=Needed;) B[i++]=0;
  MultRec(A,B,P);

  /* Restaureaza numarul de cifre */
  A[0]=SaveA;B[0]=SaveB;
  while (!P[P[0]] && P[0]>1) P[0]--;
}

Împărțirea unui vector la un scalar

Ne propunem să scriem o funcție care să împartă numărul A de tip Huge la scalarul B, să rețină valoarea câtului tot în numărul A și să întoarcă restul (care este o variabilă scalară). Să pornim de la un exemplu particular și să generalizăm apoi procedeul: să calculăm câtul și restul împărțirii lui 1997 la 7. Cu alte cuvinte, să găsim acele numere C de tip Huge și RÎÎ0, 1, 2, 3, 4, 5, 6Ș cu proprietatea că 1997 = 7xC + R.



La fiecare pas se coboară câte o cifră de la deîmpărțit alături de numărul deja existent (care inițial este 0), apoi rezultatul se împarte la împărțitor (7 în cazul nostru). Câtul este întotdeauna o cifră și se va depune la sfârșitul câtului împărțirii, iar restul va fi folosit pentru următoarea împărțire. Restul care rămâne după ultima împărțire este tocmai R pe care îl căutăm. Procedeul funcționează și atunci când deîmpărțitul are mai multe cifre. La sfârșit trebuie să decrementăm corespunzător numărul de cifre al câtului, prin neglijarea zerourilor create la începutul numărului. Numărul maxim de cifre al câtului este egal cu cel al deîmpărțitului.

unsigned long Divide(Huge A, unsigned long X) /* A <- A/X si intoarce A%X */ Î int i;

 unsigned long R=0;
 for (i=Aî0ș;i;i--)
   Î Aîiș=(R=10*R+Aîiș)/X;
     R%=X;
   Ș
 while (!AîAî0șș && Aî0ș>1) Aî0ș--;
 return R;

Ș

Dacă dorim numai să aflăm restul împărțirii, nu mai avem nevoie decât să recalculăm restul la fiecare pas, fără a mai modifica vectorul A:


unsigned long Mod(Huge A, unsigned long X) /* Intoarce A%X */ Î int i;

 unsigned long R=0;
 for (i=Aî0ș;i;i--)
   R=(10*R+Aîiș)%X;
 return R;

Ș


Împărțirea a doi vectori

Dacă se dau doi vectori A și B și se cere să se afle câtul C și restul R, etapele de parcurs sunt aceleași ca la punctul precedent, numai că operatorii „/” și „%” trebuie implementați de utilizator, ei nefiind definiți pentru vectori. Cu alte cuvinte, după ce „coborâm” la fiecare pas următoarea cifră de la deîmpărțit, trebuie să aflăm cea mai mare cifră X astfel încât împărțitorul să se cuprindă de X ori în restul de la momentul respectiv. Acest lucru se face cel mai comod prin adunări repetate: pornim cu cifra X=0 și o incrementăm, micșorând concomitent restul, până când restul care rămâne este prea mic. Să efectuăm aceeași împărțire, 1997:7, considerând că ambele numere sunt reprezentate pe tipul Huge.



Cazul cel mai defavorabil (când X=9) presupune 9 scăderi și 10 comparații, cazul cel mai favorabil (când X=0) presupune numai o comparație, deci cazul mediu presupune 4 scăderi și 5 comparații. Căutarea lui X se poate face și binar, prin înjumătățirea intervalului, ceea ce reduce timpul mediu de căutare la aproximativ 3 comparații și trei înmulțiri, dar codul se complică nejustificat de mult (de cele mai multe ori).

void DivideHuge(Huge A, Huge B, Huge C, Huge R) /* A/B = C rest R */ Î int i;

 Rî0ș=0;Cî0ș=Aî0ș;
 for (i=Aî0ș;i;i--)
   Î Shl(R,1);Rî1ș=Aîiș;
     Cîiș=0;
     while (Sgn(B,R)!=1)
       Î Cîiș++;
         Subtract(R,B);
       Ș
   Ș
 while (!CîCî0șș && Cî0ș>1) Cî0ș--;

Ș


Extragerea rădăcinii cubice

Vom sări peste prezentarea algoritmului de extragere a rădăcinii pătrate, pe care îl vom lăsa ca temă cititorului, și ne vom îndrepta atenția asupra celui de extragere a rădăcinii cubice, care este puțin mai complicat, dar care poate fi ușor extins pentru rădăcini de orice ordin. Problema este exact cea din enunț, așa că vom porni de la exemplul dat. Să notăm și . Cum se află X ? O primă variantă ar fi căutarea binară a rădăcinii, prin înjumătățirea intervalului. Inițial se pornește cu intervalul (1,A), deoarece rădăcina cubică se află undeva între 1 și A (evident, încadrarea este mai mult decât acoperitoare; ea ar putea fi mai limitativă, dar nu ar reduce timpul de lucru decât cu câteva iterații). La fiecare pas, intervalul va fi înjumătățit. Cum, probabil că știți deja; se ia jumătatea intervalului, se ridică la puterea a treia și se compară cu A. Dacă este mai mare, înseamnă că rădăcina trebuie căutată în jumătatea inferioară a intervalului. Dacă este mai mică, vom continua căutarea în jumătatea superioară a intervalului. Dacă cele două numere sunt egale, înseamnă că am găsit tocmai ce ne interesa. Prima variantă a pseudocodului este:

1 citeste A cu N cifre 2 Lo ¬ 1, Hi ¬ A, X ¬ 0 3 cat timp X=0 4 Mid ¬ (Lo+Hi)/2 5 daca Mid3<A 6 atunci Lo ¬ Mid+1 7 altfel daca Mid3>A 8 atunci Hi ¬ Mid-1 9 altfel X ¬ Mid

În cazul cel mai rău, algoritmul de mai sus efectuează log2 A înjumătățiri de interval, fiecare din ele presupunând o adunare, o împărțire la 2 și o ridicare la cub. Dintre aceste operații, cea mai costisitoare este ridicarea la cub, O( N2 ). Complexitatea totală este prin urmare O( N2 log2 A ). Deoarece A are ordinul de mărime 10N, rezultă complexitatea O( N3 log 10 )= O( N3 ), adică mai proastă decât cea cerută (de altfel, un algoritm cu această complexitate nici nu s-ar încadra în timp pentru N=1000). Dacă timpul ne permite, trebuie să căutăm altă metodă.

În exemplul ales, să observăm că 106£A<109, de unde deducem că 102£X<103. Cu alte cuvinte, X are 3 cifre. În cazul general, dacă A are N cifre, atunci X are éN/3ù cifre (prin éN/3ù se înțelege „cel mai mic întreg mai mare sau egal cu N/3”). Care ar putea fi prima cifră a lui X ? Dacă X începe cu cifra 2, atunci 200£X<300 Þ 8.000.000£2.097.152<27.000.000, ceea ce este fals. Cu atât mai puțin poate prima cifră a lui X să fie mai mare ca 2. Rezultă că prima cifră a lui X este 1. De altfel, pentru a afla acest lucru, putem să și neglijăm ultimele 6 cifre ale lui A. Ne interesează doar prima cifră, cea a milioanelor, iar prima cifră a lui X o alegem în așa fel încât cubul ei să fie mai mic sau egal cu 2.

Ce putem spune despre a doua cifră ? Dacă ar fi 3, atunci 130£X<140 Þ 2.197.000£2.097.152<2.744.000, fals (deci cifra este cel mult 2). Dacă ar fi 1, atunci 110£X<120 Þ 1.331.000£2.097.152<1.728.000, fals. Rezultă că a doua cifră este obligatoriu 2. Analog, putem neglija ultimele trei cifre ale lui A, iar a doua cifră a lui X este cel mai mare C pentru care . Pentru a afla ultima cifră, aplicăm același raționament: Dacă ar fi 9, atunci X=129 și ar rezulta 2146688=X3=2097152, absurd. Dacă considerăm că cifra este 8, atunci calculul se verifică. Am aflat așadar că X=128.

Procedeul general este următorul: dându-se un număr A cu N cifre, îl completăm cu zerouri nesemnificative până când N se divide cu 3 (poate fi necesar să adăugăm maxim două zerouri). Numărul de cifre semnificative ale rădăcinii cubice este M=N/3. Aflăm pe rând fiecare cifră, începând cu cea mai semnificativă. Să presupunem că am aflat cifrele XM, XM-1,..., XK+1. Cifra XK este cea mai mare cifră pentru care numărul . Cifra XK este unică, deoarece există, în general, mai multe cifre care verifică proprietatea cerută, dar una singură este „cea mai mare”. O a doua versiune a pseudocodului este deci:

1 citeste A cu N cifre 2 X ¬ 0; T ¬ 0 3 cat timp N%3<>0 adauga un 0 nesemnificativ 4 repeta de N/3 ori 5 adauga la T urmatoarele 3 cifre din A 6 adauga la X cea mai mare cifra astfel incat X3 £ T

Să evaluăm complexitatea acestei versiuni. Linia 1 se execută în timp liniar, O(N). Liniile 2 și 3 se execută în timp constant. Linia 5 se execută în O(N), iar linia 6 presupune adăugarea unei cifre (O(N)) și o ridicare la cub, adică două înmulțiri (O(N2)). Deoarece liniile 5 și 6 se execută de N/3 ori (linia 4), rezultă o complexitate de O(N3). Iată că nici acest algoritm nu a adus îmbunătățiri și pare și ceva mai greu de implementat. El poate fi totuși modificat pentru a-i scădea complexitatea la O(N2).

Principalul neajuns al său este efectuarea ridicării la cub, care se face în O(N2). Dacă am putea să-l aflăm la fiecare pas pe X3 fără a efectua înmulțiri, adică în timp liniar, atunci întregul algoritm ar avea complexitate pătratică. Bineînțeles, prima întrebare care vine pe buzele cititorului este „cum să ridicăm la cub fără să facem înmulțiri?”. Să nu uităm însă ceva: că noi nu-l cunoaștem numai pe X. Îl cunoaștem și pe X de la pasul anterior, care avea o cifră mai puțin (îl vom boteza OldX). Să presupunem că, printr-o metodă oarecare, am reușit să-l ridicăm pe OldX la puterile a doua și a treia (și am obținut numerele OldX2 și OldX3). Cum putem, cunoscând aceste trei numere, precum și noua cifră ce se va adăuga la sfârșitul lui X (să-i spunem C ), să-l aflăm pe X, pătratul și cubul său? Nu e prea greu:


Iată așadar că pentru a afla noile valori ale puterilor 1, 2 și 3 ale lui X, folosindu-le pe cele vechi, nu avem nevoie decât de adunări și de înmulțiri cu numere mici (de ordinul miilor). Toate aceste operații se fac în timp liniar, deci am reușit să găsim un algoritm pătratic. Iată mai jos sursa C: void FindDigit(Huge L,Huge NewL2,Huge NewL3,

    Huge OldL,Huge OldL2,Huge OldL3,Huge Target)

Î Huge Aux;

 Lî1ș=10;
 do
   Î Lî1ș--;
     /* Trebuie calculat LÂ3. Se stiu OldL (L/10)
        si noua cifra Lî1ș. Deci (OldL*10+Lî1ș)Â3=?
        Se aplica binomul lui Newton. */
     AtribHuge(NewL3,OldL3);Shl(NewL3,3);
     AtribHuge(Aux,OldL2);Mult(Aux,300*Lî1ș);
     Add(NewL3,Aux);
     AtribHuge(Aux,OldL);Mult(Aux,30*Lî1ș*Lî1ș);
     Add(NewL3,Aux);
     AtribValue(Aux,Lî1ș*Lî1ș*Lî1ș);
     Add(NewL3,Aux);
   Ș
 while (Sgn(NewL3,Target)==1);
 /* Aceeasi operatie pentru LÂ2 */
 AtribHuge(NewL2,OldL2);Shl(NewL2,2);
 AtribHuge(Aux,OldL);Mult(Aux,20*Lî1ș);
 Add(NewL2,Aux);
 AtribValue(Aux,Lî1ș*Lî1ș);
 Add(NewL2,Aux);
 /* Noile valori devin 'vechi' */
 AtribHuge(OldL2,NewL2);
 AtribHuge(OldL,L);
 AtribHuge(OldL3,NewL3);

Ș

void CubeRoot(Huge A, Huge X) Î Huge Target,OldX,OldX2,OldX3,NewX2,NewX3;

 int i;
 /* Se initializeaza vectorii cu 0 (nici o cifra) */
 OldXî0ș=OldX2î0ș=OldX3î0ș=Xî0ș=0;
 for (i=1;i<=(Aî0ș+2)/3;i++)
   Î AtribHuge(Target,A);
     Shr(Target,3*((Aî0ș+2)/3-i));
     Shl(X,1);
     FindDigit(X,NewX2,NewX3,OldX,OldX2,OldX3,Target);
   Ș

Ș

Acum nu mai avem decât să scriem rutinele de intrare/ieșire și programul principal:

  1. include <stdio.h>
  2. include <mem.h>
  3. define NMax 1000

typedef int HugeîNMax+3ș; Huge A,X; /* Aî0ș si Xî0ș indica numarul de cifre */

void ReadData(void) Î FILE *F=fopen("input.txt","rt");

 int C,i;
 Aî0ș=0;
 do Aî++Aî0șș=(C=fgetc(F))-'0';
 while (C!=EOF);
 Aî0ș--;
 fclose(F);
 /* Intoarce vectorul pe dos */
 for (i=1;i<=Aî0ș/2;i++)
   Aîiș=(AîișÂAîAî0ș+1-iș)Â(AîAî0ș+1-iș=Aîiș);

Ș

void WriteSolution(void) Î FILE *F=fopen("output.txt","wt");

 int i=Xî0ș;
 while (!Xîiș) i--;
 while (i) fputc(Xîi--ș+'0',F);
 fclose(F);

Ș

void main(void) Î

 ReadData();
 CubeRoot(A,X);
 WriteSolution();

Ș

Pentru a extinde această metodă la rădăcini de orice ordin K, trebuie numai să ținem cont de expresia binomului lui Newton:


Presupunând că avem calculate toate puterile de la 1 la p ale lui OldX, se poate calcula noua valoare a lui Xp folosind numai adunări și înmulțiri cu scalari. În felul acesta se pot calcula în timp liniar valorile lui X, X2, X3, ...., XK.

  1. Greșeală în original