Exponențierea logaritmică în C++. Generarea rapidă a șirurilor recurente

de Iulian Oleniuc | 07/07/2018 | Algoritmică și structuri de date, Matematică

Așa cum îi spune și numele, exponențierea logaritmică (exponențierea rapidă, ridicarea la putere în timp logaritmic, exponențierea prin ridicare la pătrat) este metoda optimă pentru a ridica un număr $x$ la o putere $n$ ce are complexitatea logaritmică. În acest articol voi prezenta cum se implementează exponențierea logaritmică în C++, și cum se poate utiliza acest concept pentru a genera eficient al $n$ lea termen dintr-un șir de numere ce se bazează pe o recurență liniară.

Exponențierea liniară

Soluția naivă pentru a calcula $x^n$ ce se bazează pe definiția puterii, este să inițializăm o variabilă $p$ cu $1$ iar apoi să o înmulțim cu $x$ de $n$ ori. Iată implementarea iterativă a algoritmului:

int pwr(int x, int n) {
    int p = 1;
    for (int i = 0; i < n; i++)
        p *= x;
    return p;
}

Mai jos este și o implementare recursivă ce se bazează pe recurența următoare:

x^n = \begin{cases} x \cdot x^{n - 1} & \text{pentru } n \gt 0\\ 1 & \text{pentru } n = 0 \end{cases}

int pwr(int x, int n) {
    if (!n)
        return 1;
    return x * pwr(x, n - 1);
}

Se observă că, indiferent de modul de implementare, această metodă are complexitatea liniară, mai exact $O(n)$ Bine, în cazul variantei recursive se face un apel în plus, $\mathrm{pwr}(x, 0)$ dar nu contează.

Exponențierea logaritmică

Pornim de la un alt mod de a defini ridicarea la putere:

x^n = \begin{cases} x (x^2)^{(n - 1) / 2} & \text{pentru } n \text{ impar}\\ (x^2)^{n / 2} & \text{pentru } n \text{ par} \end{cases}

De data asta, voi începe cu varianta recursivă, pentru că e mai ușor de înțeles.

Implementare recursivă

int pwr(int x, int n) {
    if (!n)
        return 1;
    if (n % 2)
        return x * pwr(x * x, n / 2);
    return pwr(x * x, n / 2);
}

Știm deja că acest algoritm are complexitatea $O(\log n)$ dar hai să o și demonstrăm. De exemplu, pentru apelul $\mathrm{pwr}(3, 618)$ lanțul de apeluri recursive arată cam așa:

\begin{align*} \mathrm{pwr}(3, 618) &= \mathrm{pwr}(3^2, 309) = 3^2 \cdot \mathrm{pwr}(3^4, 154) = 3^2 \cdot \mathrm{pwr}(3^8, 77) = 3^{10} \cdot \mathrm{pwr}(3^{16}, 38)\\ &= 3^{10} \cdot \mathrm{pwr}(3^{32}, 19) = 3^{42} \cdot \mathrm{pwr}(3^{64}, 9) = 3^{106} \cdot \mathrm{pwr}(3^{128}, 4) = 3^{106} \cdot \mathrm{pwr}(3^{256}, 2)\\ &= 3^{106} \cdot \mathrm{pwr}(3^{512}, 1) = 3^{618} \cdot \mathrm{pwr}(3^{1024}, 0) = 3^{618} \end{align*}

Se observă că la fiecare pas al doilea parametru se înjumătățește, exact cum în căutarea binară la fiecare pas se înjumătățește intervalul de căutare. Numărul de pași făcuți de exponențierea logaritmică este $\left [ \log_2 n \right ] + 1$ Acel $+ 1$ vine de la apelul final, când se returnează $1$ Iată deci, cum pentru $618$ se efectuează doar $11$ apeluri, pe când prin ridicarea la putere obișnuită, s-ar fi executat $618$ pași.

Ei bine, de obicei, pentru exponenți suficient de mari ca să se simtă diferența dintre cele două variante de ridicare la putere, rezultatul final ar depăși long long int-ul. De aceea, de cele mai multe ori, se cere doar restul lui modulo un anumit număr prim. Aplicând regulile clasice de aritmetică modulară, o funcție ce calculează acest rest arată așa:

