Mandelbrot in 67 secondi su C64

Riuscire a calcolare l'insieme di Mandelbrot su un C64 e' stata una dellle sfide (perse) di quando ero ragazzo e leggevo Dewdney su Le Scienze

Adesso dopo tanti anni sono arrivato a 67 secondi  e sotto al minuto usando ottimizzazioni un po' piu' spinta (si perde pero' la cuspide del cardioide). Ulteriori miglioramenti possono essere fatti abbassando il limite massimo della variable k e usando una lookup table dei quadrati gia' calcolata

Per il calsolo sono stati utilizzati i seguenti sistemi

1) Emulazione Vice 3.3 al 100% di velocita' per emulare un C64 reale
2) Uso della tabella di lookup dei quadrati
3) Matematica a numeri interi
4) Algoritmo di calcolo che non prevede moltiplicazioni ma solo somme e sottrazioni
5) Uso di eseguibile compilato con CC65
6) Uso della simmetria

Massima ottimizzaziona impiegata con perdita della cuspide del cardioide

L'algoritmo e' quello gia' visto qui . Ottimizzazione del controllo sul cardiode da Wikipedia


da compilare con

cl65 -O -o "mande.prg" - t c64 --static-locals fractint.c

in giallo l'ottimizzazione del controllo sul cardioide

per migliorare gli algoritmi si puo' fare riferimento al seguente link

===========fractint.c===============
#include <stdio.h>
#include <stdlib.h>
#include <cc65.h>
#include <conio.h>
#include <tgi.h>


static void Mand (void)
{

    static const unsigned char Palette[2] = { TGI_COLOR_WHITE, TGI_COLOR_ORANGE };
    register unsigned s,k;
    register int zi,zr,jr,ji,a,test;
    int q[180];

    tgi_setpalette (Palette);
    tgi_setcolor (1);

    for (s=1; s<170;s++)
{
q[s] = (int) (s*s/32);
}

    for (jr=-63;jr<18;jr++)
{
for (ji=-63; ji<1;ji++)
{
zi = zr = 0;

test = q[abs(jr-8)] + q[abs(ji)];
if (test*(test+(jr-8)) >  (q[abs(ji)]/4))
{
for (k=0; k<=25;k++)
{
a = q[abs(zr)] - q[abs(zi)] + jr;
zi = q[abs(zr)]+q[abs(zi)]-q[abs(zr-zi)]+ji;
zr = a;
if (q[abs(zr)] + q[abs(zi)] >  128 )
{
tgi_setcolor(k%2);
tgi_setpixel(jr+64,ji+64);
tgi_setpixel(jr+64,128-abs(63+ji));
k=26;
}
}
}
}
}

    cgetc ();
    tgi_clear ();
}



int main (void)
{

    tgi_load_driver (tgi_stddrv);
    tgi_init ();
    tgi_clear ();

    Mand();

    tgi_unload ();
    return EXIT_SUCCESS;
}

CC65 Tgidemo su Vice x64

Il compilatore CC65 permette di scrivere programmi in modalita' grafica mediante TGI ma non ero mai riuscito a farli poi eseguire dall'emulatore/C64 perche' l'eseguibile non riusciva a trovare la libreria grafica...la soluzione e' quella di creare un disco virtuale in cui copiare la libreria virtuale ed il file PRG

Si puo' iniziare compilando la demo con il comando

../bin/cl65 -O -o "tgidemo.prg"  -t c64 tgidemo.c

A questo punto con il comando c1541 (incluso in Vice) si crea un nuovo disco e si copiano i file tgidemo.prg e c64-hi.tgi

c1541 -format diskname,id d64 my_diskimage.d64 -attach my_diskimage.d64 -write tgidemo.prg
c1541 -attach my_diskimage.d64 -write c64-hi.tgi

dall'emulatore si seleziona Autostart Image ..si vedra' un messaggio "This program needs...." si preme y e si attende qualche secondo
A questo punto apparira' la schermata dell'immagine soprastante

Mandelbrot testuale con CC65 su C64

