Assembly Avanzato con NASM

Capitolo 19: Programmazione del coprocessore matematico

19.1 Visualizzazione di numeri in virgola mobile sullo schermo

Vediamo un primo esempio che si riferisce alla visualizzazione di numeri in virgola mobile, in base 10, sullo schermo; l'idea di fondo si basa su concetti già illustrati nel tutorial Assembly Base.
Abbiamo visto che qualsiasi informazione da visualizzare sullo schermo, compresi i numeri, deve essere prima convertita in una sequenza di codici ASCII; ciò è dovuto al fatto che l'hardware grafico del computer, in modalità testo, lavora sui simboli associati a tali codici, convertendoli in una mappa di pixel da accendere sul monitor.
Nel caso dei numeri interi, che risultano memorizzati sul computer in formato binario, si può fare ricorso al cosiddetto metodo delle divisioni successive. Dividendo ripetutamente un numero in base b1 per una diversa base b2, i resti di ogni divisione forniscono le cifre del numero stesso nella nuova base b2, a partire dalla meno significativa; il procedimento termina quando il quoziente diventa nullo. Se b2=10, ricordando che il resto di una divisione è sempre inferiore al divisore, tutti i resti che si ottengono sono singole cifre comprese tra 0 e 9; a questo punto, indicando una generica cifra in base 10 con il simbolo c, si ha: 
ASCII(c) = c + ASCII(0)
Questo metodo funziona perché i codici ASCII delle cifre da 0 a 9 sono consecutivi e contigui; infatti, essi vanno da 48 (simbolo '0') a 57 (simbolo '9'), per cui, ad esempio: 
ASCII(3) = 3 + ASCII(0) = 3 + 48 = 51 = '3'
Applicando questi concetti ai numeri in virgola mobile, possiamo schematizzare il seguente procedimento: Alla parte intera possiamo applicare direttamente il metodo delle divisioni successive; per la parte frazionaria ciò non è possibile in quanto sappiamo che, ad esempio, il numero frazionario 0.320421 equivale a: \[ 3 \cdot 10^{-1} + 2 \cdot 10^{-2} + 0 \cdot 10^{-3} + 4 \cdot 10^{-4} + 2 \cdot 10^{-5} + 1 \cdot 10^{-6} \] e non a: \[ 3 \cdot 10^5 + 2 \cdot 10^4 + 0 \cdot 10^3 + 4 \cdot 10^2 + 2 \cdot 10^1 + 1 \cdot 10^0 \] Alla parte frazionaria bisogna applicare il metodo delle moltiplicazioni successive. Moltiplicando un numero frazionario in base b1 per la sua stessa base, facciamo scorrere le sue cifre di un posto verso sinistra; la cifra che "trabocca" dalla parte frazionaria, va quindi a posizionarsi alla sinistra della virgola, producendo un numero intero.
Ad esempio, se vogliamo 4 cifre dopo la virgola, non dobbiamo fare altro che moltiplicare per 10000; nel caso della precedente parte frazionaria 0.320421, si ha: \[ 0.320421 \cdot 10000 = 3204.21 \] A questo punto, non ci resta che convertire 3204 in formato ASCII e visualizzarlo sullo schermo preceduto dal punto decimale.

Naturalmente, dobbiamo essere in grado di gestire varie situazioni, come il segno del numero o la presenza di zeri significativi subito dopo la virgola, come nel caso di 123.000326; la Figura 19.1 illustra un programma di esempio riferito ai numeri reali in formato Single Real a 32 bit. La costante FDIGITS indica quante cifre vogliamo dopo la virgola; con tale costante viene inizializzato il moltiplicatore 10FDIGITS per la parte frazionaria.

La variabile reale da visualizzare viene chiamata float_var; in tale variabile viene caricato il numero \(\pi\) attraverso le istruzioni: In ogni caso, si possono commentare tali due istruzioni e inizializzare direttamente float_var con un numero reale a piacere, purché sia un Single Real non troppo grande in valore assoluto.

Per la separazione della parte intera dalla parte frazionaria, impostiamo la Control Word in modo da avere il modo di approssimazione "truncate"; a tale proposito, i suoi due bit in posizione 10 e 11 vengono posti a 11b attraverso le istruzioni: La parte intera viene ottenuta caricando float_var in ST(0) e salvandola poi in integer_part con FISTP, in modo da ottenere un intero per troncamento della parte frazionaria.
La parte frazionaria viene ottenuta sottraendo a float_var la sua parte intera; il risultato così ottenuto viene salvato in significand. La variabile significand viene moltiplicata per 10FDIGITS; il risultato viene poi salvato in fractional_part con FISTP, in modo da ottenere un intero per troncamento della parte frazionaria. Ovviamente, dobbiamo svolgere questo compito con la FPU in quanto la CPU, non essendo in grado di gestire numeri in floating point, tratterebbe significand come un numero intero! A questo punto, il procedimento diventa molto semplice in quanto si tratta di applicare il metodo delle divisioni successive a integer_part e a fractional_part; la procedura fill_buffer esegue tale compito chiamando due volte l'altra procedura num_toascii.
Si ricordi che dividendo EAX per EBX, si ottiene il quoziente in EAX e il resto in EDX (è importantissimo quindi azzerare prima EDX); in questo modo possiamo prelevare da DL la cifra da convertire, mentre EAX contiene il dividendo pronto per la successiva divisione per EBX.
Il secondo loop in num_toascii serve per sapere se la parte frazionaria appena convertita in ASCII è formata da almeno FDIGITS cifre; ciò è necessario per evitare che numeri del tipo 0.0002315 con, ad esempio, FDIGITS=5, producano la stringa "23" anziché "00023".

Ogni cifra convertita in ASCII viene salvata nella stringa float_buff; in tale stringa vengono posizionati anche il segno e il punto decimale.
Il metodo delle divisioni successive, come sappiamo, produce una stringa numerica disposta al contrario; di conseguenza, è necessario procedere alla sua inversione. Tale compito viene svolto dalla procedura reverse, la quale ha bisogno di sapere dove inizia e dove finisce la stringa da invertire all'interno di float_buff.

Come abbiamo visto nei precedenti capitoli, un Single Real ha una precisione pari a 6-7 cifre dopo la virgola; se abbiamo bisogno di una precisione maggiore, possiamo ricorrere ai Long Real, i quali garantiscono 15-16 cifre significative dopo la virgola. Si tenga presente però che, maggiore è la parte intera, meno precisa sarà la parte frazionaria.
Per esercizio, si può provare a convertire il programma di Figura 19.1 in modo che gestisca i Long Real a 64 bit; chiaramente, in un caso del genere, la parte più impegnativa è quella relativa all'applicazione del metodo delle divisioni successive. Per svolgere tale compito, ci si può servire degli algoritmi presentati nel tutorial Assembly Base per la divisione di numeri interi di ampiezza arbitraria; in alternativa, si può sfruttare l'istruzione FPREM della FPU per calcolare il resto di una divisione intera.

Nel precedente capitolo è stato spiegato che anche l'istruzione FXTRACT può essere usata per convertire in ASCII un numero in floating point; come sappiamo, tale istruzione legge il valore presente in ST(0) e salva la caratteristica in ST(0) e la mantissa in ST(1). Conoscendo la caratteristica, sappiamo di quanti posti dobbiamo far scorrere, verso destra o sinistra del punto decimale, i bit della mantissa; a questo punto abbiamo a disposizione la parte intera e la parte frazionaria e possiamo quindi procedere come già descritto in precedenza.

19.2 Gestione delle eccezioni generate dalla FPU

