ASD2 - 02 - Algoritmi Divide et Impera

In questa lezione abbiamo ripassato il paradigma divide et impera, vedendo un modo più efficiente di quello tradizionale per moltiplicare due numeri tra loro. Abbiamo poi ripassato il teorema master per analizzare la complessità degli algoritmi di tipo divide et impera, e infine abbiamo concluso la lezione discutendo di alcuni problemi dell'artimetica modulare, come il calcolo del massimo comun divisore (MCD).

Se vogliamo moltiplicare due numeri naturali \(n_1\) e \(n_2\) possiamo usare l'algoritmo tradizionale imparato nelle scuole elementari. Una breve analisi ci mostra che tale algoritmo ha un costo di \(\mathcal{O}(\beta^2)\), dove \(\beta\) indica il numero di bit utilizzato per rappresentare il numero più grande tra \(n_1\) e \(n_2\).

Ci possiamo quindi fare la seguente domanda: esiste forse un modo più ingegnioso che ci permette di moltiplicare due interi eseguendo un numero di operazioni che è un \(o(\beta^2)\)?

Osservazione: È importante osservare che questo studio della complessità per l'algoritmo della moltiplicazione è utile solo nel caso in cui i due numeri che vogliamo moltiplicare \(n_1\) e \(n_2\) possono essere arbitrariamente grandi. I moderni processi processori con architettura a 64 bit infatti impiegano un tempo costante per la moltiplicazione di numeri rappresentabili da al più 64 bit.

Alla domanda precedentemente posta possiamo, in modo quasi sorprendente, dare una risposta affermativa. Un primo passo da eseguire per poterci avvicinare a risolvere il problema della moltiplicazione in modo più efficiente consiste nel cambiare la tecnica con cui risolviamo il problema: al posto di procedere in modo iterativo, possiamo utilizzare un approccio di tipo Divide et Impera.

Gli algoritmi di tipo Divide et Impera cercano di spezzare il problema in sotto-problemi, in modo da risolvere ricorsivamente i sotto-problemi, e infine ricombinare le soluzioni dei sotto-problemi per risolvere il problema iniziale.

Nel nostro caso, siano \(x\) e \(y\) i numeri che vogliamo moltiplicare e sia \(n\) il numero di bit richiesto per memorizzare il numero più grande tra \(x\) e \(y\). Andiamo a scrivere \(x\) e \(y\) in base 2 (ogni altra base \(\geq 1\) sarebbe equivalente rispetto alla complessità computazionale), per ottenere

\[\begin{split} x &= 2^{n-1} x_{n-1} + 2^{n-2} x_{n-2} + ... + 2^{n/2} x_{n/2}... + 2x_{1} + x_{0} \\ y &= 2^{n-1} y_{n-1} + 2^{n-2} y_{n-2} + ... + 2^{n/2} y_{n/2}... + 2y_{1} + y_{0} \\ \end{split}\]

con \(x_i \in \{0, 1\}\) e \(y_i \in \{0, 1\}\).

Per impostare un approccio divide et impera l'idea è quella di spezzare la rappresentazione in binario di \(x\) e \(y\) in due numeri. Facciamo vedere il procedimento per \(x\), sapendo che per \(y\) si procede in modo analogo

\[\begin{split} x_l &= 2^{n/2 - 1}x_{n-1} + ... 2^1 \cdot x_{n/2 + 1}+ x_{n/2} \\ x_r &= 2^{n/2 - 1}x_{n/2 - 1} + ... + 2^1 \cdot x_1 + x_{0} \\ \end{split}\]

Utilizzando \(x_r\) e \(x_l\), possiamo facilmente calcolare \(x\) nel seguente modo

\[\begin{split} 2^{n/2} x_l + x_r &= 2^{n/2}[ 2^{n/2 - 1}x_{n-1} + ... + 2x_{n/2 + 1} + x_{n/2}] + 2^{n/2 - 1}x_{n/2} + ... + 2x_1 + x_0 \\ &= 2^{n-1} x_{n-1} + 2^{n-2} x_{n-2} + ... + 2^{n/2} x_{n/2}... + 2x_{1} + x_{0} \\ &= x \end{split}\]

Ora, utilizzando questa rappresentazione possiamo scrivere la moltiplicazione tra \(x\) e \(y\) nel seguente modo:

\[\begin{split} x \cdot y & = (2^{\frac{n}{2}} \cdot x_l + x_r) \cdot (2^{\frac{n}{2}} \cdot y_l + y_r) \\ & = 2^n \cdot (x_l \cdot y_l) + 2^{\frac{n}{2}} \cdot (x_l \cdot y_r + x_r \cdot y_l) + (x_r \cdot y_r) \\ \end{split}\]

In questa forma l'applicazione della tecnica Divide et Impera è chiara: si spezza il problema della moltiplicazione in quattro moltiplicazioni di numeri più piccoli (da \(\frac{n}{2}\) bit ciascuno), si calcolano i quattro prodotti in question e si combinano in tempo lineare tramite una addizione per ottenere il prodotto ricercato. Notiamo poi che le moltiplicazioni che hanno come uno degli operandi la base in cui si sta lavorando, nel nostro caso \(2\), possono essere eseguite in tempo costante in quanto richiedono semplici shift delle cifre del numero. Dunque la maggior parte del costo è collocato nel calcolo delle quattro moltiplicazioni.

Chiediamoci ora se questo nuovo approccio ci ha fatto effettivametne guadagnare qualcosa da un punto di vista computazionale. Andando a rappresentare il numero di operazioni \(T(n)\) come il numero di operazioni eseguite dall'algoritmo su numeri a n-bit, abbiamo che

\[T(n) = 4 \cdot T\Big(\frac{n}{2}\Big) + \theta(n)\]

Risolvendo questa ricorrenza tramite il teorema master, che andremo a ricordare nel corso di questa lezione, troviamo che \(T(n) = \theta(n^2)\). Questo ci fa capire che per adesso non siamo ancora riusciti a ottenere dei miglioramenti in termine della complessità dell'algoritmo.

Il nostro problema è che ci serve un modo per ottenere i quattro prodotti: \(x_l \cdot y_l\), \(x_l \cdot x_r\), \(x_r \cdot y_l\), \(x_r \cdot y_r\), e l'unico modo per farlo sembra quello di fare esplicitamente quattro moltiplicazioni. La proprietà distributiva ci dice che

\[\begin{split} &(x_l + x_r) (y_l + y_r) = x_l y_l + (x_l y_r + x_r y_l) + x_r y_r \\ \iff \,\,\, &(x_l y_r + x_r y_l) = (x_l + x_r) (y_l + y_r) - x_l y_l - x_r y_r \end{split}\]

Per ridurre il numero di moltiplicazioni l'idea è quella di calcolare i tre prodotti \(x_r y_r\), \(x_l y_l\) e \((x_l + x_r) (y_l + y_r)\), e, una volta calcolati, ottenere l'ultima moltiplicazione mancante, ovvero \((x_l x_r + x_r y_l)\), come combinazione di somma e sottrazione dei prodotti già calcolati. Così facendo siamo in grado di eseguire solo tre moltiplcazioni, al posto delle solite quattro, trasformando quindi la ricorrenza in \(T(n) = 3 \cdot T(\frac{n}{2}) + \theta(n)\), che, se risolta, ci da una complessità di \(T(n) = \theta(n^{log_2 3})\), che è un \(o(n^2)\)!

L'osservazione che permette di ridurre 4 moltiplicazioni a 3 moltiplicazioni e qualche addizione e sottrazione fu fatta da Karatsuba nel 1960, ed è per questo che l'algoritmo che implementa questa idea viene chiamato algoritmo di Karatsuba. Andiamo quindi a proporre lo pseudocodice di tale algoritmo

Andiamo adesso a vedere un teorema fondamentale per risolvere le ricorrenze degli algoritmi di tipo Divide et Impera.

Teorema (Master Theorem): Siano \(a > 0\), \(b > 1\), \(d \geq 0\) costanti, e sia \(T\) una ricorrenza della forma

\[T(n) = a \cdot T\Big(\frac{n}{b}\Big) + \mathcal{O}(n^d)\]

allora, le soluzioni della ricorrenza sono

\[T(n) = \begin{cases} \mathcal{O}(n^d) & \,\, ,\,\, d > \log_b a \\ \mathcal{O}(n^d \log n) & \,\, ,\,\, d = \log_b a \\ \mathcal{O}(n^{\log_b a}) & \,\, ,\,\, d < \log_b a \\ \end{cases}\]

Per la dimostrazione utilizzeremo due fatti importanti, che sono:

  • La somma \(\sum_{i = 0}^k {x^i}\) può essere calcolata definendo \(S := \sum_{i = 0}^k {x^i}\) e notando che

    \[\begin{split} S - x \cdot S &= (1 + x^1 + x^2 + ... + x^k) - x \cdot (1 + x^1 + x^2 + ... + x^k) \\ &= (1 + x^1 + x^2 + ... + x^k) - (x + x^2 + x^3 + ... + x^k + x^{k+1}) \\ &= 1 - x^{k+1} \end{split}\]

    quindi, se \(x \neq 1\), abbiamo che

    \[(1 - x) \cdot S = 1 - x^{k+1} \iff S = \frac{1 - x^{k+1}}{1 - x} = \frac{x^{k+1} - 1}{x - 1}\]

    in generale vale il seguente risultato

    \[\sum_{i = 0}^k {x^i} = \begin{cases} \frac{x^{k+1} - 1}{x - 1} & \,\,,x \neq 1 \\ k & \,\,,x = 1 \\ \end{cases}\]

  • Utilizzando le proprietà del logaritmo è possibile dimostrare che \(a^{\log_b n} = n^{\log_b a}\). Notiamo infatti che

    \[\begin{split} a^{\log_b n} &= b^{\log_b(a^{\log_b n})} \\ &= b^{\log_b n \cdot \log_b a} \\ &= (b^{\log_b n})^{\log_b a} \\ &= n^{\log_b a} \end{split}\]

Procediamo quindi con la dimostrazione.

Dimostrazione: Iniziamo esplicitando la formula di \(T(n)\)

\[\begin{split} T(n) &= a \cdot T\Big(\frac{n}{b}\Big) + \mathcal{O}(n^d) \\ &\leq a \cdot T\Big(\frac{n}{b}\Big) + c \cdot n^d \\ &\leq a^2 \cdot T\Big(\frac{n}{b^2}\Big) + c \cdot n^d \cdot \Big(1 + \frac{a}{b^d}\Big) \\ &\leq \ldots \\ &\leq a^{\log_b n} \cdot T(1) + c \cdot n^d \cdot \Big( \sum_{i = 0}^{\log_b n - 1} \big(\frac{a}{b^d}\big)^i\Big) \end{split}\]

Tramite i fatti precedentemente otteniamo che

\[\sum_{i = 0}^{\log_b n - 1} \big(\frac{a}{b^d}\big)^i = \frac{ \big(\frac{a}{b^d}\big)^{\log_b n} -1}{\frac{a}{b^d} -1} = \frac{ (\frac{n^{\log_b a}}{n^d}) -1}{\frac{a}{b^d} -1}\]

Arrivati a questo punto andiamo ad analizzare i vari casi a seconda di come \(d\) si rapporta con \(\log_b{a}\). In particolare abbiamo tre casi ad analizzare, e questi sono

  • Caso 1: (\(d = \log_b{a}\))

    In questo caso abbiamo che \(b^d = a\), e \(\frac{a}{b^d} = 1\). Troviamo quindi

    \[\begin{split} T(n) &\leq a^{\log_b n} \cdot T(1) + c \cdot n^d \cdot \Big( \sum\limits_{i = 0}^{\log_b n - 1} \big(\frac{a}{b^d}\big)^i\Big) \\ &= n^{\log_b a} + c \cdot n^d \cdot \Big( \sum_{i = 0}^{\log_b n - 1} (1)^i\Big) \\ &= n^{\log_b a} + c \cdot n^d \cdot \log_b n \\ &= \mathcal{O}(n^d \cdot \log n) \\ \end{split}\]

  • Caso 2: (\(d > \log_b{a}\))

    In questo caso...

  • Caso 3: (\(d < \log_b{a}\))

    In questo caso...

La moltiplicazione matriciale tra una matrice \(A = (a_{i, j})\) di dimensione \(n \times m\) e una matrice \(B = (b_{i, j})\) di dimensione \(m \times h\) è data da \(A \cdot B = (h_{i, j})\), con

\[ h_{i, j} := \sum_{h = 1}^m a_{i, h} b_{h, j}\]

L'algoritmo tradizionale che calcola tale prodotto ha una complessità di \(\theta(n^3)\). Anche in questo caso è possibile impostare un approccio Divide et Impera per moltiplicare due matrici. Questa volta tale approccio è basato sul fatto che il prodotto matriciale si può spezzare nel seguente modo:

\[\begin{bmatrix} A & B \\ C & D \\ \end{bmatrix} \cdot \begin{bmatrix} A & B \\ C & D \\ \end{bmatrix} = \begin{bmatrix} AE + BG & AF + BH \\ CE + DG & CF + DH \\ \end{bmatrix}\]

Dove \(A, B, C, D, E, F, G, H\) possono a loro volta essere delle matrici più piccole. Con questo approccio ricorsivo, se definiamo con \(m\) la dimensione della matrice, otteniamo la seguente ricorrenza

\[T(m) = 8 \cdot T\Big(\frac{m}{4}\Big) + \theta (m)\]

che, risolta tramite il teorema master, ci da \(T(m) = \theta (n^3)\).

Anche se può sembrare nuovamente di non aver ottenuto nessun vero progresso da un punto di vista di complessità, esiste, in modo analogo a quanto visto per le moltiplicazioni, un modo ingegnioso per ridourre le moltiplicazioni matriciali da 8 a 7, ottenendo così una complessità \(o(n^3)\). Questo metodo fu scoperto da Stressen.

Al fine di studiare gli algoritmi di crittografia durante il corso necessitiamo di poter calcolare il valore \(x^y \mod N\) per valori di \(x, y\) e \(N\) che possono avere svariate centinaia di bit.

Un primo approccio per il calcolo di \(x^y \mod N\) potrebbe essere quello di eseguire le operazioni intermedie modulo \(N\), moltiplicando ogni volta per \(x\). Detto altrimenti, avremmo il seguente schema

\[\underbrace{x \mod N \to x^2 \mod N \to ... \to x^y \mod N}_{y \text{ ripetizioni }}\]

Così facendo però, se \(y\) è molto grande, andremmo ad eseguire un numero di moltiplicazioni pari a \(O(y \cdot |x|)\), dove con \(|x|\) indichiamo il numero di bit necessari per memorizzare il numero \(x\). Notiamo che tale numero è esponenziale rispetto alla grande dell'input. Dobbiamo quindi trovare un algoritmo più efficiente.

Un'idea migliore è invece quella di partire con \(x\) e elevare al quadrato in modo ripetuto andando poi a semplificare modulo \(N\). Così facendo otteniamo il seguente schema

\[x \mod N \to x^2 \mod N \to x^4 \mod N \to ... \to x^{2 \cdot \lfloor \log_2 y \rfloor }\]

Una volta che abbiamo calcolato tutte queste queste potenze, per ottenere \(x^y \mod N\) dobbiamo moltiplicare solamente un determinato numero di queste potenze tra loro: quelle che corrispondono ad un \(1\) nella rappresentazione binaria di \(y\). L'algoritmo trovato è quindi il seguente

Notiamo poi che ciascuna moltiplicazione richiede \(O(\log^2N)\) operazioni e, con questo schema, dobbiamo fare solamente \(\log y = O(\log N)\) moltiplicazioni. Segue quindi che la complessità dell'algoritmo è \(O(n^3)\), dove \(n\) è il massimo numero di bit necessari per memorizzare \(x\), \(y\) e \(N\).

È possibile trovare una implementazione di tale algoritmo nel file modular_exponentiation.py.

Infine, per terminare la lezione, consideriamo il seguente problema: dati due interi \(a\) e \(b\), vogliamo trovare il più grande intero che divide sia \(a\) che \(b\). Questo intero è conosciuto anche con il nome di massimo comun divisore di \(a\) e \(b\), o, in simboli, \(MCD(a, b)\).

Un algoritmo che risolve questo problema risale ai tempi di Euclide, intorno a 2000 anni fa. L'algoritmo si basa sulla ripetuta divisione e può essere descritto come segue

Andiamo adesso ad argomentare la correttezza e la complessità di tale algoritmo

La correttezza dell'algoritmo segue dalla seguente proposizione

Proposizione (Eucliud's Rule): Se \(x\) e \(y\) sono interi positivi con \(x \geq y\), allora

\[MCD(x, y) = MCD(x \mod y, y)\]

Al fine di dimostrare la proposizione ci basta dimostrare che

\[MCD(x, y) = MCD(x - y, y)\]

Infatti, se vale questo risultato allora, iterando il ragionamento, vale anche \(MCD(x, y) = MCD(x - k \cdot y, y)\), con \(k\) intero. La proposizione da noi richiesta segue dal fatto che

\[x \mod y = x - k^* \cdot y, \,\,\,\, k^* \in \mathbb{Z}\]

Dimostrazione: Notiamo che se \(a\) divide \(x\) e \(y\), allora \(a\) divide anche \(x - y\) e \(y\). Ma allora i divisori di \(x\) e \(y\) sono gli stessi dei divisori di \(x - y\) e \(y\). In particolare quindi anche il massimo comune divisore è lo stesso tra la coppia di numeri.

\[\tag*{$\blacksquare$}\]

Per argomentare la complessità computazionale dell'algoritmo utilizziamo il seguente lemma

Lemma: Se \(a \geq b\), allora \(a \mod b < a/2\).

Dimostrazione: Spezziamo la dimostrazione analizzando i due possibili casi

  • Se \(b \leq a/2\), allora banalmente

    \[a \mod b < b \leq a/2\]

  • Se \(b > a/2\), allora ricordando che \(b \leq a\), abbiamo che

    \[a \mod b = a - b < a/2\]

\[\tag*{$\blacksquare$}\]

Notiamo che dal lemma segue che ogni volta che eseguiamo due chiamate, almeno uno dei due interi si è dimezzato. Ma allora, se inizialmente \(a\) e \(b\) sono interi da \(n\) bit, possiamo avere al più \(2n\) chiamate ricorsive. Dato che ogni chiamata ricorsiva consiste in una divisione, che richiede \(O(n^2)\), in totale il tempo di esecuzione è \(O(n^3)\).