TI-89 Mandelbrot

 Ho scoperto per puro caso che le calcolatrici Texas Instruments TI-89 (e parenti) hanno al proprio interno un processore della classe 68000 e che possono essere programmate, oltre che in TI-Basic, anche in C...proviamo un po'

Il compilatore C che ho usato si trova a questo link https://github.com/debrouxl/gcc4ti

esistono versioni meno aggiornate ma non sono cosi' complete

Per compilare si va in /trunk/tigcc-linux/scripts e si lancia ./Install

si impostano poi le variabili di ambiente

export TIGCC=/usr/local/share/gcc4ti/

export PATH=$PATH:$TIGCC/bin

il codice di questo post e' ripreso dal precedente post 

le differenze sono :

1) la funzione main in TIGCC si chiama _main (underscore main)

2) sono necessari gli include di stdio e kbd anche se di fatto non sono richiamati dal codice

#define USE_TI89 // produce all types of files

//#define USE_TI92PLUS

//#define USE_V200

#include <stdio.h> // standard ANSI C input/output support

#include <kbd.h>

#include <tigcclib.h>

#define FIXEDPT_WBITS 4

#define WIDTH 160

#define HEIGHT 100

#include "fixedptc.h"

void _main(void) {

ClrScr();

int j,k,i;

float test;

fixedpt Re,Im,DRe,DIm;

fixedpt X,Y,DUE;

fixedpt XN,YN,YNN;

fixedpt A,B;

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.015625); //2.5/160

DIm = fixedpt_rconst(0.024); // 2.4/100

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<=127;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);

if (i%2) DrawPix(j,k,A_NORMAL);

break;

}

X = XN;

Y = YN;

}

}

}

}

tigcc -O2 -o timand timand.c 

il file binario avra' una estensione .89z (nel caso si compile di per TI-92 sara' .92z)

Per il trasferimento dell'eseguibile tramite cavo USB ho usato il programma TILP funzionante su Linux impostando Direct Link

Una volta trasferito il codice sul dispositivo si puo' eseguire digitando il nome del fie per esempio se il nome del file e' timand come se fosse una funzione ...per esempio timand()

In alcuni casi la calcolatrice puo' non mostrare la linea di comando ma le icone..in questo caso si preme il tasto MODE si scrolla in basso fino all'opzione APP Desktop e se seleziona OFF

Se non si vuole usare un dispositivo fisico si puo' usare l'emulatore TIEmu. In questo caso per caricare l'eseguibile e' sufficiente premere il tasto F10 e selezionare il flie .89z

Cuda Mandelbrot

 

 

 

#include <stdio.h>
//questa libreria crea png senza altre dipendenze
#include "libattopng.h"

#define W  5000 //width image
#define H  5000 //height image
#define TX 32 // number of threads per block along x-axis
#define TY 32 // number of threads per block along y-axis

__global__
void Kernel(int *d_out, int w, int h, int itera)
{
    //nel kernel viene lanciato un thread per ogni pixel dell'immagine
    //per calcolare la riga e la colonna si usa la formula sottostante
    int c = blockIdx.x*blockDim.x + threadIdx.x; // colonna
    int r = blockIdx.y*blockDim.y + threadIdx.y; // riga
    int i = r*w + c; // posizione del pixel nell'array lineare
    if ((c >= w) || (r >= h)) return;

    // inizia il ca
    const float x_min = -0.727;
    const float x_max = -0.730;
    const float y_min = -0.247;
    const float y_max = -0.250;

    float dx,dy;

    dx = (x_max-x_min)/w;
    dy = (y_max-y_min)/h;


    float k = x_min+(r*dx);
    float j = y_min+(c*dy);
    d_out[i] = 0;

    double x,y,x_new, y_new = 0.0;
    for (int f=0; f<itera;f++)
            {
                    x_new = (x*x) - (y*y) + k;
                    y_new = (2*x*y) +j;
                    if (((x_new*x_new)+(y_new*y_new))>4)
                    {
                        d_out[i] = f;
                        return;
                    }
                    x = x_new;
                    y = y_new;
            }
 }

