Mandelbrot su C64 in GPascal in matematica intera

Sempre per la serie "Cosa mi sono perso negli anni 80" ho cercato se esistevano compilatori Pascal per il C64 (compilatori veri...non cross-compilatori)....io ho fatto il mio primo ed unico corso di programmazione nel 1987 su TurboPascal 5.0 su un Olivetti M24 ed usavo solo Basic su C64

Con mia grande sorpresa ho trovato diversi compilatori con caratteristiche peraltro differenti

Pascal 64 : non sono riuscito ad usarlo su VICE. Si tratta comunque di un compilatore evoluto con linker esterno. Gestisce la alta risoluzione grafica ed ha il tipo dati real. Un po' macchinoso da usare perche' si deve scrivere il file in Pascal come se si stesse editando in Basic (con i numeri di linea compresi)...si salva il file, si compila, si linka e si lancia...non proprio comodo

Oxford Pascal

GPascal : ho trovato molto comodo ed intuitivo l'editor...in 10 minuti gia' gestivo l'editing con i tasti fondamentali. Si scrive nell'editor, si compila (in PCode), si lancia ed in caso di problemi RunStop/Restore e si torna all'editor. Molto comodo. Non sono riuscito a capire come funziona il debug ed il trace. Molto intuitivo ha pero' il limite di avere solo il tipo dati numerico Integer che ha peraltro dei limiti curioso (. Si trova sia il manuale che i dischi a questo indirizzo (che dovrebbe corrispondere al sito dello sviluppatore originale)

Per mettere la prova GPascal ho provato a scrivere un programma per generare l'insieme di Mandelbrot. L'assenza di un tipo dati real mi ha costretto ad usare routine di calcolo in matematica intera...inoltre ho dovuto mantenere il moltiplicatore a 256 perche' se aumentavo la precisione di calcolo sforavo il limite degli interi che devono essere compresi tra -8388608 e + 8388607.
C'e' inoltre da segnalare che il compilatore permette di disattivare l'output video in modo da velocizzare i calcoli (con il comando graphics(6,0)

Fractint per Arduino

Questo post prende spunto da Fractint, un software ormai molto datato (e' uno dei programmi open source che da piu' tempo e' sempre sviluppato) che permetteva il calcolo "veloce" di frattali nei calcolatori che erano privi di coprocessore matematico utilizzando gli interi per calcoli in virgola mobile.

Voluto vedere se riuscivo a fare qualcosa di simile con Arduino. Per fare cio' ho utilizzato la libreria AVRFix , libreria appositamente studiata per processori ATMEL

Il programma senza nessun tipo di ottimizzazione ha mostrato un tempo di esecuzione pari a 7581 ms. A parita' di condizioni di schermo la libreria AVRFix con il tipo di variabile short ha quasi dimezzato il tempo di calcolo per un totale di 4434 ms

Ovviamente questa non e' una soluzione generale. Basta infatti cambiare il tipo di varibiale lasciando intatto il codice per usa AVRFix che le prestazioni decadono rapidamente a 10590 ms per il tipo di varibile l a 18464 ms per il tipo di variabile ll

No AVRFix (7581 ms)
-----------------------------------------
#include <SPI.h>
#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>

#define OLED_RESET 4
Adafruit_SSD1306 display(OLED_RESET);



#define SCREEN_WIDTH 32
#define SCREEN_HEIGHT 128

#if (SSD1306_LCDHEIGHT != 32)
#error("Height incorrect, please fix Adafruit_SSD1306.h!");
#endif

//float re_min = -2.019531;
float a = -2.003906;
float im_min = -1.275;
float re_factor = 0.019531;
float im_factor = 0.075;

unsigned long tempo;

float b;
float x,y,x_new,y_new;

int k,j,i;

void setup()   {
  Serial.begin(115200);
  display.begin(SSD1306_SWITCHCAPVCC, 0x3C);  
  display.display();
  delay(100);
  display.clearDisplay();

  tempo = millis();
  
  for (i=0;i<SCREEN_HEIGHT;i++)
     {
     a = a + re_factor;
     b = im_min;
     for (j=0;j<SCREEN_WIDTH;j++)
      {
        b=b+im_factor;
        x = 0;
        y = 0;
      
        for (k=0;k<=64;k++)
            {
            x_new = (x*x)-(y*y)+a;
            y_new = (2*x*y)+b;
            if (((x_new*x_new)+(y_new*y_new))>4)
                {
                display.drawPixel(i, j, WHITE);
                //
                break;
                }
       x = x_new;
       y = y_new;
       } 
      }
     }
Serial.println(millis()-tempo);
display.display();
}


void loop() {
}
-----------------------------------------

AVRFix (tipo short 4728 ms)
-----------------------------------------
#include <SPI.h>
#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>

#define OLED_RESET 4

Adafruit_SSD1306 display(OLED_RESET);
extern "C"{
  #include <avrfix.h>
}

#define WIDTH 32
#define HEIGHT 128
lfix_t Im = ftosk(-1.275);
lfix_t DRe = ftosk(0.02);
lfix_t DIm = ftosk(0.075);
lfix_t A = ftosk(-2.0);
lfix_t B = ftosk(0.0);
lfix_t DUE = ftosk(2.0);
lfix_t X,Y;
lfix_t XN,YN,YNN;
lfix_t tt;
int j,k,i;
float test;

unsigned long tempo;

void setup() {
  Serial.begin(115200);
  display.begin(SSD1306_SWITCHCAPVCC, 0x3C);  
  display.display();
  delay(100);
  display.clearDisplay();
  tempo = millis();
  
  for (i=0;i<HEIGHT;i++)  
    {
      A=A+DRe;
      B=Im;
      for (j=0;j<WIDTH;j++)
        {
          B=B+DIm;
          X = ftosk(0.0); 
          Y = ftosk(0.0);
          for (k=0;k<=64;k++)
              {
               XN = smulsk(X, X)-smulsk(Y, Y)+A;
               YN = smulsk(X, Y);
               YNN = smulsk(X, Y);
               YN = smulsk(DUE, YNN) + B;
               tt = smulsk(XN,XN)+smulsk(YN,YN);
               test = sktof(tt);
               if (test > 4.0)
                    {
                        display.drawPixel(i, j, WHITE);
                        break;
                    }
               X=XN;
               Y=YN;
              }
        }
    }
  

  Serial.println(millis()-tempo);
  display.display();
  
}

void loop() {
}
-----------------------------------------

AVRFix (tipo l 10590 ms)

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

#include <SPI.h>
#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>

#define OLED_RESET 4

Adafruit_SSD1306 display(OLED_RESET);
extern "C"{
  #include <avrfix.h>
}

#define WIDTH 32
#define HEIGHT 128
//lfix_t Re = ftolk(-2.00390625); 
lfix_t Im = ftosk(-1.275);
lfix_t DRe = ftok(0.019531);
lfix_t DIm = ftok(0.075);
lfix_t A = ftok(-2.00390625);
lfix_t B = ftok(0.0);
lfix_t DUE = ftok(2.0);
lfix_t X,Y;
lfix_t XN,YN,YNN;
lfix_t tt;
int j,k,i;
double test;

unsigned long tempo;

void setup() {
  Serial.begin(115200);
  display.begin(SSD1306_SWITCHCAPVCC, 0x3C);  
  display.display();
  delay(100);
  display.clearDisplay();
  tempo = millis();
  
  for (i=0;i<HEIGHT;i++)  
    {
      A=A+DRe;
      B=Im;
      for (j=0;j<WIDTH;j++)
        {
          B=B+DIm;
          X = ftok(0.0); 
          Y = ftok(0.0);
          for (k=0;k<=64;k++)
              {
               XN = mulk(X, X)-mulk(Y, Y)+A;
               YN = mulk(X, Y);
               YNN = mulk(X, Y);
               YN = mulk(DUE, YNN) + B;
               tt = mulk(XN,XN)+mulk(YN,YN);
               test = ktod(tt);
               if (test > 4.0)
                    {
                        display.drawPixel(i, j, WHITE);
                        break;
                    }
               X=XN;
               Y=YN;
              }
        }
    }
  

  Serial.println(millis()-tempo);
  display.display();
  
}

void loop() {

}
-----------------------------------------

AVRFix (tipo ll 18464 ms)

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

#include <SPI.h>
#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>

#define OLED_RESET 4

Adafruit_SSD1306 display(OLED_RESET);
extern "C"{
  #include <avrfix.h>
}

#define WIDTH 32
#define HEIGHT 128
//lfix_t Re = ftolk(-2.00390625); 
lfix_t Im = ftolk(-1.275);
lfix_t DRe = ftolk(0.019531);
lfix_t DIm = ftolk(0.075);
lfix_t A = ftolk(-2.00390625);
lfix_t B = ftolk(0.0);
lfix_t DUE = ftolk(2.0);
lfix_t X,Y;
lfix_t XN,YN,YNN;
lfix_t tt;
int j,k,i;
double test;

unsigned long tempo;

void setup() {
  Serial.begin(115200);
  display.begin(SSD1306_SWITCHCAPVCC, 0x3C);  
  display.display();
  delay(100);
  display.clearDisplay();
  tempo = millis();
  
  for (i=0;i<HEIGHT;i++)  
    {
      A=A+DRe;
      B=Im;
      for (j=0;j<WIDTH;j++)
        {
          B=B+DIm;
          X = ftolk(0.0); 
          Y = ftolk(0.0);
          for (k=0;k<=64;k++)
              {
               XN = lmullk(X, X)-lmullk(Y, Y)+A;
               YN = lmullk(X, Y);
               YNN = lmullk(X, Y);
               YN = lmullk(DUE, YNN) + B;
               tt = lmullk(XN,XN)+lmullk(YN,YN);
               test = lktod(tt);
               if (test > 4.0)
                    {
                        display.drawPixel(i, j, WHITE);
                        
                        break;
                    }
               X=XN;
               Y=YN;
              }
        }
    }
  

  Serial.println(millis()-tempo);
  display.display();
  
}

void loop() {

}
-----------------------------------------


Disclaimer: il titolo di questo post non e' da intendersi come un plagio verso il nome del piu' noto software Fractint ma piuttosto come un tributo ad un programma che usavo sui 386 nei primi anni 90

Fixed point math e Mandelbrot

Un esempio di utilizzo della matematica a virgola fissa, un metodo di rappresentazione dei numeri reali con un numero prefissato di cifre decimali utilizzato in particolar modo per calcolatori privi di coprocessore matematico (per esempio il software Fractint usava questo tipo di approccio per il calcolo di frattali su processori 386 e 486 all'inizio degli anni 90)

Per effettuare un calcolo dell'insieme di Mandelbrot si puo' usare attualmente la libreria FIXEDPTC (Fixed Point Math Library for C)  (rispetto ad altre librerie, come FPMLib o LibFoixMath , e' possibile settare il proprio formato dati secondo le proprie necessita')

Visto che il calcolo e' tutto compreso per valori inferiori a 8 puo' essere utile, per massimizzare la precisione, usare una rappresentazione 4.28 (4 digit per la parte intera e 28 per la parte frazionaria) su un calcolatore a 32 bit. In questo modo il piu' piccolo numero calcolabile e' 0.00000000372529. Per la rappresentazione grafica e' stata usata la libreria PngWriter gia' usata qui
L'uso della libreria e' abbastanza autoesplicativo...l'unico aspetto un po' confuso e' che per effettuare dei confronti in un ciclo if...then si deve convertire il dato da formato a virgola fissa in float.

mand.c
------------------------------------------------------
#include <pngwriter.h>

#define FIXEDPT_WBITS 4
#define WIDTH 640
#define    HEIGHT 480

#include "fixedptc.h"

int main() {
int     j,k,i;
float    test;
fixedpt Re,Im,DRe,DIm;
fixedpt X,Y,DUE;
fixedpt XN,YN,YNN;
fixedpt    A,B;

pngwriter png(WIDTH,HEIGHT,0,"mandelbrot.png");

Re = fixedpt_rconst(-2.00390625); //finestra reale tra -2 e +0.5
Im = fixedpt_rconst(-1.205); //finestra immaginaria tra -1.2 e 1.2
DRe = fixedpt_rconst(0.00390625); //2.5/640
DIm = fixedpt_rconst(0.005); // 2.4/480
DUE = fixedpt_rconst(2.0);

A = Re;

for (j=0;j<WIDTH;j++)
    {
    A = fixedpt_add(A,DRe);
    B = Im;
    for (k=0;k<HEIGHT;k++)
    {
    B = fixedpt_add(B,DIm);

    X = fixedpt_rconst(0.0);
    Y = fixedpt_rconst(0.0);

    for (i=0;i<=255;i++)
        {
        XN=fixedpt_sub(fixedpt_mul(X,X),fixedpt_mul(Y,Y))+A; // (x*x) - (y*y) + A
        YN=fixedpt_mul(X,Y); // x*y
        YNN=fixedpt_mul(DUE,YN); // 2*x*y
        YN=YNN + B; // 2*x*y*+B
        test = fixedpt_tofloat(fixedpt_mul(XN,XN) + fixedpt_mul(YN,YN)); //(XN*XN)+(YN*YN)
        if (test > 4.0)
            {
            png.plot(j,k,1.0,1.0,1.0);
            break;
            }
        X = XN;
        Y = YN;
        }
    }
    }
png.close();

}
------------------------------------------------------

Per compilare il sorgente si puo' usare i seguenti comandi
------------------------------------------------------
rm mandelbrot.png
g++ mand.c -o mand  -I/usr/local/include  -L/usr/local/lib -lpng -lpngwriter
./mand
display mandelbrot.png
------------------------------------------------------