Nel Capitolo 17 è stato spiegato che la gestione delle eccezioni generate dalla FPU ha avuto una storia piuttosto tormentata; tutto questo nonostante le raccomandazioni della Intel, che invitavano i progettisti di hardware ad utilizzare il metodo classico del PIC 8259A.
Sui PC di classe XT è presente un unico PIC, con tutti gli 8 ingressi già occupati; per questo motivo, l'uscita INT di una eventuale FPU 8087, risulta collegata all'ingresso NMI della CPU. La ISR che gestisce la NMI (vettore 02h) verifica se la richiesta di interruzione è dovuta ad un problema hardware o ad una eccezione della FPU; in questo secondo caso, viene chiamata un'altra apposita ISR che rappresenta il gestore di eccezione o Exception Handler. Il meccanismo appena descritto è all'origine di quello che oggi viene definito MS-DOS Compatibility Mode.
Sui PC di classe AT, le CPU delle famiglie 80286 e 80386 sono dotate di un apposito ingresso ERROR#; a tale ingresso giungono le richieste di interruzione dalla corrispondente uscita ERROR# di una eventuale FPU esterna. Quando la CPU riceve una richiesta di interruzione dal suo ingresso ERROR#, genera una apposita IRQ16 per invocare un gestore di eccezione; questo nuovo metodo prende il nome di Native Mode.
La necessità di garantire la compatibilità con i PC di classe XT, ha però portato all'uso di un meccanismo differente, che sfrutta la presenza di due PIC sui PC di classe AT. L'ingresso ERROR# della CPU rimane inutilizzato, mentre l'uscita ERROR# della FPU risulta collegata all'ingresso 5 del PIC Slave (IRQ13); a tale IRQ viene associato il vettore 75h. La ISR del vettore 75h, sempre per motivi di compatibilità, provvede a chiamare a sua volta il vettore 02h che elabora le NMI.
Un ulteriore cambiamento si verifica con l'arrivo delle CPU 80486 DX e superiori, dotate di FPU integrata al loro interno; nel registro di controllo CR0 di tali CPU, il bit in posizione 5, denominato NE (Numeric Error), permette di scegliere tra MS-DOS Compatibility Mode (NE=0) e Native Mode (NE=1). Su queste CPU sono presenti due ulteriori pin denominati FERR# (Floating-point ERRor) e IGNNE# (IGNore Numeric Error), connessi agli omonimi pin della FPU integrata. In Native Mode, viene utilizzato lo stesso meccanismo della IRQ16 descritto in precedenza; il pin FERR# in questo caso sostituisce il pin ERROR# delle precedenti famiglie di CPU. In MS-DOS Compatibility Mode, il pin FERR# risulta connesso alla IRQ13 del PIC, associata come al solito al vettore 75h; il pin IGNNE# viene usato per permettere alla FPU di eseguire ulteriori istruzioni numeriche all'interno dell'Exception Handler senza che vengano generate altre eccezioni innestate (cosa che provocherebbe una chiamata ricorsiva allo stesso gestore di eccezione, con il rischio di un loop infinito).

Le considerazioni appena esposte ci fanno capire che la scrittura di un gestore di eccezione numerica in Modalità Reale, può risultare piuttosto complicata; infatti, il metodo da seguire viene influenzato dal modello di CPU che si sta utilizzando e dal sistema operativo. Anche gli emulatori di sistemi operativi e le macchine virtuali possono risentire di questa situazione; in certi casi, può capitare che le funzionalità descritte in precedenza non risultino neppure implementate!
Nel programma di esempio mostrato in Figura 19.2, si suppone di operare su una piattaforma hardware dotata di CPU 80486 DX o superiore; il bit NE del registro CR0 viene posto a zero (MS-DOS Compatibility Mode) e si provvede ad intercettare la IRQ13 attraverso una apposita ISR per il vettore 75h. L'esempio risulta funzionante dal prompt di MS-DOS fornito da Windows, anche quando il sistema operativo gira sotto VirtualBox; gli emulatori per MS-DOS, come DOSBox o DOSEmu, risultano inutilizzabili in quanto forniscono un supporto parziale o nullo dei meccanismi di gestione delle eccezioni numeriche. Il programma provvede innanzi tutto a porre a zero il bit NE del registro di controllo CR0; si tenga presente che tale registro è presente anche sulle CPU 80386, ma il bit NE rimane inutilizzato in quanto, come abbiamo visto in precedenza, con questo modello di microprocessore è disponibile solo il MS-DOS Compatibility Mode. Lo stesso registro CR0 dispone anche del bit a sola lettura EM (Emulation) in posizione 2, che indica la presenza (EM=0) o meno (EM=1) di una FPU.
Successivamente, il programma installa la nuova ISR per il vettore 75h; il passo successivo consiste nell'inizializzare la FPU con FINIT e nel porre a zero i primi sei bit del Control Register in modo da abilitare la generazione delle relative eccezioni.
Il programma mostra poi lo stato dei flags di eccezione nello Status Register e avvia un loop nel quale, tramite un menu, si può scegliere se generare una eccezione, azzerare i flags di eccezione o uscire.
Per la Invalid Operation Exception si effettua il calcolo \(\sqrt{-1}\) con le istruzioni: Per la Zero Divide Exception si effettua il calcolo \({+1.0} / {+0.0}\) con le istruzioni: Per la Precision Exception si salva \(\pi\) come Word Integer con le istruzioni: Per esercizio si può provare a far generare alla FPU ulteriori eccezioni; ad esempio, la divisione tra un numero molto piccolo e uno molto grande produce un risultato denormalizzato.
Si ricordi che l'eccezione viene effettivamente generata solo quando la FPU incontra l'istruzione successiva (a quella che ha prodotto l'errore) o una istruzione FWAIT (che è il metodo utilizzato nell'esempio di Figura 19.2).

Il gestore di eccezione non fa altro che mostrare lo stato dei flags nello Status Register e azzerare i flags stessi con l'istruzione FNCLEX; come sappiamo, la N del mnemonico sta per No Wait e ciò è fondamentale per evitare che venga generata una ulteriore eccezione (come accadrebbe con FCLEX) con conseguente chiamata ricorsiva dell'Exception Handler.
Un ulteriore importante compito da svolgere attraverso il gestore di eccezione consiste nell'usare l'istruzione FNSAVE per salvare l'intero stato della FPU, compreso il contenuto dello stack; in questo modo abbiamo a disposizione tutti i dettagli relativi al problema che ha provocato l'eccezione. Si ricordi che FNSAVE provvede anche a reinizializzare la FPU, con conseguente mascheramento di tutti i flags di eccezione nel Control Register.

All'interno del gestore di eccezione si può notare anche l'istruzione: 
   out      FERR, al                     ; clear FERR# interrupt
Scrivendo un byte nella porta FP_IRQ=00F0h, si azzera il segnale FERR# che ha inviato la richiesta di interruzione a causa di una eccezione della FPU; tale istruzione comunque è necessaria in Native Mode, mentre in MS-DOS Compatibility Mode dobbiamo procedere come al solito ad inviare un End Of Interrupt al PIC (in questo caso, trattandosi della IRQ13, dobbiamo inviare l'EOI al PIC Slave e al PIC Master).

19.3 Grafici di funzioni matematiche con la libreria VESALIB

Vediamo un esempio che consiste nel tracciare il grafico della funzione \(y = \sin(x)\); la FPU ci facilita il compito mettendoci a disposizione l'istruzione FSIN, ma il vero problema che dobbiamo affrontare consiste nell'adattare il grafico stesso alla modalità video grafica del computer.
La libreria VESALIB utilizza il pixel come unità di misura; di conseguenza, se non prendiamo le opportune precauzioni, otteniamo un grafico pressoché invisibile. Si ricordi, infatti, che \(\forall \ x \in \Bbb R\) si ha \(-1 \le \sin(x) \le +1\); di conseguenza, il grafico che si ottiene risulta completamente schiacciato, in quanto oscilla tra -1 pixel e +1 pixel!
La soluzione a questo problema consiste nell'applicare un opportuno fattore di scala, sia all'asse x, sia all'asse y.

