Buddhabrot konečně v C

Tak jsem se konečně rozhoupal k tomu, abych v C napsal generátor buddhabrota, který je zobrazen skrze knihovnu SDL. Je to lepší než (pomalým) pythonem sázet do 3D světa blenderu vertexy, které mají efekt halo a dostatečně velkou alfu a pak to celé ještě renderovat jako klasickou 3D scénu blenderu... No prostě neřešil jsem to úplně čistě a optimálně. Navíc po několika stech milionů vertexů se Blenderu přestal můj záměr líbit a leckdy už ani nebyl schopen nastartovat renderer...

První, co jsem tedy po svém programu chtěl byla prakticky neomezená doba běhu, dále vysoká preciznost a pokud možno i rychlost (kterou nejsem schopen zhodnotit, protože neznám úplně optimální algoritmy). Toho jsem docílil tím, že všechny buňky rastru, do kterého se sází iterované body, jsou typu "unsigned long long int" (opravdu hodně hodně dlouhý integer bez znaménka), komplexní čísla mají "long double" preciznost a při hledání bodů nenáležících Mandelbrotově množině se jednoduchým matematickým vzorcem odhazují všechny body z tzv. "hlavního kardioidu", což je ta největší část M-množiny, ta ve tvaru srdce.
Zde je tedy zdrojový kód:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <sdl.h>
#include <time.h>

/* Render plane */
#define WIDTH    1600
#define HEIGHT   1200

#define R_MIN   -2.0
#define R_MAX    1.0
#define I_MIN   -1.25
#define I_MAX    1.25

/* Z pool */
#define ZR_MIN  -4.0
#define ZR_MAX   4.0
#define ZI_MIN  -4.0
#define ZI_MAX   4.0

/* Math */
#define MIN_ITER_R 50000
#define MAX_ITER_R 100000
#define MIN_ITER_G 100
#define MAX_ITER_G 50000
#define MIN_ITER_B 1
#define MAX_ITER_B 100
#define BAILOUT  64

/* SDL */
#define DEPTH   32
#define BMP_FILE "C:\\buddha.bmp"

/* Color */
#define FACTOR_R 0.25
#define FACTOR_G 0.35
#define FACTOR_B 0.85

/* Global variables for everyday use */
unsigned int width=WIDTH;
unsigned int height=HEIGHT;
long double r_min=R_MIN;
long double r_max=R_MAX;
long double i_min=I_MIN;
long double i_max=I_MAX;
long double zr_min=ZR_MIN;
long double zr_max=ZR_MAX;
long double zi_min=ZI_MIN;
long double zi_max=ZI_MAX;
unsigned long long int min_iter_r=MIN_ITER_R;
unsigned long long int max_iter_r=MAX_ITER_R;
unsigned long long int min_iter_g=MIN_ITER_G;
unsigned long long int max_iter_g=MAX_ITER_G;
unsigned long long int min_iter_b=MIN_ITER_B;
unsigned long long int max_iter_b=MAX_ITER_B;
float bailout=BAILOUT;

char bmp_file[]=BMP_FILE;

/* Used for determining color */
unsigned long long int maximum_r=1;
unsigned long long int maximum_g=1;
unsigned long long int maximum_b=1;
unsigned int minimum_r=0;
unsigned int minimum_g=0;
unsigned int minimum_b=0;
float factor_r=FACTOR_R;
float factor_g=FACTOR_G;
float factor_b=FACTOR_B;


/* Cycler accepts a starting point "c" and 2D array where the changes will be made */
int cycler(long double complex c, unsigned long long int **array[3]);
/* printf_array is plain-text output of generated array */
void printf_array(unsigned long long int **array[3]);
/* redraw accepts a 2D array and vizualize its values on SDL surface provided */
void redraw(unsigned long long int **array[3],SDL_Surface *screen);
/* Function for freeing array */
void free_array(unsigned long long int **array[3]);
/* SDL function copied from demo */
void setpixel(SDL_Surface *screen, int x, int y, Uint8 r, Uint8 g, Uint8 b);
/* Checks if c is inside M-set; if no, returns n of iters needed to escape */
unsigned long long int inside(long double complex c);
/* Coloring function accepts number of hits */