int main()
{
 // alloca e mette a zero (malloc alloca e basta) un array di int di dimensione u
 // uguale all'immagine desiderata sulla ram della CPU
 int *out = (int*)calloc(W*H, sizeof(int));

 // alloca un array nella GPU delle stesse dimensioni
 // la GPU viene definita come device mentre la CPU e' l'host
 int *d_out; // pointer for device array
 cudaMalloc(&d_out, W*H*sizeof(int));

 // definisce i threads
 const dim3 blockSize(TX, TY);
 const int bx = (W + TX - 1)/TX;
 const int by = (W + TY - 1)/TY;
 const dim3 gridSize = dim3(bx, by);
 //lancia il kernel sulla GPU passando come parametro l'array sulla GPU
 Kernel<<<gridSize, blockSize>>>(d_out, W, H,4096);

 //trasferisce il contenuto dell'array dalla GPU alla RAM della CPU
 cudaMemcpy(out, d_out, W*H*sizeof(int), cudaMemcpyDeviceToHost);

 //crea una PNG partendo dall'array
 libattopng_t* png = libattopng_new(W, H, PNG_GRAYSCALE);
 
 int x, y;
 for (y = 0; y < W; y++) {
   for (x = 0; x < H; x++) {
     libattopng_set_pixel(png, y, x, out[x+y*H]%255);
   }
 }
 libattopng_save(png, "mandelbrot.png");
 libattopng_destroy(png);


 cudaFree(d_out);
 free(out);
 return 0;
 }
 

 

Webassembly Golang Mandelbrot

Era da tempo che volevo provare WebAssembly (abbreviato wasm).. ho scoperto che sia GoLang che Qt permettono di compilare in Wasm e cosi' ho fatto una prova

Per eseguire applicazioni WAsm si deve configurare il webserver a riconoscere il Mime Type wasm . In Apache e' sufficiente aggiungere la linea in /etc/mime.types mentre in NGinx /etc/nginx/mime.types con application/wasm  wasm;

Il file compilato deve avere l'estensione .wasm

Per compilare il file go.wasm si usa la linea di comando

GOOS=js GOARCH=wasm go build -o wasm.wasm

ATTENZIONE: i file wasm rimangono in cache del browser. In fase di sviluppo e' quindi necessario sincerarsi di pulire la cache con CTRL+SHIFT+R (od usare una sessione anonima)

Per creare il file wasm_exec.js e' sufficiente (attenzione al punto finale)

cp "$(go env GOROOT)/misc/wasm/wasm_exec.js" .

Il programma genera una immagine dell'insieme di Mandelbrot che viene codificata in PNG e successivamente in Base64 per divenire la src del tag img della pagina web

https://github.com/c1p81/gowasm_mand

wasm.go

package main

// cp "$(go env GOROOT)/misc/wasm/wasm_exec.js" .

import (

"bytes"

"encoding/base64"

"fmt"

"image"

"image/color"

"image/png"

"io/ioutil"

"syscall/js"

)

func funzione() js.Func {

jsonFunc := js.FuncOf(func(this js.Value, args []js.Value) interface{} {

//////////////////////////////////////////////////////////////////////////////

var SCREEN_WIDTH int = 400

var SCREEN_HEIGHT int = 400

var re_min float32 = -2.0

var im_min float32 = -1.2

var re_max float32 = 1.0

var im_max float32 = 1.2

var iterazioni int = 1024

var a, b float32

var x, y, x_new, y_new, somma float32

var k, i, j int

var re_factor float32 = (re_max - re_min)

var im_factor float32 = (im_max - im_min)

m := image.NewRGBA(image.Rect(0, 0, SCREEN_WIDTH, SCREEN_HEIGHT))

for i = 0; i < SCREEN_HEIGHT; i++ {

for j = 0; j < SCREEN_WIDTH; j++ {

a = re_min + (float32(j) * re_factor / float32(SCREEN_WIDTH))

b = im_min + (float32(i) * im_factor / float32(SCREEN_HEIGHT))

x = 0

y = 0

for k = 0; k < iterazioni; k++ {

x_new = (float32(x) * float32(x)) - (float32(y) * float32(y)) + float32(a)

y_new = (float32(2) * float32(x) * float32(y)) + float32(b)

somma = (x_new * x_new) + (y_new * y_new)

if somma > 4 {

if k%2 == 0 {

m.Set(j, i, color.RGBA{0, 0, 0, 255})

} else {

m.Set(j, i, color.RGBA{255, 255, 255, 255})

}

break

}

x = x_new

y = y_new

}

}

}

buf := new(bytes.Buffer)

if err := png.Encode(buf, m); err != nil {

return ("unable to encode png")

}

readBuf, _ := ioutil.ReadAll(buf)

enc := base64.StdEncoding.EncodeToString([]byte(readBuf))

return "data:image/png;base64," + enc

})

return jsonFunc

}

func main() {

fmt.Println("Inizio")

js.Global().Set("go_function", funzione())

fmt.Println(funzione())

<-make(chan bool)

}

 

Index.html

<html>

<head>

<meta charset="utf-8"/>

<script src="wasm_exec.js"></script>

<script>

const go = new Go();

WebAssembly.instantiateStreaming(fetch("wasm.wasm"), go.importObject).then((result) => {

go.run(result.instance);

});

</script>

</head>

<body>

<center> 2

<img id="immagine" name="immagine" src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAUAAAAFCAYAAACNbyblAAAAHElEQVQI12P4

//8/w38GIAXDIBKE0DHxgljNBAAO9TXL0Y4OHwAAAABJRU5ErkJggg==" alt="Red dot" /><br><br>

<input id="button" type="submit" name="button" onclick="ff()"/> <br>

<textarea id="imgbase64" name="imgbase64" cols="80" rows="200"></textarea>

</center>

<script>

var ff = function() {

immagine.src = go_function(ff)

imgbase64.value = go_function(ff)

}

</script>

</body>

</html>

Mandelbrot M1 Metal Python

 Dopo il test di Cuda e' arrivato il momento di Metal, il linguaggio per la GPU di Apple

Gli ersempi in C++ sono un po' complicati per una semplice prova ed ho optato per la versione in Python in cui il kernel e' comunque scritto in C++ MetalCompute (negli esempi della libreria c'e' gia' una implementazione di Mandelbrot, l'ho scritta da zero per imparare...sicuramente la mia versione e' meno efficiente)

La logica seguita e' quella di Cuda, ogni thread si occupa di una riga

Si puo' solo passare un solo parametro come array di float al kernel e si ha in uscita un solo array di float. L'array di input e' generato con NumPy con le coordinate X,Y dei punti dell'insieme 

import metalcompute as mc

from PIL import Image as im

mc.init()

code = """

#include <metal_stdlib>

using namespace metal;

kernel void mandelbrot(

const device float* arr [[ buffer(0) ]],

device float* out [[ buffer(1) ]],

uint id [[ thread_position_in_grid ]]) {

float a = arr[id];

float b;

float x_new,y_new,x,y;

int iterazioni = 255;

for (int s=0; s<640;s++)

{

x=0.0;

y=0.0;

b = arr[640+s];

for (int k=0;k<iterazioni;k++)

{

x_new = (x*x)-(y*y)+a;

y_new = (2.0*x*y)+b;

if (((x_new*x_new)+(y_new*y_new))>4)

{

out[id*640+s] = (float)(k%2)*255;

k=iterazioni;

}

x = x_new;

y = y_new;

}

}

}

"""

mc.compile(code, "mandelbrot")

import numpy as np

from time import time as now

dimensioni = 640

x_start = -2.0

x_stop = 0.75

x_res = (x_stop-x_start)/dimensioni

x = np.arange(start=x_start, stop=x_stop,step = x_res,dtype='f')

y_start = -1.5

y_stop = 1.5

y_res = (y_stop-y_start)/dimensioni

y = np.arange(start=y_start, stop=y_stop,step = y_res,dtype='f')

arr = np.concatenate((x,y))

print(arr.shape)

out_buf = np.empty(dimensioni*dimensioni,dtype='f')

start = now()

mc.run(arr, out_buf, dimensioni)

end = now()

immagine = out_buf.reshape(dimensioni,dimensioni)

print(end-start)

immagine = np.int_(immagine)

np.savetxt("dati.txt",immagine.astype(int),fmt='%i')

img = im.fromarray(np.uint8(immagine) , 'L')

img.save('m1.png')

img.show()

Mandelbrot Cuda su Nvidia Jetson Nano

Ho ripreso in mano la scheda Jetson Nano per iniziare a programmare in CUDA 

Gli esempi di CUDA hanno gia' un sorgente per l'insieme di Mandelbrot ma visto che era a scopo didattico sono partito da zero

L'idea e' quella di usare la GPU per i calcoli. Ogni thread calcola una riga dell'immagine. Visto che ogni block della Nvidia puo' gestire al massimo 1024 thread la massima dimensione che il sorgente puo' generare e' 1024x1024

Per condividere una memoria tra la GPU e la CPU (ognuna puo' accere solo alle proprie risorse) si usa  cudaMallocManaged. CudaMemset port i valori della matrice tutti a zero

la funzione kernel e'e quella che viene eseguita dalla GPU. Una volta lanciato il kernel la GPU restituisce il controllo alla CPU solo quando tutti i threads sono terminati. La variabile threadIdx viene utilizzate come indice di colonna dell'immagine

al termine l'array viene convertito in una immagine ppm a scala di grigio (nello specifico si sono solo due colori)

per compilare si usa nvcc con gli switch per indicano la compilazione sull'architettura 50 (relativa alla GPU sulla Jetson) con l'ottimizzazione fast_math

nvcc -use_fast_math -arch=sm_50 mand.cu -o mand

Ho usato VSCode con il terminale ma per programmare Cuda si puo' usare NSight (una versione di Eclipse con plugin NVidia) contenuto nei Cuda Tools

questo il codice

#include <stdio.h>

#include <time.h>

__global__ void kernel (int max, char *matrice)

{

int id = threadIdx.x;

/*zoom1

const double CxMin=-0.3041;

const double CxMax=-0.1874;

const double CyMin=-0.8867;

const double CyMax=-0.7699;*/

/*zoom2*/

const double CxMin=-0.2;

const double CxMax=-0.074;

const double CyMin=-1.058;

const double CyMax=-0.933;

/*insieme completo

const double CxMin=-2.5;

const double CxMax=1.5;

const double CyMin=-2.0;

const double CyMax=2.0;*/

const int iterazioni = 4096;

double x_new,y_new,x,y;

double a,b;

double PixelWidth=(CxMax-CxMin)/max;

double PixelHeight=(CyMax-CyMin)/max;

for (int s=0;s<max;s++)

{

x = 0;

y = 0;

a = CxMin + (PixelWidth*id);

b = CyMin + (PixelHeight*s);

for (int k=0;k<iterazioni;k++)

{

x_new = (x*x)-(y*y)+a;

y_new = (2*x*y)+b;

if (((x_new*x_new)+(y_new*y_new))>4)

{

// colora il punto

matrice[id*max+s] =(k%2)*255;

k=iterazioni;

}

x = x_new;

y = y_new;

}

}

}

int main(void)

{

clock_t tic = clock();

// crea la matrice in GPU

int dimensione = 1000;

char *matcuda;

// crea un array Cuda

cudaMallocManaged(&matcuda,dimensione*dimensione*sizeof(char));

// azzera tutti i valori dell'array

cudaMemset(matcuda,0,dimensione*dimensione*sizeof(char));

// 1024 e' il numero di threads che vengono eseguiti in contemporanea

// e' il valore massimo per ogni block della NVidia

// tra parentesi si possono passare variabili al kernel

kernel <<<1, dimensione>>>(dimensione,matcuda);

cudaDeviceSynchronize();

clock_t toc = clock();

/*

for (int s=0;s<dimensione;s++)

{

for (int t=0;t<dimensione;t++)

{

printf("%i",matcuda[(dimensione*s)+t]);

}

printf("\n");

}

*/

printf("Elapsed: %f seconds\n", (double)(toc - tic) / CLOCKS_PER_SEC);

FILE *fpi = fopen("mand.ppm", "wb");

fprintf(fpi, "P5\n%d %d\n255\n", dimensione, dimensione);

fwrite(matcuda, 1, dimensione*dimensione, fpi);

fclose(fpi);

cudaFree(matcuda);

return 0;

}

QML Mandelbrot in UBPorts / Ubuntu Touch (2)

Avevo gia' provato qui a creare l'insieme di Mandelbrot su UBPorts in QML ma il programma era scandalosamente lento ed ho provato a migliorare

In questo codice il risultato di calcolo per ogni pixel viene salvato in un array imgData che poi 

trasferito in un colpo sul Canvas. In questo modo il tempo di calcolo e' divenuto inferiore ad un terzo di quello della precedente prova

==============================================================
/*
 * Copyright (C) 2020  Luca Innocenit
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; version 3.
 *
 * ubuntu-calculator-app is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

import QtQuick 2.7
import Ubuntu.Components 1.3
//import QtQuick.Controls 2.2
import QtQuick.Layouts 1.3
import Qt.labs.settings 1.0
import QtQuick.Window 2.2

MainView {
    id: root
    objectName: 'mainView'
    applicationName: 'qtcpp.com.luca.innocenti'
    automaticOrientation: true

    visible: true
      width: Screen.width
      height: Screen.width
      Canvas {
          id: mycanvas
          anchors.centerIn: parent
          width: Screen.width
          height: Screen.width
          onPaint: {
              console.log("width "+ width.toString())
              console.log("height "+ height.toString())

              var ctx = getContext("2d");
              ctx.fillStyle = Qt.rgba(0, 0, 0, 1);
              ctx.fillRect(0, 0, width, height);
              var imgData = ctx.createImageData(height, width);

              ////////
              var re_min = -2.0;
             var im_min = -1.2;
             var re_max = 0.7;
             var im_max = 1.2;
             var iterazioni = 50;

             var r;
             var a,b;
             var x,y,x_new,y_new;
             var test;
             var re_factor = (re_max-re_min);
             var im_factor = (im_max-im_min);
             var indice = 0;

              for (var i=0;i<width;i++)
                          {
                          for (var j=0;j<height;j++)
                           {
                               a = re_min+(j*re_factor/height);
                               b = im_min+(i*im_factor/width);
                               x = 0;
                               y = 0;
                               test = 0;
                              imgData.data[indice+3]= 255;
                              //console.log(indice)

                               for (var k=0;k<iterazioni;k++)
                                    {
                                    x_new = (x*x)-(y*y)+a;
                                    y_new = (2*x*y)+b;
                                    if (((x_new*x_new)+(y_new*y_new))>4)
                                          {
                                            if (k%2)
                                                {
                                                  // colora il punto
                                                  imgData.data[indice]=255;
                                                  imgData.data[indice+1]=255;
                                                  imgData.data[indice+2]=255;  
                                                }

                                            break;
                                          }
                                    x = x_new;
                                    y = y_new;
                                    }
                               indice = indice + 4;
                           }
              }
              ///////
              console.log("terminato")
              ctx.drawImage(imgData, 0, 0);

          }
      }
}

QML Mandelbrot in UBPorts / Ubuntu Touch

Ovviamente, visto che mi sono messo a studiare UBPorts, non poteva mancare un porting di Mandelbrot. La prima scelta e' stata QML che in pratica e' un sottoinsieme di HTML5 per quanto riguarda il codice

Il progetto completo in Clickable si trova a questo link
https://github.com/c1p81/lucainnoc-gmail.com

Piccola nota: il progetto e' dannatamente lento

Le icone nel progetto sono in formato .svg (Gimp non supporta direttamente questo formato)

===================================================================

/*

 * Copyright (C) 2020  Your FullName

 *

 * This program is free software: you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; version 3.

 *

 * ubuntu-calculator-app is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with this program.  If not, see <http://www.gnu.org/licenses/>.

 */

import QtQuick 2.7

import Ubuntu.Components 1.3

//import QtQuick.Controls 2.2

import QtQuick.Layouts 1.3

import Qt.labs.settings 1.0

import QtQuick.Window 2.2

MainView {

    id: root

    objectName: 'mainView'

    applicationName: 'pixel.luca.innocenti'

    automaticOrientation: true

    width: Screen.width

    height: Screen.height

    Page {

        anchors.fill: parent

        header: PageHeader {

            id: header

            title: i18n.tr('Mandelbrot QML')

        }

Canvas {

id: canvas

width:  Screen.width

height:  Screen.height

  onPaint: {

      var context=canvas.getContext("2d")

           context.fillStyle="rgb(0,0,0)";

           context.fillRect(0, 0, width, height);

           var re_min = -2.0;

           var im_min = -1.2;

           var re_max = 0.7;

           var im_max = 1.2;

           var iterazioni = 25;

           var colori = ["#FFFFFF","#FF0000","#00FF00","#0000FF","#FF00FF","#FFFF00","#00FFFF","#FFaa00","#abcedf","#fedc0a","#ab16dd","#d00fba","#edabcc","#ddacff"];

           var r;

           var a,b;

           var x,y,x_new,y_new;

           var test;

           var k,j,i;

           var re_factor = (re_max-re_min);

           var im_factor = (im_max-im_min);

           for (var i=0;i<width;i++)

            {

            for (var j=0;j<height;j++)

             {

                 a = re_min+(j*re_factor/height);

                 b = im_min+(i*im_factor/width);

                 x = 0;

                 y = 0;

                 test = 0;

                 for (var k=0;k<iterazioni;k++)

                      {

                      x_new = (x*x)-(y*y)+a;

                      y_new = (2*x*y)+b;

                      if (((x_new*x_new)+(y_new*y_new))>4)

                            {

                              // colora il punto

                              r = k%12;

                              context.beginPath();

                              context.fillRect(i-1,j-1,1,1);

                              context.fillStyle=colori[r];

                              context.stroke();

                              break;

                            }

                      x = x_new;

                      y = y_new;

                      }

             }

           }

      }

}

}

}

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;
}