Ero partito con l'idea ambiziosa di creare l'insieme di Mandelbrot in Assembly 6502 per Commodore 64 usando un po' l'esperienza di questo precedente post
In sintesi volevo un calcolo rapido utilizzando solo gli interi e la modalita' testo (un po' per riprendere Tequila Virus per DOS)


Alla fine ho barato....usando CC65, il crosscompilatore per 6502 e' possibile ottenere, oltre all'eseguibile. anche un listato Assembly....e tra gli esempi c'e' pure un calcolo di Mandelbrot che rispetto l'algoritmo che mi ero prefisso

Il tempo di esecuzione e' stato 1 minuto e 37 secondi

Versione in C
--------------------------------------------------------------------------------------------
include <stdlib.h>
#include <stdio.h>
#include <conio.h>

#define SCREEN_X        40
#define SCREEN_Y        22

#define maxiterations   30
#define fpshift         (10)
#define tofp(_x)        ((_x)<<fpshift)
#define fromfp(_x)      ((_x)>>fpshift)
#define fpabs(_x)       (abs(_x))

#define mulfp(_a,_b)    ((((signed long)_a)*(_b))>>fpshift)
#define divfp(_a,_b)    ((((signed long)_a)<<fpshift)/(_b))

#pragma static-locals (1);


int main (void)
{
    register unsigned char count;
    register signed short r, r1, i;
    register signed short xs, ys, xx, yy;
    register signed short x, y;
    register signed short x1,y1,x2,y2;

    x1 = -2048;
    y1 = -1024;
    x2 = 1024;
    y2 = 1024;
    

    xs = 77;
    ys = 93;

    yy = y1;
    for (y = 0; y < (SCREEN_Y); y++) {
        yy += ys;
        xx = x1;
        for (x = 0; x < (SCREEN_X); x++) {
            xx += xs;
            /* Do iterations */
            r = 0;
            i = 0;
            for (count = 0; (count < maxiterations) &&
                 (fpabs (r) < 2048) && (fpabs (i) < 2048);
                 ++count) {
                r1 = (mulfp (r, r) - mulfp (i, i)) + xx;
                i = (((signed long) r * i) >> 9) + yy;
                r = r1;
            }
    textcolor(count/2);
    printf("O");
        }
    }

    return EXIT_SUCCESS;
}
--------------------------------------------------------------------------------------------

Traduzione in Assembly
--------------------------------------------------------------------------------------------

;
; File generated by cc65 v 2.17 - Git N/A
;
.fopt compiler,"cc65 v 2.17 - Git N/A"
.setcpu "6502"
.smart on
.autoimport on
.case on
.debuginfo off
.importzp sp, sreg, regsave, regbank
.importzp tmp1, tmp2, tmp3, tmp4, ptr1, ptr2, ptr3, ptr4
.macpack longbranch
.forceimport __STARTUP__
.import _abs
.import _printf
.export _main

.segment "RODATA"

L0067:
.byte $25,$63,$00

; ---------------------------------------------------------------
; int __near__ main (void)
; ---------------------------------------------------------------

.segment "CODE"

.proc _main: near

.segment "BSS"

L0003:
.res 1,$00
L0004:
.res 2,$00
L0005:
.res 2,$00
L0006:
.res 2,$00
L0007:
.res 2,$00
L0008:
.res 2,$00
L0009:
.res 2,$00
L000A:
.res 2,$00
L000B:
.res 2,$00
L000C:
.res 2,$00
L000D:
.res 2,$00
L000E:
.res 2,$00
L000F:
.res 2,$00
L0010:
.res 2,$00

.segment "CODE"

;
; x1 = -2048;
;
ldx     #$F8
lda     #$00
sta     L000D
stx     L000D+1
;
; y1 = -1024;
;
ldx     #$FC
lda     #$00
sta     L000E
stx     L000E+1
;
; x2 = 1024;
;
ldx     #$04
lda     #$00
sta     L000F
stx     L000F+1
;
; y2 = 1024;
;
ldx     #$04
lda     #$00
sta     L0010
stx     L0010+1
;
; xs = 77;
;
ldx     #$00
lda     #$4D
sta     L0007
stx     L0007+1
;
; ys = 93;
;
ldx     #$00
lda     #$5D
sta     L0008
stx     L0008+1
;
; yy = y1;
;
lda     L000E
ldx     L000E+1
sta     L000A
stx     L000A+1
;
; for (y = 0; y < (SCREEN_Y); y++) {
;
ldx     #$00
lda     #$00
sta     L000C
stx     L000C+1
L001F: lda     L000C
ldx     L000C+1
cmp     #$16
txa
sbc     #$00
bvc     L0027
eor     #$80
L0027: asl     a
lda     #$00
ldx     #$00
rol     a
jne     L0022
jmp     L0020
;
; yy += ys;
;
L0022: lda     L0008
ldx     L0008+1
clc
adc     L000A
sta     L000A
txa
adc     L000A+1
sta     L000A+1
tax
lda     L000A
;
; xx = x1;
;
lda     L000D
ldx     L000D+1
sta     L0009
stx     L0009+1
;
; for (x = 0; x < (SCREEN_X); x++) {
;
ldx     #$00
lda     #$00
sta     L000B
stx     L000B+1
L002D: lda     L000B
ldx     L000B+1
cmp     #$28
txa
sbc     #$00
bvc     L0035
eor     #$80
L0035: asl     a
lda     #$00
ldx     #$00
rol     a
jne     L0030
jmp     L0021
;
; xx += xs;
;
L0030: lda     L0007
ldx     L0007+1
clc
adc     L0009
sta     L0009
txa
adc     L0009+1
sta     L0009+1
tax
lda     L0009
;
; r = 0;
;
ldx     #$00
lda     #$00
sta     L0004
stx     L0004+1
;
; i = 0;
;
ldx     #$00
lda     #$00
sta     L0006
stx     L0006+1
;
; for (count = 0; (count < maxiterations) &&
;
ldx     #$00
lda     #$00
sta     L0003
L003D: ldx     #$00
lda     L0003
cmp     #$1E
jsr     boolult
jeq     L0045
;
; (fpabs (r) < 2048) && (fpabs (i) < 2048);
;
lda     L0004
ldx     L0004+1
jsr     _abs
cmp     #$00
txa
sbc     #$08
bvc     L0049
eor     #$80
L0049: asl     a
lda     #$00
ldx     #$00
rol     a
jeq     L0045
lda     L0006
ldx     L0006+1
jsr     _abs
cmp     #$00
txa
sbc     #$08
bvc     L004D
eor     #$80
L004D: asl     a
lda     #$00
ldx     #$00
rol     a
jne     L0043
L0045: ldx     #$00
lda     #$00
jeq     L004E
L0043: ldx     #$00
lda     #$01
L004E: jne     L0040
jmp     L003E
;
; r1 = (mulfp (r, r) - mulfp (i, i)) + xx;
;
L0040: lda     L0004
ldx     L0004+1
jsr     axlong
jsr     pusheax
lda     L0004
ldx     L0004+1
jsr     axlong
jsr     tosmuleax
txa
ldx     sreg
ldy     sreg+1
sty     sreg
cpy     #$80
ldy     #$00
bcc     L0058
dey
L0058: sty     sreg+1
jsr     asreax2
jsr     pusheax
lda     L0006
ldx     L0006+1
jsr     axlong
jsr     pusheax
lda     L0006
ldx     L0006+1
jsr     axlong
jsr     tosmuleax
txa
ldx     sreg
ldy     sreg+1
sty     sreg
cpy     #$80
ldy     #$00
bcc     L005E
dey
L005E: sty     sreg+1
jsr     asreax2
jsr     tossubeax
jsr     pusheax
lda     L0009
ldx     L0009+1
jsr     axlong
jsr     tosaddeax
sta     L0005
stx     L0005+1
;
; i = (((signed long) r * i) >> 9) + yy;
;
lda     L0004
ldx     L0004+1
jsr     axlong
jsr     pusheax
lda     L0006
ldx     L0006+1
jsr     axlong
jsr     tosmuleax
txa
ldx     sreg
ldy     sreg+1
sty     sreg
cpy     #$80
ldy     #$00
bcc     L0063
dey
L0063: sty     sreg+1
jsr     asreax1
jsr     pusheax
lda     L000A
ldx     L000A+1
jsr     axlong
jsr     tosaddeax
sta     L0006
stx     L0006+1
;
; r = r1;
;
lda     L0005
ldx     L0005+1
sta     L0004
stx     L0004+1
;
; ++count) {
;
ldx     #$00
inc     L0003
lda     L0003
jmp     L003D
;
; printf("%c",count+50);
;
L003E: lda     #<(L0067)
ldx     #>(L0067)
jsr     pushax
ldx     #$00
lda     L0003
ldy     #$32
jsr     incaxy
jsr     pushax
ldy     #$04
jsr     _printf
;
; for (x = 0; x < (SCREEN_X); x++) {
;
lda     L000B
ldx     L000B+1
sta     regsave
stx     regsave+1
jsr     incax1
sta     L000B
stx     L000B+1
lda     regsave
ldx     regsave+1
jmp     L002D
;
; for (y = 0; y < (SCREEN_Y); y++) {
;
L0021: lda     L000C
ldx     L000C+1
sta     regsave
stx     regsave+1
jsr     incax1
sta     L000C
stx     L000C+1
lda     regsave
ldx     regsave+1
jmp     L001F
;
; return EXIT_SUCCESS;
;
L0020: ldx     #$00
lda     #$00
jmp     L0002
;
; }
;
L0002: rts

.endproc

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

Mandelbrot (per pigri) modalita' grafica con CC65

Sono riuscito ad attivare la modalita' grafica su CC65. Compilando il codice ed usando gli esempi senza modifiche avevo sempre l'errore

Error #2 initialising graphics

Cannot load driver

quando cercavo di usare TGI. Grazie al suggerimento degli sviluppatore basta flaggare a zero

#  define DYN_DRV       1

e si usa il driver in modalita' statica. La prima prova e' andata bene ma non troppo....Su VICE in modalita' velocita' 100% ho impiegato 3659.9 secondi !!!!!! (quasi quasi ritorno a Simons Basic)...non capisco ....CC65 e' compilato mentre Simons Basic e' interpretato e non c'e' un guadagno sensibile e c'e' pure il trucco della integer math contro i float si Simons...gia' la integer math...se si prova a disegnare le bande di fuga si vede che queste vengono a cuspide mentre dovrebbero essere tondeggianti e cio' deriva dalla approssimazione apportata da matematica per interi

Poi capito su questo sito http://retro64.altervista.org/blog/the-mandelbrot-set-c64-assembly-program-to-draw-it-fast-using-integer-math-part-1/ con un esempio in Assembler che genera Mandelbrot in  12 minuti ....tempo per me di attaccare la tastiera al chiodo

cl65 -O -l man.lst -vm -m man -o man.prg -t c64 man.c

man.c
--------------------------------------------------

#include <stdlib.h>
#include <conio.h>
#include <cc65.h>
#include <tgi.h>


#define SCREEN_X        320
#define SCREEN_Y        200
#define MAXCOL          2


#define maxiterations   25
#define fpshift         (14)
#define tofp(_x)        ((_x)<<fpshift) // bitwise shift left
#define fromfp(_x)      ((_x)>>fpshift)
#define fpabs(_x)       (abs(_x))

#define mulfp(_a,_b)    ((((signed long)_a)*(_b))>>fpshift)
#define divfp(_a,_b)    ((((signed long)_a)<<fpshift)/(_b))


#pragma static-locals (1);



void mandelbrot (signed short x1, signed short y1, signed short x2,
                 signed short y2)
{
    // le variabili register sono quelle che il compilatore mette in pagina zero
    // e quindi quelle a piu' rapido accesso
    register unsigned char count;
    register signed short r, r1, i;
    register signed short xs, ys, xx, yy;
    register signed short x, y;

    /* Calc stepwidth */
    xs = ((x2 - x1) / (SCREEN_X));
    ys = ((y2 - y1) / (SCREEN_Y));

    yy = y1;
    for (y = 0; y < (SCREEN_Y); y++) {
        yy += ys;
        xx = x1;
        for (x = 0; x < (SCREEN_X); x++) {
            xx += xs;
            /* Do iterations */
            r = 0;
            i = 0;
            for (count = 0; (count < maxiterations) &&
                 (fpabs (r) < tofp (2)) && (fpabs (i) < tofp (2));
                 ++count) {
                r1 = (mulfp (r, r) - mulfp (i, i)) + xx;
                /* i = (mulfp(mulfp(r,i),tofp(2)))+yy; */
                i = (((signed long) r * i) >> (fpshift - 1)) + yy;
                r = r1;
            }
            if (count == maxiterations) {
                 tgi_setcolor (0);
            } else {
                 tgi_setcolor (count%2);
          }
            tgi_setpixel (x, y);
        }
    }
}


int main (void)
{
    tgi_install (tgi_static_stddrv);
    tgi_init ();
    tgi_clear ();
    
    mandelbrot (tofp (-2), tofp (-2), tofp (2), tofp (2));
    cgetc ();
    tgi_uninstall ();
    return 0;
}

CC65 Mandelbrot (per pigri) modalita' testo

Devo ammettere che scoprire CC65 e' stato un bel passo avanti...prima avevo due opzioni sul C64: o fare i programmi in BASIC (molto molto lenti) oppure farli in Assembler. CC65 e' invece un cross compiler indirizzato al 6503 che usa la sintassi C. Il vantaggio e' che mentre il BASIC e' un linguaggio interpretato (quindi il sistema deve ogni volta effettuare la conversione in linguaggio macchina nonostante il codice sia in ogni caso tokenizzato), con CC65 si ha gia' codice macchina (anche se ovviamente non ottimizzato)

Mandelbrot modalita' testuale su C64Debugger

per esempio si puo' fare il classico Hello World e compilarlo in questo modo

cl65 -O -l hello.lst -vm -m hello -o hello.prg -t c64 hello.c

il codice generato non e' esattamente compatto (sono circa 3Kb). Inoltre e' da notare che non e' disponibile il tipo dati float

In ogni caso con poco sforzo e' possibile generare un insieme di Mandelbrot (non sono ancora riuscito a gestire la modalita' grafica via TGI)

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

#include <stdlib.h>
#include <conio.h>
#include <cc65.h>

#define SCREEN_X        40
#define SCREEN_Y        25

#define maxiterations   25
#define fpshift         (10)
#define tofp(_x)        ((_x)<<fpshift)  bitwise shift left
#define fromfp(_x)      ((_x)>>fpshift)
#define fpabs(_x)       (abs(_x))

#define mulfp(_a,_b)    ((((signed long)_a)*(_b))>>fpshift)
#define divfp(_a,_b)    ((((signed long)_a)<<fpshift)/(_b))


#pragma static-locals (1);



void mandelbrot (signed short x1, signed short y1, signed short x2,
                 signed short y2)
{
    // le variabili register sono quelle che il compilatore mette in pagina zero
    // e quindi quelle a piu' rapido accesso
    register unsigned char count;
    register signed short r, r1, i;
    register signed short xs, ys, xx, yy;
    register signed short x, y;

    /* Calc stepwidth */
    xs = ((x2 - x1) / (SCREEN_X));
    ys = ((y2 - y1) / (SCREEN_Y));

    yy = y1;
    for (y = 0; y < (SCREEN_Y); y++) {
        yy += ys;
        xx = x1;
        for (x = 0; x < (SCREEN_X); x++) {
            xx += xs;
            /* Do iterations */
            r = 0;
            i = 0;
            for (count = 0; (count < maxiterations) &&
                 (fpabs (r) < tofp (2)) && (fpabs (i) < tofp (2));
                 ++count) {
                r1 = (mulfp (r, r) - mulfp (i, i)) + xx;
                /* i = (mulfp(mulfp(r,i),tofp(2)))+yy; */
                i = (((signed long) r * i) >> (fpshift - 1)) + yy;
                r = r1;
            }
            if (count == maxiterations) {
                 cprintf("*");
            } else {
                 cprintf(" ");
          }
            /* Set pixel */
            //tgi_setpixel (x, y);
        }
    }
}

int main (void)
{
    clrscr ();
    mandelbrot (tofp (-2), tofp (-2), tofp (2), tofp (2));
    cgetc ();
    return 0;
}