int main(int argc,char *argv[])
{ 
 /* Create array (with little GOTo rebellion :D ) */
 unsigned long long int **array[3];
 array[0] = malloc(width * sizeof(unsigned long long int *));
 if(array[0]==NULL) goto malloc_fail_1;
 array[1] = malloc(width * sizeof(unsigned long long int *));
 if(array[1]==NULL) goto malloc_fail_2;
 array[2] = malloc(width * sizeof(unsigned long long int *));
 if(array[2]==NULL) goto malloc_fail_3;
 goto malloc_nofail;
 malloc_fail_3:
  free(array[1]);
 malloc_fail_2:
  free(array[0]);
 malloc_fail_1:
  fprintf(stderr,"Out of memory! Can't create output array.\n");
  return 1;
 malloc_nofail:
 /* Let's fill the R channel with pointers to actual cols */
 for(unsigned int i=0;i<width;i++)
 {
  array[0][i]=malloc(height * sizeof(unsigned long long int));
  if(array[0][i] == NULL)
  {
   fprintf(stderr,"Out of memory! Can't create output array.\n");
   /* Rollback */
   for(;i>0;free(array[0][--i]));
   free(array[0]);
   return 1;
  }
 }
 /* Do the same with G channel */
 for(unsigned int i=0;i<width;i++)
 {
  array[1][i]=malloc(height * sizeof(unsigned long long int));
  if(array[1][i] == NULL)
  {
   fprintf(stderr,"Out of memory! Can't create output array.\n");
   /* Rollback */
   for(;i>0;free(array[1][--i]));
   free(array[1]);
   for(i=0;i<width;free(array[0][i++]));
   free(array[0]);
   return 1;
  }
 }
 /* And B channel */
 for(unsigned int i=0;i<width;i++)
 {
  array[2][i]=malloc(height * sizeof(unsigned long long int));
  if(array[2][i] == NULL)
  {
   fprintf(stderr,"Out of memory! Can't create output array.\n");
   /* Rollback */
   for(;i>0;free(array[2][--i]));
   free(array[2]);
   for(i=0;i<width;free(array[1][i++]));
   free(array[1]);
   for(i=0;i<width;free(array[0][i++]));
   free(array[0]);
   return 1;
  }
 }

 
 /* Clear array (fill with zeros) */
 for(unsigned int x=0;x<width;x++)
 {
  for(unsigned int y=0;y<height;y++)
  {
   array[0][x][y]=0;
   array[1][x][y]=0;
   array[2][x][y]=0;
  }
 }

 /* Now lets prepare our output */
 SDL_Surface *screen;
 if (SDL_Init(SDL_INIT_VIDEO) < 0 )
 {
  fprintf(stderr,"Can't initialize graphics!\n");
  free_array(array);
  return 1;
 }
 if (!(screen = SDL_SetVideoMode(width,height,DEPTH,/*SDL_FULLSCREEN|*/SDL_HWSURFACE)))
 {
  SDL_Quit();
  free_array(array);
  return 1;
 }

 /* and prepare SDL quit handling */
 SDL_Event event;
 int keypress=0;

 /* If we don't want always the same set of starting c's, then we need to init rand */
 srand(time(NULL));
 
 /* Everything is prepared, we can start shooting */
 for(unsigned long long int i=0;!keypress;i++)
 {
  //long double complex c=((zr_max-zr_min)*((long double)rand()/RAND_MAX)+zr_min)+((zi_max-zi_min)*((long double)rand()/RAND_MAX)+zi_min)*I;
  cycler(((zr_max-zr_min)*(rand()/(long double)RAND_MAX)+zr_min)+((zi_max-zi_min)*(rand()/(long double)RAND_MAX)+zi_min)*I,array);
  /* Time to time, we show out results */
  if(i%2000==0)
  {
   printf("%llu random complex numbers cycled through the cycler.\n",i);
   redraw(array,screen);
   /* May we end? (Was SDL Window closed?) */
   SDL_PollEvent(&event);
   switch (event.type) 
   {
    case SDL_QUIT:
     keypress = 1;
     break;
   }
  }
 }
 
 /* When we reached this, main loop ended and we're ready to see the result */
 //printf_array(array);

 /* Saving of generaded picture, just some SDL function */
 SDL_SaveBMP(screen,bmp_file);

 SDL_Quit();
 /* At the end, free the memory */
 free_array(array);
 /* Issue success */
 return 0;
}

int cycler(long double complex c, unsigned long long int **array[3])
{
 /* Is c outside of M-set? */
 int j;
 if((j=inside(c))==0)
  return 0;
 /* So now we know c is outside and interesting enough, start shooting live */
 /* Decide which channel we will be shooting to */
 int color=0;
 if((min_iter_r<j)&&(j<max_iter_r)) color+=1;
 if((min_iter_g<j)&&(j<max_iter_g)) color+=2;
 if((min_iter_b<j)&&(j<max_iter_b)) color+=4;
 /* If it is not in any channel, we don't want this point */
 if(color==0) return 0;

 long double complex Z=c;
 for(unsigned int i=0;i<j;i++)
 {
  Z=Z*Z+c;
  if(cabsl(Z)>bailout)
   return i;
  signed int x=(creall(Z)-r_min)/(r_max-r_min)*width;
  signed int y=(cimagl(Z)-i_min)/(i_max-i_min)*height;
  if(0<=x&&x<width&&0<=y&&y<height)
  {
   if(color&1)
   {
    if(++array[0][x][y]>maximum_r)
    {
     maximum_r++;
    }
   }
   if(color&2)
   {
    if(++array[1][x][y]>maximum_g)
    {
     maximum_g++;
    }
   }
   if(color&4)
   {
    if(++array[2][x][y]>maximum_b)
    {
     maximum_b++;
    }
   }
  }
 }

 return -1; /* This shouldn't happen */
}