int pwr(int x, int n) {
    if (!n)
        return 1;
    if (n % 2)
        return 1LL * x * pwr(1LL * x * x % MOD, n / 2) % MOD;
    return pwr(1LL * x * x % MOD, n / 2);
}

Implementare iterativă

Exponențierea logaritmică poate fi implementată și iterativ, dar aici este nevoie de un pic mai multă atenție.

int pwr(int x, int n) {
    int p = 1;
    while (n > 0) {
        if (n % 2) {
            p *= x;
            n--;
        }
        x *= x;
        n /= 2;
    }
    return p;
}

În variabila $p$ rețin rezultatul. Parametrul $x$ se ridică la pătrat la fiecare intrare în while, iar exponentul $n$ se înjumătățește. Când $n$ este impar, se actualizează rezultatul printr-o înmulțire cu $x$ și o decrementare a lui $n$ Întotdeauna $n$ va fi impar măcar o dată, pentru că la ultima intrare în while, $n$ va fi $1$

Generarea rapidă a șirurilor recurente prin exponențiere logaritmică

Cele mai importante două aplicații ale exponențierii rapide sunt calcularea inversului modular (despre care voi vorbi în alt articol) și generarea celui de-al $n$ lea termen dintr-un șir ce se bazează pe o recurență liniară, ambele, desigur, în timp logaritmic.

Voi lua ca exemplu Șirul lui Fibonacci, care este definit astfel:

f(n) = \begin{cases} f(n - 1) + f(n - 2) & \text{pentru } n \gt 1\\ n & \text{pentru } 0 \le n \le 1 \end{cases}

Pentru a genera al $n$ lea termen din Șirul lui Fibonacci (tot modulo un număr dat, pentru că numerele lui Fibonacci cresc foarte repede), ar trebui, în mod obișnuit, să generăm mai întâi toți termenii mai mici decât $n$ Complexitatea este $O(n)$ Există totuși o soluție mult mai bună, care necesită niște matematică de a 11-a, ce se bazează pe înmulțirea matricelor. Voi da mai întâi soluția și o explic după. Am notat al $n$ lea termen din Șirul lui Fibonacci cu $f(n)$

\begin{pmatrix} 0 & 1\\ 1 & 1 \end{pmatrix}^n = \begin{pmatrix} f(n - 1) & f(n)\\ f(n) & f(n + 1) \end{pmatrix}

Folosim inducție matematică. Presupunând că ceea ce am scris mai sus este adevărat, arătăm că:

\begin{pmatrix} 0 & 1\\ 1 & 1 \end{pmatrix} \cdot \begin{pmatrix} f(n - 1) & f(n)\\ f(n) & f(n + 1) \end{pmatrix} = \begin{pmatrix} f(n) & f(n + 1)\\ f(n + 1) & f(n + 2) \end{pmatrix}

Păi, respectând regulile de înmulțire a matricelor, obținem următoarele relații:

\begin{align*} 0 \cdot f(n - 1) + 1 \cdot f(n) &= f(n)\\ 0 \cdot f(n) + 1 \cdot f(n + 1) &= f(n + 1)\\ 1 \cdot f(n - 1) + 1 \cdot f(n) &= f(n + 1)\\ 1 \cdot f(n) + 1 \cdot f(n + 1) &= f(n + 2) \end{align*}

Care, se vede clar că sunt adevărate. Ei bine, știind că înmulțirea matricelor este asociativă, putem folosi exponențierea logaritmică pentru ridicarea matricei de mai sus la puterea $n$ Iată o sursă elegantă, de 100 de puncte, pentru problema Kfib de pe InfoArena (care cere generarea celui de-al $k$ lea termen Fibonacci):

Problema Kfib
#include <fstream>
using namespace std;

const int MOD = 666013;

ifstream fin("kfib.in");
ofstream fout("kfib.out");

struct Mat {
    int mat[2][2];
};

const Mat nullMat = {
    {{1, 0},
     {0, 1}}
};