Partiamo dall'asse x e supponiamo di volere un grafico la cui ascissa vada da -300 pixel a +300 pixel; inoltre, vogliamo che la variabile indipendente \(x\) di \(\sin(x)\) vada da \(-{2\pi}\) a \(+{2\pi}\). Si vede allora che l'incremento di 1 pixel corrisponde al seguente incremento \(\Delta{x}\) di \(x\): \[ \Delta{x} = \frac {{+2\pi} - (-2\pi)} {+300 - (-300)} = \frac {4\pi} {600} \approx +0.02093 \] In sostanza, all'interno del loop che traccia il grafico di \(y = \sin(x)\), ad ogni iterazione l'ascissa del pixel si incrementa di 1, mentre l'ascissa \(x\) si incrementa di circa +0.02093.

Per l'asse y il procedimento è ancora più semplice in quanto la FPU ci mette a disposizione l'istruzione FSCALE; come sappiamo, questa istruzione calcola ST(0)*2ST(1).
Supponiamo allora di volere per l'asse y un fattore di scala pari a 128, che corrisponde a 27 (ovviamente, dobbiamo scegliere un numero esprimibile come potenza intera di 2); dopo aver caricato il valore 7 in ST(0), facciamo partire il loop e, ad ogni iterazione, carichiamo in ST(0) l'ascissa \(x\) (con il fattore di scala 7 che quindi passa in ST(1)) e poi eseguiamo le istruzioni: Il risultato così ottenuto rappresenta l'ordinata del pixel da visualizzare; trattandosi di un numero in virgola mobile, lo dobbiamo arrotondare ad un intero con l'istruzione FISTP. Tale istruzione è importantissima in quanto estrae dallo stack il contenuto di ST(0) e lo salva in una variabile intera; in questo modo, quando alla successiva iterazione viene caricata in ST(0) la nuova ascissa \(x + \Delta{x}\), il fattore di scala 7 verrà spostato sempre in ST(1).

Mettendo assieme le considerazioni appena esposte, otteniamo il codice sorgente FNCGRAPH.ASM, illustrato in Figura 19.3 (questo esempio e quelli che seguono, utilizzano come scheletro il programma VLIBTEST.ASM presentato alla fine del Capitolo 11). Il grafico di \(y = \sin(x)\) viene tracciato dalla procedura plot_sinx, la quale segue le convenzioni del linguaggio C ed ha la seguente dichiarazione: 
plot_sinx(int pcol, int x1i, int x2i, float x1f, float x2f, float ysf)
Supponendo allora di aver definito le variabili: la chiamata a plot_sinx averrà in questo modo: In questo esempio, la procedura plot_sinx si occupa anche del calcolo di \(\Delta{x}\); in alternativa, possiamo svolgere tale calcolo nel programma principale e poi passare il risultato come parametro a plot_sinx (dopo aver modificato opportunamente la struttura della procedura).

La Figura 19.4 mostra l'output di FNCGRAPH.EXE. Il programma è interamente parametrizzato, per cui può essere eseguito in tutte le modalità video grafiche a 256 colori supportate dall'hardware; a tale proposito, basta modificare la costante simbolica VESA_MODE in FNCGRAPH.ASM, assegnandole uno tra i codici 0101h, 0103h, 0105h, 0107h.

Al posto di FSIN, nella procedura plot_sinx possiamo anche usare FCOS per avere il grafico di \(y = \cos(x)\). Il programma funziona anche con l'istruzione FPTAN, nonostante la presenza di punti di discontinuità (dove \(y = \tan(x)\) tende a \(\infty\)); la libreria VESALIB, infatti, salta i pixel che cadono fuori dallo schermo.
Se si decide di usare l'istruzione FPTAN, si deve ricordare che la FPU, dopo aver calcolato \(y = \tan(x)\), carica il valore +1.0 in ST(0), in modo da facilitare l'eventuale calcolo di \(\cot(x) = 1 / \tan(x)\); prima di eseguire FSCALE, bisogna quindi estrarre il +1.0 dallo stack e salvarlo da qualche parte.

19.3.1 Grafico di un segmento di retta