void printf_array(unsigned long long int **array[3])
{
 /* We will start with x, so the buddhabrot will be "standing" and walking through array faster */
 printf("Red:\n");
 for(int x=0;x<width;x++)
 {
  for(int y=0;y<height;y++)
   printf("%llu ",array[0][x][y]);
  printf("\n");
 }
 printf("Maximum: %llu\n",maximum_r);
 printf("Green:\n");
 for(int x=0;x<width;x++)
 {
  for(int y=0;y<height;y++)
   printf("%llu ",array[1][x][y]);
  printf("\n");
 }
 printf("Maximum: %llu\n",maximum_g);
 printf("Blue:\n");
 for(int x=0;x<width;x++)
 {
  for(int y=0;y<height;y++)
   printf("%llu ",array[2][x][y]);
  printf("\n");
 }
 printf("Maximum: %llu\n",maximum_b);
}

void redraw(unsigned long long int **array[3],SDL_Surface *screen)
{
 /* Find global minimum */
 minimum_r=maximum_r;
 for(int x=0;x<width;x++)
  for(int y=0;y<height;y++)
   minimum_r=(array[0][x][y]<minimum_r)?array[0][x][y]:minimum_r;
 minimum_g=maximum_g;
 for(int x=0;x<width;x++)
  for(int y=0;y<height;y++)
   minimum_g=(array[1][x][y]<minimum_g)?array[1][x][y]:minimum_g;
 minimum_b=maximum_b;
 for(int x=0;x<width;x++)
  for(int y=0;y<height;y++)
   minimum_b=(array[2][x][y]<minimum_b)?array[2][x][y]:minimum_b;
 /* fill screen with pixels */
 for(int x=0;x<width;x++)
 {
  for(int y=0;y<height;y++)
  {
   int r=255*(pow( ((array[0][x][y]-minimum_r)/(float)(maximum_r-minimum_r)) ,factor_r));
   int g=255*(pow( ((array[1][x][y]-minimum_g)/(float)(maximum_g-minimum_g)) ,factor_g));
   int b=255*(pow( ((array[2][x][y]-minimum_b)/(float)(maximum_b-minimum_b)) ,factor_b));
   setpixel(screen,x,y,r,g,b);
  }
 }
    if(SDL_MUSTLOCK(screen)) SDL_UnlockSurface(screen);
    SDL_Flip(screen); 
}

void free_array(unsigned long long int **array[3])
{
 for(int i=0;i<width;i++)
 {
  free(array[0][i]);
  free(array[1][i]);
  free(array[2][i]);
 }
 free(array[0]);
 free(array[1]);
 free(array[2]);
}

void setpixel(SDL_Surface *screen, int x, int y, Uint8 r, Uint8 g, Uint8 b)
{
 Uint32 *pixmem32;
 Uint32 colour;  

 colour = SDL_MapRGB( screen->format, r, g, b );
 pixmem32 = (Uint32*) screen->pixels  + y*width + x;
 *pixmem32 = colour;

}

unsigned long long int inside(long double complex c)
{
 /* First, little optimization */
 if((creall(c)<sqrtl(powl(creall(c)-0.25,2)+powl(cimagl(c),2))-2*(powl(creall(c)-0.25,2)+powl(cimagl(c),2))+0.25)||((powl(creall(c)+1,2)+powl(cimagl(c),2))<(1/16)))
  return 0;
 /* If still in game, check conventionally */
 long double complex Z=c;
 /* Find max_iter */
 unsigned long long int max_iter=(max_iter_r>max_iter_g)?(max_iter_b>max_iter_r)?max_iter_b:max_iter_r:max_iter_g;
 /* Iterate! */
 for(int i=0;i<max_iter;i++)
 {
  Z=Z*Z+c;
  /* We're just waiting if Z exits after a few iterations or not */
  if(cabsl(Z)>bailout)
   return i;
 }
 return 0;
}

Tak a tady jsou výsledky. Někdy to budu muset nechat běžet pěkně dlouho, aby to něco zajímavého udělalo.

Všechny hodnoty iterací nastaveny stejně, takže červená, želená i modrá mají stejný obsah -> obrázek je černobílý

Tady má modrá nejmenší počet iterací, zelená trochu větší a červená největší.

Tohle bych chtěl na tričko ;)

Teď ještě plánuju přidat parsování command-line argumentů a bude hotovo :) Kdybyste měli jakékoliv připomínky k programu, neváhejte mi je sdělit, budu jedině rád.