mercoledì 27 aprile 2016

Mandelbrot Assembler su Linux in FPU x87 32 bit

Calcolare l'insieme di Mandelbrot in ASM e' stata per me una specie di ossessione agli inizi degli anni 90, quando i miei amici erano presi dalla creazione di demo 4K (una competizione tra programmatori per fare vedere cosa si riusciva a creare nel limite di un eseguibile da 4K)
Grazie all'aiuto di internet (il nucleo del codice di calcolo e' stato ripreso riadattandolo da questo link) sono riuscito nel mio scopo utilizzando codice x87 per i calcoli in virgola mobile ed il terminale per la visualizzazione (l'accesso diretto alla memoria video del modo 13h della VGA e' possibile solo sotto DOS)




Il codice e' piuttosto generale e permette di generare immagini di qualsiasi grandezza (basta agire sui registri di e dx) e con una finestra di visualizzazione qualsiasi (basta agire sui parametri a,b,bb,da e db) e con una profondita' di calcolo modificabile (variabile cicli)

Il codice e' totalmente commentato
Per calcolare il valore di da e di db si prende la dimensione della finestra per esempio -1.5 e +1.5, si somma (dimensione 3) e si divide per la dimensione in pixel dello schermo da quindi assume il valore di 3/80=0.0375


Esempio con finestra (-1.5,-2/1.5,2) 80x80 pixels
------------------------------------------------------
section .data
    a:        dq    -1.5     ; minimo reale della finestra
    b:        dq    -2.0     ; minimo immaginario della finestra
    bb:      dq    -2.0    ; minimo immaginario della finestra (serve per resettare la variabile)
    da:      dq    0.0375     ; incremento parte reale ((2-(-2))/80px delta_a
    db:      dq    0.05    ; incremento parte immaginaria ((1.5-1.5))/80px delta_b
    limite:  dq    4.0    ; condizione di fuga
    cicli:    dd    255     ; numero di cicli di iterazione per ogni punto
    stella:  db '*',0    ; stella (dentro insieme Mandelbrot)
    spazio: db ' ',0    ; spazio (fuori insiemeMandelbrot)
    acapo: db 0x0a,0   

section .text

global _start                   ;serve a gcc

_start:

    mov     di,80 ; dimensione dell'asse reale in pixel sullo schermo

loop_reale:

    mov    dx,80  ; dimensione dell'asse immaginario in pixel sullo schermo

    ; calcola a = a + delta_a
    fstp    st0
    fld    qword     [a]    ; mette a in st0
    fld    qword    [da]    ; mette l'incremento di a in st1
    fadd    st0,st1        ; a+da
    fstp    qword    [a]    ; mette il risultato in a

loop_immaginario:
    ; calcola b = b + delta_b
    fld    qword     [b]    ; mette b in st0
    fld    qword    [db]    ; mette l'incremento di b in st1
    fadd    st0,st1        ; b+db
    fstp    qword     [b]    ; mette il risultato in b

    ;inizia la fase di calcolo
    fstp    st0
    fld    qword    [a]     ; mette a nello stack
    fld    qword    [b]     ; mette b nello stack

    push     ecx        ; salva il registro ecx
    mov     ecx,[cicli]    ; in ecx ci sara' il contatore dei cicli di iterazioni di calcolo. Si parte da 255

    fldZ            ; mette stack st3 = 0. Sara' utilizzato come parte reale z.r = 0
    fldZ            ; mette stack st2 = 0. Sara' utilizzato come parte immaginaria z.i = 0
    fldZ            ; mette stack st1 = 0. Sara' utilizzato come parte reale quadrata z.r^2 = 0
    fldZ            ; mette stack st0 = 0. Sara' utilizzato come parte immaginaria quadrata z.i^2 = 0

l1:                ; inizia il loop di calcolo piu' interno
    fsubp    st1,st0        ; sottrae z.i^2- e lo mette nello stack come st0.
    fadd    qword     [a]    ; aggiunge a ad st0
    fstp    st3        ; prende il valore di st3 con pop
    fmulp    st1,st0        ; moltiplica st0 ed st1 z.r*z.i
    fadd    st0,st0        ; moltiplica 2*st0 facendolo come una somma st0+st0
    fadd    qword    [b]    ; aggiunge ad st0 b ovvero 2*z.r+z.i+b

    fld    st1        ; mette st1 in st0 ovvero z.r
    fmul    st0,st0        ; quadrato di z.r
    fld    st1        ; mette st1 in st0 ovvero z.i
    fmul     st0,st0        ; quadrato di z.i
    fld    st1
    fadd    st0,st1        ; somma i due quadrati z.r^2+z.i^2
    fcomp    qword     [limite]; compara st0 con il limite di calcolo 4
    fstsw    ax        ; passa il valore di st0 al registro ax
    sahf            ; mette AH nelle flags
    jae    fine        ; se ax e' sopra il valore di 4 allora esci dal ciclo di calcolo
    dec    ecx        ; decrementa il contatore della iterazioni di calcolo
    jnz    l1        ; se ecx non e' zero continua il calcolo saltando alla label l1

fine:
    ffree    st5 ; ripulisce lo stack della FPU
    ffree    st4
    ffree    st3
    ffree    st2
    ffree    st1
    ffree     st0

    ;fase di stampa sul terminale
    push edx
    mov eax,4            ; chiamata di sistema (sys_write)
    mov ebx,1            ; File descriptor 1 - standard output
    mov edx,1          ;
    cmp ecx,0        ; controlla se si e' raggiunto il limite dei cicli di iterazione
    je stampa_stella   ; se e' zero stampa uno spazio
    mov ecx,spazio     ; punta alla variabile spazio
    jmp salta
stampa_stella:
    mov ecx,stella       ; punta alla variabile stella
salta:   
    int 80h              ; Chiama il kernel
    pop edx

    pop ecx


ritorna:
    ;controlla se siamo alla fine di un loop immaginario
    dec     dx
    cmp    dx,0
    jnz    loop_immaginario

    ; resetta la variabile immaginaria ovvero la riporta al valore minimo
    fld     qword    [bb]
    fstp    qword    [b]


    ; stampa un ritorno a capo se siamo arrivati al margine della finestra
    mov eax,4            ; (sys_write)
    mov ebx,1            ; File descriptor 1 - standard output
    mov ecx,acapo        ;
    mov edx,1
    int 80h              ; Chiama il kernel


    ;controlla se siamo alla fine di un loop reale
    dec    di
    cmp    di,0
    jnz    loop_reale

    ; Ritorna al sistema operativo
    mov     eax,1
    mov    ebx,0
    int     80h

------------------------------------------------------

L'eseguibile compilato e privato del codice di debug e' lungo attorno ai 650 bytes (sicuramente si puo' fare di meglio sia come scrittura del codice che come compattezza dell'eseguibile)

E per  concludere il ricordo del primo insieme di Mandelbrot mai visualizzato (direttamente sul calcolatore di Benoit Mandelbrot)