In una libreria di procedure grafiche, si ha spesso bisogno di un algoritmo per tracciare un segmento di retta. Il caso delle rette orizzontali o verticali è banale, in quanto basta far variare una delle coordinate tenendo fissa l'altra; la libreria VESALIB fornisce le due procedure vesa_puthline256 e vesa_putvline256 per svolgere tale compito. Occupiamoci quindi del caso più interessante, che riguarda una retta obliqua (o, più in generale, non verticale).
Come sappiamo dalla Geometria, una retta non verticale del piano ha equazione cartesiana: \[ y = mx + q \] dove \(m\) è il coefficiente angolare e \(q\) è una costante (pari a zero se la retta passa per l'origine); se \(P_1(x_1, y_1)\) e \(P_2(x_2, y_2)\) sono due punti distinti della retta, si ha, come noto: \[ m = \frac {y_2 - y_1} {x_2 - x_1} \] Si vede subito che questa formula non è applicabile ad una retta verticale \((x_1 = x_2)\) in quanto si otterrebbe \(m = \infty\); per una retta orizzontale \((y_1 = y_2)\), invece, si ha \(m = 0\) (pendenza nulla) e equazione \(y = q\).

Nel caso dell'equazione \(y = \sin(x)\), per il fattore di scala, si è fatto ricorso alla moltiplicazione ST(0)*2ST(1), che con la FPU è piuttosto pesante in termini di cicli macchina; naturalmente, esistono appositi algoritmi che permettono di seguire una strada più efficiente.
Nel caso dell'equazione \(y = mx + q\), per fortuna, è possibile eliminare la moltiplicazione \(mx\); a tale proposito, si può osservare che quando l'ascissa \(x\) passa da \(x_i\) a \(x_{i+1} = x_i + \Delta{x}\) (con \(\Delta{x} = 1\) pixel), l'ordinata \(y\) passa da \(y_i = mx_i + q\) a \(y_{i+1} = mx_{i+1} + q\). Si ha allora: \[ y_{i+1} = mx_{i+1} + q = m(x_i + \Delta{x}) + q = (mx_i + q) + m\Delta{x} \] Ma \(mx_i + q = y_i\) e \(m\Delta{x} = m \cdot 1 = m\); di conseguenza: \[ y_{i+1} = y_i + m \] In sostanza, partendo dal punto \(P_1(x_1, y_1)\), si avvia un loop nel quale, ad ogni iterazione, la \(x\) si incrementa di \(\Delta{x} = 1\), mentre la \(y\) si incrementa di \(m\); ovviamente, \(y\) è un numero in floating point che deve essere arrotondato ad un intero con FISTP. Si ottiene allora la seguente procedura (sintassi NASM), la quale richiede obbligatoriamente \(x_1 \lt x_2\): Chiamando questa procedura secondo le convenzioni C, con P1(-200,+180), P2(+80,-200) e colore 9, si ottiene l'output mostrato in Figura 19.5. Si tenga presente che la procedura put_line usa il pixel come unità di misura (scala 1:1); se si vuole usare put_line in combinazione con plot_sinx (ad esempio, per tracciare la tangente a \(y = \sin(x)\) in un suo punto), si deve applicare ad entrambe le funzioni matematiche lo stesso fattore di scala.

La procedura put_line può essere notevolmente velocizzata eliminando del tutto il ricorso alla FPU e utilizzando solamente numeri interi con la CPU; a tale proposito, si può ricorrere al famoso algoritmo di Bresenham (si veda, ad esempio, Wikipedia).

19.4 Approssimazione di un'onda quadra con lo sviluppo in serie di Fourier

Affrontiamo ora un problema più complesso che consiste nell'approssimare un'onda quadra con lo sviluppo in serie di Fourier; si tratta di un procedimento che richiede una notevole capacità di calcolo e quindi rappresenta un test molto impegnativo per la FPU.

Una funzione reale di variabile reale \(y = f(x)\) si definisce "periodica di periodo \(T\)" se \(f(x + kT) = f(x)\), con \(k \in \Bbb Z\); in sostanza, i valori \(y\) che la funzione assume si ripetono identicamente ogni \(T\) unità di misura (radianti, secondi, etc).
La funzione \(y = \sin(x)\), ad esempio, è periodica di periodo \(2\pi\) radianti in quanto si ha \(\sin(x + 2k\pi) = \sin(x)\); la funzione \(y = \tan(x)\), invece, è periodica di periodo \(\pi / 2\) radianti in quanto si ha \(\tan(x + (k\pi) / 2) = \tan(x)\).
Dalle considerazioni appena esposte risulta evidente che, una funzione periodica di periodo \(T\), può essere studiata in modo completo esaminando il suo comportamento in un qualunque intervallo di ampiezza \(T\).

Consideriamo ora le seguenti coppie di funzioni: \[ \cos\bigg(\frac {2k\pi x} {T} \bigg), \quad \sin\bigg(\frac {2k\pi x} {T} \bigg), \qquad k = 0, 1, 2, \dots, n \qquad (n \in \Bbb N) \] Per \(k = 0\) si ottiene \(\cos(0) = 1\) e \(\sin(0) = 0\), qualunque sia \(x\); una funzione che assume un valore costante \(c\) dappertutto può essere considerata periodica di periodo \(T\) qualsiasi, in quanto in ogni intervallo di ampiezza \(T\) si comporta sempre allo stesso modo (vale sempre \(c\)).
Per \(k = 1\) otteniamo due funzioni periodiche di periodo \(T\); infatti, si ha ad esempio: \[ \cos\bigg[\frac {2\pi (x + T)} {T}\bigg] = \cos\bigg(\frac {2\pi x} {T} + \frac {2\pi T} {T}\bigg) = \cos\bigg(\frac {2\pi x} {T} + 2\pi\bigg) = \cos\bigg(\frac {2\pi x} {T}\bigg) \] Analogamente, per \(k = 2\) otteniamo due funzioni periodiche di periodo \(T / 2\), per \(k = 3\) otteniamo due funzioni periodiche di periodo \(T / 3\) e così via, sino a \(k = n\) che ci fornisce due funzioni periodiche di periodo \(T / n\).
Ricordando che la frequenza con cui si ripete un evento periodico è l'inverso del periodo, possiamo affermare che per \(k = 2\) la frequenza è doppia rispetto al caso \(k = 1\); analogamente, per \(k = 3\) si ha una frequenza tripla, per \(k = 4\) si ha una frequenza quadrupla e così via.

Dato che \(n\) può assumere un valore grande a piacere, abbiamo a che fare con infinite funzioni le quali hanno in comune il periodo \(T\), per cui anche una loro combinazione lineare a coefficienti costanti sarà una funzione periodica di periodo \(T\); si ottiene allora il seguente polinomio trigonometrico di ordine n: \[ \frac {a_0} {2} + \sum\limits_{k=1}^n \bigg[a_k \cos\bigg(\frac {2k\pi x} {T}\bigg) + b_k \sin\bigg(\frac {2k\pi x} {T}\bigg)\bigg] \qquad a_0, a_k, b_k \in \Bbb R \] (dove \(a_0 / 2\) rappresenta il coefficiente del termine associato a \(k = 0\)).

Facendo tendere \(n\) a \(+\infty\) si ottiene una serie trigonometrica.

Se la serie trigonometrica converge in un intervallo del tipo \([x, x + T]\), cioè, se per ogni \(x\) appartenente a tale intervallo si ottiene un numero \(y \in \Bbb R\), allora essa converge in tutto \(\Bbb R\) (perché è periodica) e la sua somma è una funzione periodica di periodo \(T\).

Il problema che ci interessa però è quello inverso: vogliamo sapere se, data una funzione \(y = f(x)\), periodica di periodo \(T\), sia possibile determinare i coefficienti \(a_0, a_k, b_k \in \Bbb R\) in modo da avere: \[ f(x) = \frac {a_0} {2} + \sum\limits_{k=1}^{+\infty} \bigg[a_k \cos\bigg(\frac {2k\pi x} {T}\bigg) + b_k \sin\bigg(\frac {2k\pi x} {T}\bigg)\bigg] \] Si dimostra che la risposta è affermativa sotto condizioni piuttosto larghe.

Se ciò accade, si dice che (escludendo \(a_0 / 2\)) il primo termine della serie \((k = 1)\) rappresenta la prima armonica o armonica fondamentale di \(f(x)\), il secondo termine \((k = 2)\) rappresenta la seconda armonica di \(f(x)\) e così via.

Osserviamo ora che, se una funzione \(y = f(x)\) è periodica di periodo \(T\), possiamo sempre assumere \(T = 2\pi\); infatti, basta effettuare il cambiamento di variabile \(x = (Tu) / (2\pi)\), per cui: \[ f(x) = f\bigg(\frac {Tu} {2\pi}\bigg) = \psi(u) \] A questo punto si vede subito che: \[ \psi(u + 2\pi) = f\bigg[\frac {T(u + 2\pi)} {2\pi}\bigg] = f\bigg(\frac {Tu} {2\pi} + T\bigg) = f\bigg(\frac {Tu} {2\pi}\bigg) = \psi(u) \] Abbiamo quindi una funzione \(y = f(x)\) periodica di periodo \(2\pi\) che può essere studiata in un qualsiasi intervallo di ampiezza \(2\pi\) come, ad esempio, \([{-\pi}, {+\pi}]\); ciò porta ad una prima semplificazione in quanto, per \(T = 2\pi\), la serie trigonometrica diventa: \[ \frac {a_0} {2} + \sum\limits_{k=1}^{+\infty} \bigg[a_k \cos(kx) + b_k \sin(kx)\bigg] \qquad a_0, a_k, b_k \in \Bbb R \] Data una funzione \(y = f(x)\), periodica di periodo \(2\pi\), se tale funzione è sviluppabile in serie trigonometrica nell'intervallo \([{-\pi}, {+\pi}]\) e in tale intervallo la serie stessa converge uniformemente (e quindi la \(f\) è ivi continua), allora si dimostra che i coefficienti \(a_0, a_k, b_k\) si ottengono dalle formule seguenti (dette formule di Eulero-Fourier): \[ a_k = \frac {1} {\pi} \int_{-\pi}^{+\pi} f(x) \cos(kx) dx, \qquad k = 0, 1, 2, \dots \] \[ b_k = \frac {1} {\pi} \int_{-\pi}^{+\pi} f(x) \sin(kx) dx, \qquad k = 1, 2, 3, \dots \] Il problema che vogliamo affrontare non è ancora risolto, in quanto abbiamo dato per scontato che la \(f(x)\) sia sviluppabile in serie trigonometrica; si nota subito però che, se la \(f(x)\) è generalmente continua e integrabile in valore assoluto nell'intervallo \([{-\pi}, {+\pi}]\), allora i due precedenti integrali esistono sicuramente e quindi si può costruire la serie trigonometrica associata ai coefficienti \(a_0, a_k, b_k\). Se ciò accade, i numeri \(a_0, a_k, b_k\) prendono il nome di coefficienti di Fourier, mentre la serie ad essi associata viene chiamata serie di Fourier.

Quando è possibile sviluppare una funzione in serie di Fourier?

Esiste un teorema dovuto a Dirichlet in base al quale, una funzione continua e monotòna a tratti nell'intervallo \([{-\pi}, {+\pi}]\) e ivi limitata, è sviluppabile in serie di Fourier in tutti i punti interni all'intervallo stesso, nei quali essa è continua; nei punti di discontinuità interni all'intervallo la serie converge verso la media aritmetica dei limiti sinistro e destro relativi a quei punti, mentre agli estremi \(-\pi\) e \(+\pi\) converge verso la media aritmetica dei limiti relativi a quei punti.
Una funzione è "continua e monotòna a tratti" in un intervallo, quando tale intervallo può essere decomposto in un numero finito di intervalli, privi a due a due di punti interni comuni, in ciascuno dei quali la funzione è continua e monotòna.

Il problema che vogliamo affrontare, ricade in uno dei casi particolari che portano ad ulteriori semplificazioni; a tale proposito, definiamo pari una funzione \(f(x)\) tale che \(f(-x) = f(x)\), mentre definiamo dispari una funzione \(f(x)\) tale che \(f(-x) = -f(x)\). Si può notare subito che \(\cos(kx)\) è pari, mentre \(\sin(kx)\) è dispari.

Se \(f(x)\) è pari, allora lo è anche \(f(x) \cdot \cos(kx)\), mentre \(f(x) \cdot \sin(kx)\) è dispari; in tal caso si ha: \[ a_k = \frac {1} {\pi} \int_{-\pi}^{+\pi} f(x) \cos(kx) dx = \frac {1} {\pi} \int_{-\pi}^{0} f(x) \cos(kx) dx + \frac {1} {\pi} \int_{0}^{+\pi} f(x) \cos(kx) dx = \frac {2} {\pi} \int_{0}^{+\pi} f(x) \cos(kx) dx \] mentre: \[ b_k = \frac {1} {\pi} \int_{-\pi}^{+\pi} f(x) \sin(kx) dx = \frac {1} {\pi} \int_{-\pi}^{0} f(x) \sin(kx) dx + \frac {1} {\pi} \int_{0}^{+\pi} f(x) \sin(kx) dx = 0 \] Ciò accade perché, nel caso di \(a_k\), i due integrali della somma sono uguali e dello stesso segno, mentre per \(b_k\) sono uguali ma di segno opposto; si ottiene allora una serie di Fourier di soli coseni.

Analogamente, se \(f(x)\) è dispari, allora lo è anche \(f(x) \cdot \cos(kx)\), mentre \(f(x) \cdot \sin(kx)\) è pari; in tal caso si ha: \[ a_k = \frac {1} {\pi} \int_{-\pi}^{+\pi} f(x) \cos(kx) dx = \frac {1} {\pi} \int_{-\pi}^{0} f(x) \cos(kx) dx + \frac {1} {\pi} \int_{0}^{+\pi} f(x) \cos(kx) dx = 0 \] mentre: \[ b_k = \frac {1} {\pi} \int_{-\pi}^{+\pi} f(x) \sin(kx) dx = \frac {1} {\pi} \int_{-\pi}^{0} f(x) \sin(kx) dx + \frac {1} {\pi} \int_{0}^{+\pi} f(x) \sin(kx) dx = \frac {2} {\pi} \int_{0}^{+\pi} f(x) \sin(kx) dx \] Ciò accade perché, nel caso di \(a_k\), i due integrali della somma sono uguali ma di segno opposto, mentre per \(b_k\) sono uguali e dello stesso segno; si ottiene allora una serie di Fourier di soli seni.

Vediamo ora come si possono applicare questi concetti al caso che vogliamo affrontare. La Figura 19.6 mostra un'onda quadra di periodo \(T = 2\pi\) che oscilla tra i due valori costanti \(-d\) e \(+d\); per quanto visto prima, se \(T \ne 2\pi\), valgono ugualmente tutti i concetti esposti (si tratta solamente di effettuare un numero maggiore di calcoli). Nell'intervallo \([-\pi, +\pi]\) l'equazione di questa funzione è: \[ f(x) = \left\{ \begin{array}{ll} 0 & \quad x = -\pi \\ -d & \quad -\pi \lt x \lt 0 \\ 0 & \quad x = 0 \\ +d & \quad 0 \lt x \lt +\pi \\ 0 & \quad x = +\pi \end{array} \right. \] Tutti i punti \(x = k\pi\), con \(k \in \Bbb Z\), sono di discontinuità per la funzione e vi è un numero finito di essi (c'è solo lo zero) internamente all'intervallo \([-\pi, +\pi]\); in tutti gli altri punti (cioè, nei punti interni agli intervalli \([-\pi, 0]\) e \([0, +\pi]\)), la funzione è continua e monotòna, per cui risulta soddisfatto il teorema di Dirichlet.
In sostanza, la serie di Fourier converge a \(-d\) nei punti interni all'intervallo \([-\pi, 0]\) e a \(+d\) nei punti interni all'intervallo \([0, +\pi]\); nel punto di discontinuità \(0\) converge a: \[ \frac {\lim_{x \to 0^-} f(x) + \lim_{x \to 0^+} f(x)} {2} = \frac {-d + d} {2} = 0 \] Analogamente, in ciascuno degli estremi \(-\pi\) e \(+\pi\) converge a: \[ \frac {\lim_{x \to -\pi^+} f(x) + \lim_{x \to +\pi^-} f(x)} {2} = \frac {-d + d} {2} = 0 \] Da queste considerazioni si deduce che otterremo un'ottima approssimazione nei punti di continuità, mentre in prossimità dei punti \(0\), \(-\pi\) e \(+\pi\) ci saranno dei "problemi" che analizzeremo più avanti.

La nostra onda quadra è una funzione dispari \((f(-x) = - f(x))\), per cui la sua serie di Fourier sarà formata da soli seni.
Considerando il fatto che nei punti interni a \([0, +\pi]\) si ha \(f(x) = d\), risulta semplicissimo calcolare i coefficienti \(b_k\); si ha, infatti: \[ b_k = \frac {2} {\pi} \int_{0}^{+\pi} d \sin(kx) dx = -\frac {2d}{k\pi} \bigg | \cos(kx) \bigg |_{0}^{+\pi} = -\frac {2d}{k\pi} \bigg [ \cos(k\pi) - 1 \bigg ] = \frac {2d}{k\pi} \bigg [ 1 - (-1)^k \bigg ] \] Si vede subito che il termine \([1 - (-1)^k]\) vale \(2\) per \(k = 1, 3, 5, \dots\), mentre vale \(0\) per \(k = 2, 4, 6, \dots\); sono presenti quindi solo i coefficienti di posto dispari, i quali assumono i seguenti valori: \[ b_1 = \frac {4d} {\pi}, \quad b_3 = \frac {4d} {3\pi}, \quad b_5 = \frac {4d} {5\pi}, \quad b_7 = \frac {4d} {7\pi}, \quad \dots \] La serie di Fourier della nostra onda quadra è quindi la seguente: \[ f(x) = \frac {4d} {\pi} \bigg [\sin(x) \bigg] + \frac {4d} {\pi} \bigg [\frac {\sin(3x)} {3} \bigg] + \frac {4d} {\pi} \bigg [\frac {\sin(5x)} {5} \bigg] + \frac {4d} {\pi} \bigg [\frac {\sin(7x)} {7} \bigg] + \dots \] Ovviamente, non possiamo pretendere di calcolare infiniti termini, perché ciò richiederebbe un tempo infinito; dobbiamo limitarci a n termini che, sommati tra loro, ci forniranno una approssimazione dell'onda quadra, tanto migliore quanto più grande è lo stesso n.
Si può notare che, per n=1, si ottiene una semplice sinusoide di periodo \(2\pi\) e frequenza \(1 / (2\pi)\); si tratta dell'armonica fondamentale della \(f(x)\). Per n=2, la precedente sinusoide viene sommata ad una sinusoide di frequenza tripla che rappresenta la terza armonica; per n=3 abbiamo la somma tra la prima, la terza e la quinta armonica e così via.
Al crescere di n, le elaborazioni aumentano notevolmente; in pratica, i termini da calcolare sono dati da n moltiplicato per il numero di pixel che formano il grafico. Supponiamo, ad esempio, di voler associare all'intervallo \([-\pi, +\pi]\) una ampiezza pari a 600 pixel sullo schermo, con n=400 armoniche; per ciascuna delle 600 ascisse x, dobbiamo sommare 400 termini in modo da ricavare la relativa ordinata y e quindi il punto P(x,y). Ripetendo il procedimento per tutti i 600 punti, dobbiamo calcolare ben 600*400=240000 termini!

Data la pesantezza dei calcoli, non conviene visualizzare i pixel man mano che vengono ottenuti; è meglio salvare le relative coordinate in un apposito buffer e mostrare il grafico tutto in una volta.
Per rendere il programma più interessante, possiamo anche creare una animazione che mostra l'evolversi del grafico in seguito alla somma delle varie armoniche; a tale proposito, ci serve un buffer per le somme parziali delle armoniche e un buffer per i pixel dello schermo. Creiamo un loop principale da ripetere un numero di volte pari al numero desiderato di armoniche; all'interno di tale loop, ripristiniamo lo sfondo dello schermo, calcoliamo la nuova armonica (da sommare a quella precedente) e visualizziamo i pixel del grafico.
La Figura 19.7 mostra il programma FOURIER.ASM che mette in pratica tutti i concetti appena esposti. L'ampiezza in pixel dell'intervallo \([-\pi, +\pi]\) è data dalla costante XINTERVAL; tale costante è modificabile a piacere in quanto il programma è interamente parametrizzato.
Il vettore delle armoniche vharmonics è formato da XINTERVAL DWORD, tutte inizializzate a +0.0. Il vettore dei pixel vpixelgraph è formato da XINTERVAL coppie di WORD; la prima WORD contiene le ordinate dei pixel, tutte inizializzate a zero, mentre la seconda WORD contiene i colori dei pixel. Per inizializzare i colori dei pixel viene usata la procedura vesa_getpixel256 della libreria VESALIB; tale procedura restituisce il colore del pixel in posizione (x,y).

La procedura plot_squarewave visualizza sullo schermo l'onda quadra da approssimare; a tale proposito, vengono chiamate le procedure vesa_puthline256 (linee orizzontali) e vesa_putvline256 (linee verticali). Ovviamente, l'onda quadra deve avere lo stesso fattore di scala delle armoniche che la approssimano; come si può notare nel listato, il fattore di scala è proprio la costante d che rappresenta l'ampiezza di oscillazione dell'onda stessa (come sappiamo, FSCALE richiede che d debba essere un numero esprimibile come potenza intera di 2).

La procedura plot_fourier svolge il lavoro più importante. Prima di tutto, vengono inizializzati vari parametri necessari per i calcoli, come \(4 / \pi\), \(1 / k\) (con \(k\) iniziale pari a \(1\)), \(\Delta k = 2\) e \(\Delta x\); in seguito viene avviato il loop principale, da ripetere un numero di volte pari a quello delle armoniche da sommare.
All'interno del loop principale sono presenti tre loop secondari. Il primo ripristina lo sfondo dello schermo usando vesa_putpixel256; i dati dei pixel vengono letti dal vettore vpixelgraph. Il secondo loop calcola una nuova armonica, la somma a quelle precedenti (memorizzate in vharmonics), ricava le coordinate (x,y) dei nuovi pixel (con i rispettivi colori ottenuti tramite vesa_getpixel256) e salva tutto in vpixelgraph. Il terzo loop disegna il grafico con vesa_putpixel256; anche in questo caso, i dati dei pixel vengono prelevati da vpixelgraph.

Da un punto di vista matematico, si può notare che al crescere di n si ottengono armoniche che offrono un contributo sempre più piccolo al risultato finale; ciò è dovuto al fatto che \(\sin(kx)\) risulta diviso per k, il quale viene incrementato di 2 ad ogni loop.
Un altro aspetto importante, come era stato anticipato, riguarda il fatto che in prossimità dei punti di discontinuità non possiamo aspettarci una buona approssimazione dell'onda quadra; ciò accade in quanto, in ciascuno di tali punti, la serie di Fourier converge, non alla \(f(x)\), ma alla media aritmetica dei limiti sinistro e destro della \(f(x)\).
La Figura 19.8 mostra il livello di approssimazione per n=6 armoniche. Si possono già notare le maggiori oscillazioni in prossimità dei punti di discontinuità; tale comportamento prende il nome di fenomeno di Gibbs.
Per n molto grande, si ottiene un'ottima approssimazione dell'onda quadra nei suoi punti di continuità, mentre il fenomeno di Gibbs non scompare; la Figura 19.9 illustra tali aspetti per n=400. In realtà, se fosse possibile ingrandire la Figura 19.9, nei tratti dove l'onda quadra è continua si noterebbero 400 micro oscillazioni per ogni semi onda; le ampiezze di tali oscillazioni tendono a zero solo quando n diventa infinitamente grande.

19.5 Insieme frattale di Mandelbrot

Nel Capitolo 16 è stato spiegato che, nell'insieme \(\Bbb R\) dei numeri reali, non tutte le operazioni sono possibili; ad esempio, non possiamo calcolare \(\sqrt{-1}\) in quanto non esiste alcun numero reale che elevato al quadrato dia come risultato un numero negativo. Si rende necessario quindi un ampliamento di \(\Bbb R\) che porta all'insieme \(\Bbb C\) dei numeri complessi.

Si definisce numero complesso una coppia ordinata di numeri reali, rappresentata dal simbolo \((a, b)\); per i motivi che vedremo più avanti, \(a\) prende il nome di parte reale, mentre \(b\) è il coefficiente dell'immaginario. Una coppia del tipo \((a, 0)\) prende il nome di numero complesso reale, mentre nel caso della coppia \((0, b)\) si parla di numero complesso puramente immaginario.
Due numeri complessi \((a, b)\) e \((c, d)\) si dicono uguali quando \(a = c\) e \(b = d\); non ha alcun senso invece la definizione di maggiore o minore tra due numeri complessi.

La somma tra due numeri complessi è un nuovo numero complesso che ha per parte reale la somma delle parti reali e per coefficiente dell'immaginario la somma dei coefficienti dell'immaginario; si ha quindi: \[ (a, b) + (c, d) = (a + c, b + d) \] Questa regola si estende al caso di tre o più numeri complessi. Si può notare che si ha: \[ (a, b) + (0, 0) = (a, b) \quad e \quad (a, b) + (-a, -b) = (0, 0) \] Per questo motivo, \((0, 0)\) viene chiamato numero complesso zero, mentre \((-a, -b)\) si chiama opposto di \((a, b)\) e lo si indica come \(-(a, b)\).

La differenza tra due numeri complessi non è altro che la somma del primo con l'opposto del secondo; si ha quindi: \[ (a, b) - (c, d) = (a, b) + (-c, -d) = (a - c, b - d) \] Il prodotto tra due numeri complessi è dato dalla seguente regola (che verrà chiarita più avanti): \[ (a, b) \cdot (c, d) = (ac - bd, ad + bc) \] Si può notare che: \[ (a, b) \cdot (1, 0) = (a \cdot 1 - b \cdot 0, a \cdot 0 + b \cdot 1) = (a, b) \] Per questo motivo, \((1, 0)\) viene chiamato numero complesso uno.

Se \((c, d)\) è diverso da \((0, 0)\), il suo reciproco \((x, y)\) deve essere tale che: \[ (c, d) \cdot (x, y) = (cx - dy, cy + dx) = (1, 0) \] Risolvendo il sistema formato dalle due equazioni \(cx - dy = 1\) e \(cy + dx = 0\) si ottiene: \[ (x, y) = \bigg ( \frac {c} {c^2 + d^2}, \frac {-d} {c^2 + d^2} \bigg ) \] A questo punto, definiamo divisione tra due numeri complessi (di cui il secondo diverso da \((0, 0)\)), il prodotto del primo per il reciproco del secondo; si ha quindi: \[ \frac {(a, b)} {(c, d)} = \bigg ( \frac {ac + bd} {c^2 + d^2}, \frac {bc - ad} {c^2 + d^2} \bigg ) \] Si definisce unità immaginaria il numero complesso \((0, 1)\) e si pone \(i = (0, 1)\); ne consegue che: \[ i^2 = i \cdot i = (0, 1) \cdot (0, 1) = (0 - 1, 0 + 0) = (-1, 0) = -(1, 0) \] Quindi, il quadrato dell'unità immaginaria è l'opposto del numero complesso uno.