const Mat initMat = {
    {{0, 1},
     {1, 1}}
};

Mat prod(Mat a, Mat b) {
    Mat ret;
    ret.mat[0][0] = (1LL * a.mat[0][0] * b.mat[0][0] + 1LL * a.mat[0][1] * b.mat[1][0]) % MOD;
    ret.mat[0][1] = (1LL * a.mat[0][0] * b.mat[0][1] + 1LL * a.mat[0][1] * b.mat[1][1]) % MOD;
    ret.mat[1][0] = (1LL * a.mat[1][0] * b.mat[0][0] + 1LL * a.mat[1][1] * b.mat[1][0]) % MOD;
    ret.mat[1][1] = (1LL * a.mat[1][0] * b.mat[0][1] + 1LL * a.mat[1][1] * b.mat[1][1]) % MOD;
    return ret;
}

Mat pwr(Mat mat, int n) {
    if (!n)
        return nullMat;
    if (n % 2)
        return prod(mat, pwr(prod(mat, mat), n / 2));
    return pwr(prod(mat, mat), n / 2);
}

int main() {
    int k; fin >> k;
    fout << pwr(initMat, k).mat[0][1] << '\n';
    return 0;
}

Constanta din spatele complexității logaritmice din această problemă este dată de numărul de înmulțiri pe scalari ce se efectuează la înmulțirea a două matrice. Care, pentru matricele pătratice este dimensiunea lor la puterea a treia; în cazul nostru $8$

Generalizare

Mai rămâne să găsim generalizarea acestei metode pentru generarea celui de-al $n$ lea termen al șirului recurent:

f(n) = \begin{cases} a_1 \cdot f(n - 1) + a_2 \cdot f(n - 2) + \cdots + a_k \cdot f(n - k) & \text{pentru } n \gt k\\ x_1 & \text{pentru } n = 1\\ x_2 & \text{pentru } n = 2\\ \text{ }\vdots\\ x_k & \text{pentru } n = k \end{cases}

Fără alte introduceri, relația este:

\begin{pmatrix} 0 & 1 & 0 & 0 & \cdots & 0 & 0\\ 0 & 0 & 1 & 0 & \cdots & 0 & 0\\ 0 & 0 & 0 & 1 & \cdots & 0 & 0\\ 0 & 0 & 0 & 0 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & 0 & 0 & \cdots & 0 & 1\\ a_k & a_{k - 1} & a_{k - 2} & a_{k - 3} & \cdots & a_2 & a_1 \end{pmatrix}^{n - k} \cdot \begin{pmatrix} x_1\\ x_2\\ x_3\\ x_4\\ \text{ }\vdots\\ x_{k - 1}\\ x_k \end{pmatrix} = \begin{pmatrix} x_{n - k + 1}\\ x_{n - k + 2}\\ x_{n - k + 3}\\ x_{n - k + 4}\\ \text{ }\vdots\\ x_{n - 1}\\ x_n \end{pmatrix}

Luați-vă un exemplu mic și veți vedea că este chiar o relație intuitivă; un antrenament bun este problema Iepuri de pe InfoArena. Recurența asta se poate aplica la fel de bine și la Șirul lui Fibonacci. Însă acolo, termenii șirului sunt de așa natură încât merge să lucrăm doar cu matrice. Mai exact, prima coloană din matricea de tranziție este aceeași cu vectorul de valori inițiale.

Probleme recomandate

Kfib clasică
Iepuri clasică
2șah observație faină
Cameleoni dinamică pe biți
Ikebana dinamică tractor
Recurența2 o recurență mai complexă
Power un șmen foarte util în problemele cu multe query-uri
Once again… poți calcula în timp logaritmic și alte legi de compoziție asociative

Înainte de problemele astea, un exercițiu bun ar fi să modificați formula pe matrice pentru a genera șiruri a căror recurență conține și un termen liber.

Sper că am explicat suficient de clar cum stă treaba cu exponențierea logaritmică și cu șirurile recurente. Dacă aveți vreo întrebare despre acest subiect, nu ezitați să o adresați printr-un comentariu mai jos :smile: