Vai al contenuto
Topos Uranos
Italiano

Algoritmi numerici

Riepilogo

Le lezioni precedenti dell’unità hanno insegnato a scrivere cicli e a ragionare su di essi con un invariante e una variante; questa lezione applica quel ragionamento a quattro algoritmi numerici classici, più antichi dei computer e ancora in uso. Il test di primalità per divisioni di prova, che si ferma a n\sqrt{n} per una ragione che si dimostra, e la cui condizione di arresto va in overflow se la si scrive con trascuratezza. L’algoritmo di Euclide per il massimo comune divisore, con il suo invariante, la sua terminazione, la sua trappola con i numeri negativi e il suo confronto con std::gcd. Il cambio di base per divisioni successive e lo schema di Horner, che ricostruisce un numero a partire dalle sue cifre e valuta un polinomio con lo stesso schema. E la bisezione, che approssima una radice dividendo un intervallo in due, garantita dal teorema dei valori intermedi, con un numero di iterazioni che si calcola in anticipo. Ogni algoritmo è presentato come un oggetto matematico: con una dimostrazione che fa ciò che dice, una dimostrazione che termina e un primo conteggio dei suoi passi.

Obiettivi di apprendimento

  1. Implementare il test di primalità per divisioni fino a n\sqrt{n} e giustificare perché basta arrivare a n\sqrt{n}.
  2. Implementare l’algoritmo di Euclide, enunciarne l’invariante (il massimo comune divisore non cambia) e confrontarlo con std::gcd.
  3. Convertire un numero in un’altra base per divisioni successive e valutare una rappresentazione posizionale, o un polinomio, con lo schema di Horner.
  4. Approssimare una radice con la bisezione fino a una tolleranza, e calcolare quante iterazioni richiede.

Divisori e numeri primi

Divisioni di prova

Un intero n≥2n \geq 2 è primo se i suoi unici divisori positivi sono 1 e nn; se ne ha un altro, è composto; 0 e 1 non sono né l’uno né l’altro. La definizione suggerisce subito un algoritmo: provare ogni candidato dd tra 2 e n−1n - 1 e dichiarare composto il numero non appena uno lo divide. È il ciclo del minimo divisore della lezione Cicli, ed è corretto. Tuttavia, per un primo di dieci cifre esegue quasi dieci miliardi di divisioni, e quasi tutte sono inutili, per la ragione seguente.

ProposizioneBasta provare fino a n\sqrt{n}

Se n≥2n \geq 2 è composto, ha un divisore dd con 2≤d≤n2 \leq d \leq \sqrt{n}.

Dimostrazione

  1. n=a⋅b,2≤a≤bn = a \cdot b, \quad 2 \leq a \leq b

    Se nn è composto, si scrive come prodotto di due fattori, nessuno uguale a 1 né a nn. Si chiama aa il minore dei due.

  2. a2≤a⋅b=na^2 \leq a \cdot b = n

    Moltiplicando a≤ba \leq b per il numero positivo aa si ottiene a2≤ab=na^2 \leq ab = n, da cui a≤na \leq \sqrt{n}: il fattore minore di ogni fattorizzazione sta a sinistra di n\sqrt{n}.

  3. d≤n  ⟺  nd≥nd \leq \sqrt{n} \iff \frac{n}{d} \geq \sqrt{n}

    I divisori si presentano, dunque, in coppie (d,n/d)(d, n/d) riflesse attorno a n\sqrt{n}: con n=36n = 36, le coppie sono (2,18)(2, 18), (3,12)(3, 12), (4,9)(4, 9) e (6,6)(6, 6). Se nessun candidato fino a n\sqrt{n} divide nn, non lo divide nemmeno nessuno oltre, perché il suo compagno sarebbe comparso prima.

Ne segue che il numero di divisioni scende da circa nn a circa n\sqrt{n}: per nn vicino a due miliardi, da due miliardi a quarantaseimila. Il programma usa lo schema del flag della lezione Progettare cicli corretti: isPrime comincia vero e il ciclo lo spegne se trova un divisore.

Programma in C++
#include <iostream>#include <print>​int main(){    int n{};    std::cin >> n;    bool isPrime{n >= 2};    for (int d{2}; isPrime && d <= n / d; ++d) {        if (n % d == 0) {            isPrime = false;        }    }    if (isPrime) {        std::println("{} è primo", n);    } else {        std::println("{} non è primo", n);    }}
Input del programma
2147483647
Output
2147483647 è primo

Tre decisioni meritano un commento. Il flag si inizializza con n >= 2, in modo che 0, 1 e i negativi siano dichiarati non primi senza entrare nel ciclo: dimenticare questo caso è l’errore tipico, perché con n=1n = 1 il ciclo non esegue nessuna iterazione e il programma dichiarerebbe primo l’1. La condizione del for ha due parti: il flag, che ferma il ciclo al primo divisore, e il limite. E il limite è scritto d <= n / d, che, per interi positivi, equivale a d2≤nd^2 \leq n (la divisione intera tronca, ma dd non supera il quoziente troncato esattamente quando non supera n/dn/d, perché dd è intero). Il numero dell’input, 231−12^{31} - 1, è il massimo int ed è primo; il programma lo decide con 46 33946\,339 resti, quelli dei candidati da 2 a 46 34046\,340.

L’invariante del ciclo è: nessun intero da 2 a d−1d - 1 divide nn, finché il flag resta vero. Vale all’ingresso, quando l’intervallo da 2 a 1 è vuoto; si conserva, perché se dd non divide nn l’intervallo cresce di un numero che non lo divide; e all’uscita con il flag vero, l’invariante e la condizione falsa, d2>nd^2 > n, dicono insieme che nessun candidato fino a n\sqrt{n} divide nn, che è ciò che la proposizione richiede. La variante è la quantità di candidati che restano da provare, cioè di interi tra dd e n\sqrt{n}: decresce di uno a ogni iterazione e non scende sotto zero.

La condizione che va in overflow

Il modo naturale di scrivere il limite è d * d <= n, ed è un errore. Con l’input precedente, il ciclo arriva a d=46 341d = 46\,341, e allora il prodotto 46 3412=2 147 488 28146\,341^2 = 2\,147\,488\,281 supera il massimo di int.

Programma in C++
#include <iostream>#include <print>​int main(){    int n{};    std::cin >> n;    bool isPrime{n >= 2};    for (int d{2}; isPrime && d * d <= n; ++d) {        if (n % d == 0) {            isPrime = false;        }    }    std::println("Primo? {}", isPrime);}
Input del programma
2147483647
Output
main.cpp:9:33: runtime error: signed integer overflow: 46341 * 46341 cannot be represented in type 'int'

Il programma è terminato con il codice 1.

Il sanitizer del comportamento indefinito ferma il programma sul prodotto. Senza di esso, l’overflow con segno è comportamento indefinito, e lo standard non dice che cosa accade. Ciò che si osserva con GCC su Linux, con e senza ottimizzazione, è istruttivo: il prodotto fa il giro e risulta negativo o minore di nn, e poiché nessun int supera il massimo, la condizione d * d <= n non smette più di essere vera; il ciclo continua a provare candidati per circa due miliardi di iterazioni, finché arriva a d=nd = n, che divide nn, e il programma dichiara composto un numero primo. Un errore di overflow non si manifesta come un overflow, ma come una risposta falsa. La forma d <= n / d non può traboccare, perché divide invece di moltiplicare; la lezione Bit e libreria matematica offriva un’altra via d’uscita, calcolare il limite una volta con std::sqrt, ma allora l’esattezza del limite dipende dall’arrotondamento della virgola mobile, e la divisione intera non ha questo problema.

La lezione «Numeri primi e teorema fondamentale dell’aritmetica», del corso Algebra e Geometria I, dimostra che ogni intero maggiore di 1 è prodotto di primi in un solo modo; le divisioni di prova sono il modo più diretto di trovare quel prodotto, dividendo nn per ogni divisore trovato e proseguendo con il quoziente.