Le considerazioni appena esposte mostrano che l'insieme \(\Bbb R\) è totalmente contenuto in \(\Bbb C\); infatti, se consideriamo solo numeri complessi reali come \((a, 0)\), vediamo che essi comprendono tutto \(\Bbb R\) e, inoltre, se ripetiamo le quattro operazioni con tali numeri, otteniamo gli stessi risultati che già conosciamo in \(\Bbb R\). \[ (a, 0) + (c, 0) = (a + c, 0 + 0) = (a + c, 0) \\ (a, 0) - (c, 0) = (a - c, 0 - 0) = (a - c, 0) \\ (a, 0) \cdot (c, 0) = (ac - 0, 0 + 0) = (ac, 0) \\ \frac {(a, 0)} {(c, 0)} = \bigg ( \frac {ac + 0} {c^2 + 0}, \frac {0 - 0} {c^2 + 0} \bigg ) = \bigg ( \frac {ac} {c^2}, 0 \bigg ) = \bigg ( \frac {a} {c}, 0 \bigg ) \qquad (c \ne 0) \] Per convenzione si pone allora \((a, 0) = a\), per qualunque numero complesso reale; nel caso di \(i^2 = (-1, 0)\) si ha quindi: \[ i^2 = (-1, 0) = -1 \implies i = \sqrt {-1} \] Grazie ai concetti appena esposti, possiamo ottenere una forma dei numeri complessi molto più adatta per i calcoli; a tale proposito, osserviamo che risulta: \[ (a, b) = (a, 0) + (0, b) = (a, 0) + (0, 1) \cdot (b, 0) = (a, 0) + i \cdot (b, 0) \] Essendo \((a, 0)\) e \((b, 0)\) due numeri complessi reali che possiamo indicare come \(a\) e \(b\), si ottiene la cosiddetta forma algebrica dei numeri complessi: \[ (a, b) = (a, 0) + i \cdot (b, 0) = a + ib \] \(a\) è la parte reale, \(ib\) è la parte immaginaria, mentre \(b\) è il coefficiente dell'immaginario.

Il numero complesso \(a - ib\) prende il nome di coniugato di \(a + ib\); moltiplicando \(a + ib\) per il suo coniugato si ottiene: \[ (a + ib) \cdot (a - ib) = a^2 - iab + iab -i^2 b^2 = a^2 + b^2 \] Il numero \(a^2 + b^2\) è la norma del numero complesso \(a + ib\) ed è chiaramente un numero reale non negativo.
Il numero \(\sqrt{a^2 + b^2}\) è il modulo del numero complesso \(a + ib\) ed è anch'esso un numero reale non negativo.

Ripetendo le quattro operazioni con numeri complessi in forma algebrica e ricordando che \(i^2 = -1\), si giustificano i risultati ottenuti in precedenza: \[ (a + ib) + (c + id) = a + ib + c + id = (a + c) + i(b + d) \\ (a + ib) - (c + id) = a + ib - c - id = (a - c) + i(b - d) \\ (a + ib) \cdot (c + id) = ac + iad + ibc +i^2bd = (ac - bd) + i(ad + bc) \\ \frac {a + ib} {c + id} = \frac {a + ib} {c + id} \cdot \frac {c - id} {c - id} = \frac {ac + bd} {c^2 + d^2} + i \frac {bc - ad} {c^2 + d^2} \] Un numero complesso \(a + ib\) può avere anche una rappresentazione geometrica; la Figura 19.10 illustra tale aspetto. In sostanza, consideriamo un sistema di riferimento cartesiano ortogonale e monometrico Oxy nel piano; al numero complesso \(a + ib\) associamo il punto del piano che si ottiene prendendo a come ascissa e b come ordinata.
L'asse x prende il nome di asse reale, l'asse y è detto asse immaginario, mentre il piano così ottenuto viene definito piano di Gauss.
Applicando il teorema di Pitagora al triangolo rettangolo di lati a, b, d in Figura 19.10, si vede subito che il modulo del numero complesso \(a + ib\) non è altro che la distanza d del punto P dall'origine O.

I numeri complessi vengono largamente impiegati nella geometria dei frattali; si tratta di oggetti geometrici autosomiglianti, nel senso che mantengono sempre la stessa forma qualunque sia il livello di ingrandimento a cui vengono sottoposti.
Uno dei primi esempi di frattale risale all'inizio del 1900 ed è la famosa curva di Von Koch; il suo scopo è quello di dimostrare l'esistenza di curve prive di tangente. La Figura 19.11 mostra il procedimento che si segue per costruire questa particolare curva. Si parte da un segmento di lunghezza L, che viene diviso in 3 parti uguali, ciascuna delle quali ha quindi lunghezza L/3; la parte centrale viene cancellata e sostituita con due lati (entrambi di lunghezza L/3) di un triangolo equilatero. I 4 nuovi lati così ottenuti vengono sottoposti allo stesso procedimento; ciascun lato di lunghezza L/3 dà quindi origine a 4 lati di lunghezza L/9, per un totale di 4*4=16 nuovi lati. Ripetendo il procedimento all'infinito, si ottiene una curva priva di tangente, in quanto dotata di infinite punte; infatti, per ogni punta si hanno due semitangenti con pendenza diversa.
Qualunque sia il livello di ingrandimento a cui viene sottoposta la curva di Von Koch, si ottiene sempre la stessa forma del profilo.

Un altro esempio di frattale è legato al cosiddetto insieme di Mandelbrot, diventato famosissimo grazie alle spettacolari immagini a cui dà luogo; si tratta di un caso importante in quanto solo con l'avvento della computer grafica è stato possibile visualizzare in dettaglio la sua forma, ad opera del matematico Benoît Mandelbrot.
A dispetto della sua complessità, l'insieme di Mandelbrot viene generato attraverso una formula matematica incredibilmente semplice; infatti, si tratta dell'insieme dei numeri complessi c per i quali risulta convergente la successione: \[ z_n = z_{n-1}^2 + c \qquad n = 1, 2, 3, \dots \qquad z_0 = 0 \qquad c, z_n \in \Bbb C \] In pratica, prendiamo un numero complesso c, poniamo n=1 e calcoliamo il nuovo numero complesso \(z_1 = z_0^2 + c\), con \(z_0 = 0\); poniamo poi n=2 e calcoliamo il nuovo numero complesso \(z_2 = z_1^2 + c\). Proseguendo in questo modo, otteniamo la successione di numeri complessi \(z_1, z_2, z_3, \dots\); rappresentando tutto sul piano di Gauss, tale successione costituisce la cosiddetta orbita relativa al numero complesso c.
Al crescere di n, può capitare che l'orbita diventi sempre più "stretta", convergendo verso un punto preciso del piano di Gauss; si dice allora che la successione è convergente e il punto c appartiene all'insieme di Mandelbrot. Viceversa, può capitare che al crescere di n l'orbita diventi sempre più "larga", divergendo verso l'infinito; si dice allora che la successione è divergente e il punto c non appartiene all'insieme di Mandelbrot.

Si dimostra che la successione è divergente quando, al crescere di n il modulo di zn diventa maggiore di 2; ricordando il significato geometrico del modulo, possiamo dire che l'orbita diverge quando supera la distanza 2 dall'origine.

Poniamo, ad esempio, c=-0.5+i; si ottiene allora una successione di numeri complessi i cui primi 4 elementi sono elencati in Figura 19.12. Come si può notare, per n=4 si ottiene \(z_4 = -0,37109 + 3,125i\) e \(|z_4| = 3,14696\); il punto P(-0.5,1), associato a c=-0.5+i, non appartiene quindi all'insieme di Mandelbrot.

Poniamo ora c=0.2+0.5i; si ottiene allora una successione di numeri complessi i cui primi 8 elementi sono elencati in Figura 19.13. Se fosse possibile proseguire all'infinito, si scoprirebbe che non si ottiene mai \(|z_n| \gt 2\); il punto P(0.2,0.5), associato a c=0.2+0.5i, appartiene quindi all'insieme di Mandelbrot.

I due casi appena trattati vengono illustrati graficamente in Figura 19.14; in questo modo appare chiaro il significato di "orbita divergente" e "orbita convergente". In Figura 19.14a si vede che l'orbita relativa a c=-0.5+i si allontana sempre più dall'origine, sino a superare la distanza 2; in Figura 19.14b si vede che l'orbita relativa a c=0.2+0.5i si stringe a spirale, convergendo verso un punto che si trova a distanza inferiore a 2 dall'origine.

Le considerazioni appena esposte ci fanno capire che la scrittura di un programma per la visualizzazione dell'insieme di Mandelbrot, comporta parecchie difficoltà legate al criterio da adottare per stabilire se un punto appartiene o no all'insieme stesso; ciò è dovuto al fatto che, per sapere se \(|z_n| \gt 2\) o meno, potrebbe essere necessario calcolare un elevato numero n di termini della successione. Non possiamo certo pensare di ricorrere ad un valore di n enorme, perché ciò comporterebbe tempi di elaborazione insostenibili; basti pensare al fatto che, ad esempio, alla risoluzione di 640x480 pixel, con n=10000, se tutti i punti appartenessero all'insieme di Mandelbrot, dovremmo calcolare sino a 10000 termini della successione per ognuno dei 307200 pixel dello schermo!
Il criterio più seguito consiste nell'utilizzare come valore limite di n, il numero di colori della tavolozza (palette) che abbiamo a disposizione; nel nostro caso, stiamo sfruttando le modalità video grafiche a 256 colori, per cui possiamo porre n=256. In pratica, per ogni pixel dello schermo calcoliamo i vari termini della successione e, se otteniamo \(|z_n| \gt 2\) con \(n \lt 256\), siamo sicuri che quel punto non appartiene all'insieme di Mandelbrot e gli assegniamo proprio il colore n; se, invece, arriviamo a 256 con \(|z_n| \le 2\), allora siamo costretti a supporre (anche se ciò potrebbe essere un errore) che quel punto appartenga all'insieme di Mandelbrot e gli assegniamo un colore apposito (in genere, si usa il nero).
Un altro aspetto importante riguarda la possibilità di effettuare degli zoom su determinate parti dell'insieme di Mandelbrot; chiaramente, le parti più interessanti sono quelle relative al bordo frattale dell'insieme stesso. Per affrontare tale aspetto, dobbiamo tenere presente che la FPU ha una precisione limitata e non sarà quindi possibile spingersi oltre un certo livello di ingrandimento; possiamo alleviare questo problema sfruttando i numeri in floating point di tipo Long Real a 64 bit, che ci offrono una precisione di circa 15-16 cifre dopo la virgola.

La Figura 19.15 mostra il programma MANDELBR.ASM che mette in pratica tutti i concetti appena esposti. Molte delle impostazioni presenti in MANDELBR.ASM seguono i criteri già esposti per i precedenti esempi; l'unica differenza sta nel fatto che, per fare posto ad un apposito menu sulla destra, il grafico risulta traslato verso sinistra.
Considerando il fatto che l'insieme di Mandelbrot è totalmente contenuto all'interno di una circonferenza di raggio 2, sono stati scelti come estremi reali -2.2 e +2.2, sia lungo l'asse x, sia lungo l'asse y; come al solito, i punti reali vengono convertiti in pixel sfruttando il fatto che, ad ogni incremento di 1 pixel corrisponde un proporzionale incremento \(\Delta\) in floating point (si noti che il grafico è quadrato, per cui gli incrementi sono identici su entrambi gli assi).

Il grafico viene tracciato a partire dal punto più in basso a sinistra, di coordinate (-2.2,-2.2); alla risoluzione 640x480, tale punto convertito in pixel corrisponde a (-320,-240). Si pone quindi c=-2.2-2.2i e si effettuano le varie iterazioni per sapere se il punto appartiene o meno all'insieme; ad ogni iterazione, si calcola un nuovo \(z_n = z_{n-1}^2 + c\) tenendo conto che: \[ z_{n-1}^2 + c = (x + iy)^2 + (a + ib) = (x^2 - y^2 + a) + i(2xy + b) \] Quindi, la parte reale è \(x^2 - y^2 + a\), mentre il coefficiente dell'immaginario è \(2xy + b\).
A questo punto dobbiamo verificare se \(\sqrt{x_n^2 + y_n^2} \gt 2\); a tale proposito, si può rendere il calcolo più veloce osservando che tale comparazione equivale a \(x_n^2 + y_n^2 \gt 4\) (in questo modo risparmiamo parecchi cicli macchina dovuti all'istruzione FSQRT).
Come è stato spiegato nei precedenti capitoli, anche la comparazione tra numeri in floating point comporta l'uso del registro FLAGS/EFLAGS della CPU; supponendo di avere \(x_n^2 + y_n^2\) in ST(0) e il valore 4 in una variabile shm_cmp di tipo QWORD, possiamo scrivere: Terminate le iterazioni, visualizziamo il pixel corrispondente (-320,-240) con un colore opportuno; se il punto appartiene all'insieme (256 iterazioni), viene usato il nero, altrimenti si usa il colore corrispondente proprio al numero di iterazioni (tra 0 e 255) eseguite. Il passaggio al punto successivo lungo l'asse x consiste nell'incrementare di 1 l'ascissa del pixel e di \(\Delta\) quella del punto reale
Tutto ciò viene effettuato in un loop interno al termine del quale si passa ad un loop esterno che incrementa l'ordinata e dà inizio così all'elaborazione di una nuova riga di punti; una volta esaurite tutte le righe, il tracciamento del grafico è terminato.
La Figura 19.16 mostra l'intero insieme frattale di Mandelbrot prodotto dal programma MANDELBR.EXE; sulla destra si può notare un menu con varie scelte. La Figura 19.16 rende evidente il cerchio di raggio 2 che include tutto l'insieme. La parte esterna al cerchio è colorata di nero in quanto tutti i suoi punti si trovano sicuramente a distanza maggiore di 2 dall'origine (e quindi non appartengono all'insieme); di conseguenza, i calcoli producono 0 iterazioni e all'indice 0 corrisponde proprio il colore nero nella tavolozza standard.
Il bordo dell'insieme di Mandelbrot viene definito frattale non in senso stretto; infatti, più che il ripetersi delle stesse forme, si ha una varietà infinita di forme. I differenti colori al di fuori dell'insieme indicano la velocità con la quale la successione diverge.

Il programma permette anche di effettuare degli zoom, finché la precisione della FPU lo consente; sono possibili quasi 30 zoom, ciascuno dei quali ingrandisce l'area selezionata di 3 volte, per cui si può arrivare ad ingrandimenti pari a quasi 30*3=90 volte. Come si nota nel menu a destra, premendo il tasto [G] viene mostrata una griglia che suddivide il grafico in 9 parti e rende più semplice la scelta dell'area che si vuole ingrandire; le varie parti possono essere selezionate premendo i tasti numerici da [1] a [9], disposti nel menù in modo da corrispondere alla griglia.
Per gestire lo zoom, vengono modificate le coordinate reali, in modo che rappresentino il quadrato selezionato nella griglia; le coordinate in pixel ovviamente non vengono modificate in quanto si riferiscono all'intero schermo. A causa di queste modifiche, è necessario chiaramente ricalcolare il \(\Delta\), il quale, ad ogni zoom diventerà via via più piccolo; si può continuare a zoomare quindi finché non si va oltre la precisione dei numeri in virgola mobile di tipo Long Real a 64 bit.
La Figura 19.17 mostra un dettaglio ottenuto dopo 5 zoom, pari a 5*3=15 ingrandimenti; per tornare all'intero insieme basta premere il tasto [F]. I colori mostrati sono quelli della tavolozza standard VGA da 256 colori in formato RGB e potrebbero differire da un sistema operativo all'altro o anche da un emulatore all'altro; per evitare problemi di questo genere (ma anche per rendere il grafico più suggestivo), si può sfruttare la procedura vesa_setpalette256 della libreria VESALIB, la quale permette di installare una propria tavolozza. Per tutti i dettagli su tale procedura e sulla struttura dei colori RGB, si veda il Capitolo 11 (La memoria video in modalità grafica VESA).

Bibliografia

Intel® 64 and IA-32 Architectures Software Developer’s Manual - Volume 2 - Instruction Set Reference
(disponibile nella sezione Documentazione tecnica di supporto al corso assembly dell’ Area Downloads di questo sito - 325383-sdm-vol-2abcd.pdf)

Testi di Analisi Matematica I e II