/*       Adam Majewski   adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero             Structure of a program or how to analyze the program          ============== Image X ========================      DrawImageOfX -> DrawPointOfX -> ComputeColorOfX       first 2 functions are identical for every X   check only last function =  ComputeColorOfX   which computes color of one pixel !              ==========================================       ---------------------------------   indent d.c    default is gnu style    -------------------      c console progam       export  OMP_DISPLAY_ENV="TRUE"	   gcc d.c -lm -Wall -march=native -fopenmp   time ./a.out > b.txt     gcc d.c -lm -Wall -march=native -fopenmp     time ./a.out    time ./a.out >i.txt   time ./a.out >e.txt                     convert -limit memory 1000mb -limit disk 1gb dd30010000_20_3_0.90.pgm -resize 2000x2000 10.png           */  #include <stdio.h> #include <stdlib.h>		// malloc #include <string.h>		// strcat #include <math.h>		// M_PI; needs -lm also #include <complex.h> #include <omp.h>		// OpenMP #include <limits.h>		// Maximum value for an unsigned long long int    // https://sourceforge.net/p/predef/wiki/Standards/  #if defined(__STDC__) #define PREDEF_STANDARD_C_1989 #if defined(__STDC_VERSION__) #if (__STDC_VERSION__ >= 199409L) #define PREDEF_STANDARD_C_1994 #endif #if (__STDC_VERSION__ >= 199901L) #define PREDEF_STANDARD_C_1999 #endif #endif #endif     /* --------------------------------- global variables and consts ------------------------------------------------------------ */    // virtual 2D array and integer ( screen) coordinate // Indexes of array starts from 0 not 1  //unsigned int ix, iy; // var static unsigned int ixMin = 0;	// Indexes of array starts from 0 not 1 static unsigned int ixMax;	// static unsigned int iWidth;	// horizontal dimension of array  static unsigned int iyMin = 0;	// Indexes of array starts from 0 not 1 static unsigned int iyMax;	//  static unsigned int iHeight = 10000;	//   // The size of array has to be a positive constant integer  static unsigned long long int iSize;	// = iWidth*iHeight;   // memmory 1D array  unsigned char *data; unsigned char *edge; //unsigned char *edge2;  // unsigned int i; // var = index of 1D array //static unsigned int iMin = 0; // Indexes of array starts from 0 not 1 static unsigned int iMax;	// = i2Dsize-1  =  // The size of array has to be a positive constant integer  // unsigned int i1Dsize ; // = i2Dsize  = (iMax -iMin + 1) =  ;  1D array with the same size as 2D array    // see SetPlane  double radius = 1.2;  complex double center = 0.0; double  DisplayAspectRatio  = 1.3; // https://en.wikipedia.org/wiki/Aspect_ratio_(image) // dx = dy compare setup : iWidth = iHeight; double ZxMin; //= -1.3;	//-0.05; double ZxMax;// = 1.3;	//0.75; double ZyMin;// = -1.3;	//-0.1; double ZyMax;// = 1.3;	//0.7; double PixelWidth;	// =(ZxMax-ZxMin)/ixMax; double PixelHeight;	// =(ZyMax-ZyMin)/iyMax;  double ratio;    /*   ER = pow(10,ERe);   AR = pow(10,-ARe); */ //int ARe ;			// increase ARe until black ( unknown) points disapear  //int ERe ; double ER; double ER2;			//= 1e60; double AR; // bigger values do not works double AR2; double AR12;    int IterMax = 100000;   /* colors = shades of gray from 0 to 255      unsigned char colorArray[2][2]={{255,231},    {123,99}};    color = 245;  exterior  */ unsigned char iColorOfExterior = 245; unsigned char iColorOfInterior1 = 99; unsigned char iColorOfInterior2 = 183; unsigned char iColorOfBoundary = 0; unsigned char iColorOfUnknown = 5;  // pixel counters unsigned long long int uUnknown = 0; unsigned long long int uInterior = 0; unsigned long long int uExterior = 0;    // critical points complex double zc1a = -0.47068779553447764874 + 0.0026098248687148155323*I; //period 1 attract complex double zc1p = 0.24496023578261891251 + 0.0050657319837705428595*I; // period 1  attracting from parabolic  // periodic points = attractors complex double zp1a =-0.33036439123272171026-0.035156244692189524137*I ; //period 1 attract complex double zp1p = 0.024367377494104072722 +0.043208283893844831591*I ; // period 1 attracting from parabolic    /*    F(z,C):=z(1+z)(1+z)(C-(2C+1)z+(3+4C)zz/4)     Then for any C, F(z,C) has the following properties:    1. z=0 is a fixed point with multiplier C.   2. z=-1 is a critical point and F(-1,C)=0.   3. z=1 is another critical point, and F(1,C)=-1, thus F(F(1,C),C)=0       C*z^5+(3*z^5)/4+z^4/2-2*C*z^3-(5*z^3)/4-z^2+C*z    C*z^5   +(3*z^5)/4   +z^4/2   -2*C*z^3   -(5*z^3)/4   -z^2   +C*z    ==============   (%o12) C*z^5+(3*z^5)/4+z^4/2-2*C*z^3-(5*z^3)/4-z^2+C*z    (%i13) coeff(f,z,5);   (%o13) C+3/4 = 0.058*%i+1.77     (%i14) coeff(f,z,4);   (%o14) 1/2    (%i15) coeff(f,z,3);   (%o15) (-2*C)-5/4 = (-0.116*%i)-3.29     (%i16) coeff(f,z,2);   (%o16) -1    (%i17) coeff(f,z,1);   (%o17) C =  1.02 + 0.05*I;     (%i18) coeff(f,z,0);   (%o18) 0     =============================    coefficients read from input file kawahira_sc_c3.txt   degree 5 coefficient = ( +1.7700000000000000 +0.0580000000000000*i)    degree 4 coefficient = ( +0.5000000000000000 +0.0000000000000000*i)    degree 3 coefficient = ( -3.2900000000000000 -0.1160000000000000*i)    degree 2 coefficient = ( -1.0000000000000000 +0.0000000000000000*i)    degree 1 coefficient = ( +1.0200000000000000 +0.0500000000000000*i)    degree 0 coefficient = ( +0.0000000000000000 +0.0000000000000000*i)     Input polynomial p(z)=(1.7700000000000000178+0.058000000000000002942i)*z^5+(0.5+0i)*z^4+(-3.2900000000000000355-0.11600000000000000588i)*z^3+(-1+0i)*z^2+(1.0200000000000000178+0.050000000000000002776i)*z^1    3 critical points found    cp#0: -0.47068779553447764874,0.0026098248687148155323 . It's critical orbit is bounded and enters cycle #0 length=1 and it's stability = |multiplier|=0.66225 =attractive    internal angle = 0.97054734997537162045   cycle = {   -0.33036439123272171026,-0.035156244692189524137 ; }    cp#1: 0.24496023578261891251,0.0050657319837705428595 . It's critical orbit is bounded and enters cycle #1 length=1 and it's stability = |multiplier|=0.98594 =attractive    internal angle = 0.99082387912800862217   cycle = {   0.024367377494104072722,0.043208283893844831591 ; }    cp#2: -1.0000287546513304537,-0.00068489313584920526092 . It's critical orbit is bounded  and enters cycle #1    */    // C=1.02+0.05i complex double C = 1.02 + 0.05*I;   /* ------------------------------------------ functions -------------------------------------------------------------*/      //------------------complex numbers -----------------------------------------------------      // from screen to world coordinate ; linear mapping // uses global cons double GiveZx (int ix) {   return (ZxMin + ix * PixelWidth); }  // uses globaal cons double GiveZy (int iy) {   return (ZyMax - iy * PixelHeight); }				// reverse y axis   complex double GiveZ (int ix, int iy) {   double Zx = GiveZx (ix);   double Zy = GiveZy (iy);    return Zx + Zy * I;     }    double cabs2(complex double z){    return creal(z)*creal(z)+cimag(z)*cimag(z);   }     // ===================== int IsPointInsideTrap1(complex double  z){  	 	  	   if ( cabs2(z - zp1a) < AR2) {return 1;} // circle with prabolic point zp on it's boundary   return 0; // outside    }    // ===================== int IsPointInsideTrap2(complex double  z){  	   if (cabs2(z - zp1p) <AR2) {return 1;} // circle around periodic point 	   return 0; // outside    }   complex double F(complex double z, complex double C){     return (z*(1+z)*(1+z)*(C-(2*C+1)*z+(3+4*C)*z*z/4));   };       // ****************** DYNAMICS = trap tests ( target sets) ****************************   /* -----------  array functions = drawing -------------- */  /* gives position of 2D point (ix,iy) in 1D array  ; uses also global variable iWidth */ unsigned int Give_i (unsigned int ix, unsigned int iy) {   return ix + iy * iWidth; }    // f(z)=1+z−3z2−3.75z3+1.5z4+2.25z5 unsigned char ComputeColor_Fatou (complex double z, int IterMax) {    	 	   double r2;     int i;			// number of iteration   for (i = 0; i < IterMax; ++i)     {   		        z = F(z,C);		// complex iteration f(z)=z^6+A*z+c       r2 =cabs2(z); 		       if (r2 > ER2) // esaping = exterior 	{ 	  uExterior += 1; 	  return iColorOfExterior; 	}			 	       // solid color for each Fatou components 	       if ( IsPointInsideTrap1(z)) { 	uInterior +=1; 	return iColorOfInterior1;       } // 50 + (i % 114); } 	       if (IsPointInsideTrap2(z)){ 	uInterior +=1; 	return iColorOfInterior2;}  	      }    uUnknown += 1;   return iColorOfUnknown;   }      // plots raster point (ix,iy)  int DrawFatouPoint (unsigned char A[], int ix, int iy, int IterMax) {   int i;			/* index of 1D array */   unsigned char iColor = 0;   complex double z;     i = Give_i (ix, iy);		/* compute index of 1D array from indices of 2D array */   z = GiveZ (ix, iy);   iColor = ComputeColor_Fatou (z, IterMax);   A[i] = iColor;		// interior    return 0; }     // fill array  // uses global var :  ... // scanning complex plane  int DrawFatouImage (unsigned char A[], int IterMax) {   unsigned int ix, iy;		// pixel coordinate     fprintf (stdout, "compute Fatou image \n");   // for all pixels of image  #pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax, uUnknown, uInterior, uExterior)   for (iy = iyMin; iy <= iyMax; ++iy)     {       fprintf (stderr, " %d from %d \r", iy, iyMax);	//info        for (ix = ixMin; ix <= ixMax; ++ix) 	DrawFatouPoint (A, ix, iy, IterMax);	//       }    return 0; }   //=========    int IsInside (int x, int y, int xcenter, int ycenter, int r){  	   double dx = x- xcenter;   double dy = y - ycenter;   double d = sqrt(dx*dx+dy*dy);   if (d<r)      return 1;   return 0; 	    }   int PlotBigPoint(complex double z, unsigned char A[]){  	   unsigned int ix_seed = (creal(z)-ZxMin)/PixelWidth;   unsigned int iy_seed = (ZyMax - cimag(z))/PixelHeight;   unsigned int i; 	 	   /* mark seed point by big pixel */   int iSide =1.0*iWidth/4000.0 ; /* half of width or height of big pixel */   int iY;   int iX;   for(iY=iy_seed-iSide;iY<=iy_seed+iSide;++iY){      for(iX=ix_seed-iSide;iX<=ix_seed+iSide;++iX){        if (IsInside(iX, iY, ix_seed, iy_seed, iSide)) { 	i= Give_i(iX,iY); /* index of _data array */ 	A[i]= 255-A[i];}}} 	 	   return 0; 	 }   // fill array  // uses global var :  ... // scanning complex plane  int MarkAttractors (unsigned char A[]) {    	 	 	   fprintf (stderr, "mark attractors \n");      PlotBigPoint(zp1a, A); // period 114  cycle   PlotBigPoint(zp1p, A);	// period 2 attracting cycle     		        	    return 0; }   // ===================== int IsPointInsideTraps(unsigned int ix, unsigned int iy){  	   complex double  z = GiveZ (ix, iy); 	   if ( IsPointInsideTrap1(z)) {return 1;} // circle with prabolic point on it's boundary 	   if (IsPointInsideTrap2(z)) {return 1;} 	   return 0; // outside    }      int MarkTraps(unsigned char A[]){    unsigned int ix, iy;		// pixel coordinate    unsigned int i;     fprintf (stderr, "Mark traps \n");   // for all pixels of image  #pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax, uUnknown, uInterior, uExterior)   for (iy = iyMin; iy <= iyMax; ++iy)     {       fprintf (stderr, " %d from %d \r", iy, iyMax);	//info        for (ix = ixMin; ix <= ixMax; ++ix){ 	if (IsPointInsideTraps(ix, iy)) { 	  i= Give_i(ix,iy); /* index of _data array */ 	  A[i]= 255-A[i]; // inverse color 	}}}   return 0; }       int PlotPoint(complex double z, unsigned char A[]){  	   unsigned int ix = (creal(z)-ZxMin)/PixelWidth;   unsigned int iy = (ZyMax - cimag(z))/PixelHeight;   unsigned int i = Give_i(ix,iy); /* index of _data array */ 	 	   A[i]= 255-A[i]; // Mark point with inveres color 	 	   return 0; 	 }   int DrawForwardOrbit(complex double z, unsigned long long int iMax,  unsigned char A[] ) {      unsigned long long int i; /* nr of point of critical orbit */        PlotBigPoint(z, A);      /* forward orbit of critical point  */   for (i=1;i<iMax ; ++i)     {       z  = F(z,C);       if (cabs2(z - zp1p) > 2.0) {return 1;} // escaping       PlotBigPoint(z, A);     }                return 0;   }        // *********************************************************************************************** // ********************** edge detection usung Sobel filter *************************************** // ***************************************************************************************************  // from Source to Destination int ComputeBoundaries(unsigned char S[], unsigned char D[]) {     unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */   unsigned int i; /* index of 1D array  */   /* sobel filter */   unsigned char G, Gh, Gv;    // boundaries are in D  array ( global var )     // clear D array   memset(D, iColorOfExterior, iSize*sizeof(*D)); // for heap-allocated arrays, where N is the number of elements = FillArrayWithColor(D , iColorOfExterior);     // printf(" find boundaries in S array using  Sobel filter\n");    #pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax)   for(iY=1;iY<iyMax-1;++iY){      for(iX=1;iX<ixMax-1;++iX){        Gv= S[Give_i(iX-1,iY+1)] + 2*S[Give_i(iX,iY+1)] + S[Give_i(iX-1,iY+1)] - S[Give_i(iX-1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[Give_i(iX+1,iY-1)];       Gh= S[Give_i(iX+1,iY+1)] + 2*S[Give_i(iX+1,iY)] + S[Give_i(iX-1,iY-1)] - S[Give_i(iX+1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[Give_i(iX-1,iY-1)];       G = sqrt(Gh*Gh + Gv*Gv);       i= Give_i(iX,iY); /* compute index of 1D array from indices of 2D array */       if (G==0) {D[i]=255;} /* background */       else {D[i]=0;}  /* boundary */     }   }           return 0; }    // copy from Source to Destination int CopyBoundaries(unsigned char S[],  unsigned char D[]) {     unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */   unsigned int i; /* index of 1D array  */       //printf("copy boundaries from S array to D array \n");   for(iY=1;iY<iyMax-1;++iY)     for(iX=1;iX<ixMax-1;++iX)       {i= Give_i(iX,iY); if (S[i]==0) D[i]=0;}         return 0; }                 // ******************************************************************************************* // ********************************** save A array to pgm file **************************** // *********************************************************************************************  int SaveArray2PGMFile (unsigned char A[], int a, int b,  int c, char *comment) {    FILE *fp;   const unsigned int MaxColorComponentValue = 255;	/* color component is coded from 0 to 255 ;  it is 8 bit color file */   char name[100];		/* name of file */   snprintf (name, sizeof name, "%d_%d_%d", a, b, c );	/*  */   char *filename = strcat (name, ".pgm");   char long_comment[200];   sprintf (long_comment, "fc(z)=z*(1+z)*(1+z)*(C-(2C+1)*z+(3+4*C)*z*z/4)  %s", comment);        // save image array to the pgm file    fp = fopen (filename, "wb");	// create new file,give it a name and open it in binary mode    fprintf (fp, "P5\n # %s\n %u %u\n %u\n", long_comment, iWidth, iHeight, MaxColorComponentValue);	// write header to the file   size_t rSize = fwrite (A, sizeof(A[0]), iSize, fp);	// write whole array with image data bytes to the file in one step    fclose (fp);    // info    if ( rSize == iSize)    	{   		printf ("File %s saved ", filename);   		if (long_comment == NULL || strlen (long_comment) == 0)     		printf ("\n");   			else { printf (". Comment = %s \n", long_comment); }   	}   	else {printf("wrote %zu elements out of %llu requested\n", rSize,  iSize);}    return 0; }     int PrintCInfo () {    printf ("gcc version: %d.%d.%d\n", __GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__);	// https://stackoverflow.com/questions/20389193/how-do-i-check-my-gcc-c-compiler-version-for-my-eclipse   // OpenMP version is displayed in the console : export  OMP_DISPLAY_ENV="TRUE"    printf ("__STDC__ = %d\n", __STDC__);   printf ("__STDC_VERSION__ = %ld\n", __STDC_VERSION__);   printf ("c dialect = ");   switch (__STDC_VERSION__)     {				// the format YYYYMM      case 199409L:       printf ("C94\n");       break;     case 199901L:       printf ("C99\n");       break;     case 201112L:       printf ("C11\n");       break;     case 201710L:       printf ("C18\n");       break;       //default : /* Optional */      }    return 0; }   int PrintProgramInfo () {     // display info messages   printf ("Numerical approximation of Julia set for F(z,C) =  z*(1+z)*(1+z)*(C-(2C+1)*z+(3+4*C)*z*z/4) \n");   //printf ("iPeriodParent = %d \n", iPeriodParent);   //printf ("iPeriodOfChild  = %d \n", iPeriodChild);      printf ("parameter C = ( %.16f ; %.16f ) \n", creal (C), cimag (C));          printf ("Image Width = %f in world coordinate\n", ZxMax - ZxMin);   printf ("PixelWidth = %.16f \n", PixelWidth);   printf ("AR = %.16f = %f *PixelWidth\n", AR, AR / PixelWidth);     printf("pixel counters\n");   printf ("uUnknown = %llu\n", uUnknown);   printf ("uExterior = %llu\n", uExterior);   printf ("uInterior = %llu\n", uInterior);   printf ("Sum of pixels  = %llu\n", uInterior+uExterior + uUnknown);   printf ("all pixels of the array = iSize = %llu\n", iSize);     // image corners in world coordinate   // center and radius   // center and zoom   // GradientRepetition   printf ("Maximal number of iterations = iterMax = %d \n", IterMax);   printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);   //       return 0; }    int SetPlane(complex double center, double radius, double a_ratio){    ZxMin = creal(center) - radius*a_ratio;	   ZxMax = creal(center) + radius*a_ratio;	//0.75;   ZyMin = cimag(center) - radius;	// inv   ZyMax = cimag(center) + radius;	//0.7;   return 0;  }   // ***************************************************************************** //;;;;;;;;;;;;;;;;;;;;;;  setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; // **************************************************************************************  int setup () {    fprintf (stderr, "setup start\n");         /* 2D array ranges */    iWidth = iHeight* DisplayAspectRatio ;   iSize = iWidth * iHeight;	// size = number of points in array    // iy   iyMax = iHeight - 1;		// Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].   //ix    ixMax = iWidth - 1;    /* 1D array ranges */   // i1Dsize = i2Dsize; // 1D array with the same size as 2D array   iMax = iSize - 1;		// Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].     SetPlane( center, radius,  DisplayAspectRatio );	   /* Pixel sizes */   PixelWidth = (ZxMax - ZxMin) / ixMax;	//  ixMax = (iWidth-1)  step between pixels in world coordinate    PixelHeight = (ZyMax - ZyMin) / iyMax;   ratio = ((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((double) iWidth / (double) iHeight);	// it should be 1.000 ...    ER = 2.0; //    ER2 = ER*ER;   AR = 18.5*PixelWidth*iWidth/2000.0 ; // adjust first number    AR2 = AR * AR;   //AR12 = AR/2.0;         	            /* create dynamic 1D arrays for colors ( shades of gray ) */   data = malloc (iSize * sizeof (unsigned char));    edge = malloc (iSize * sizeof (unsigned char));   if (data == NULL || edge == NULL)     {       fprintf (stderr, " Could not allocate memory");       return 1;     }            fprintf (stderr, " end of setup \n");    return 0;  }				// ;;;;;;;;;;;;;;;;;;;;;;;;; end of the setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;     int end () {     fprintf (stderr, " allways free memory (deallocate )  to avoid memory leaks \n");	// https://en.wikipedia.org/wiki/C_dynamic_memory_allocation   free (data);   free(edge);     PrintProgramInfo ();   PrintCInfo ();   return 0;  }  // ******************************************************************************************************************** /* -----------------------------------------  main   -------------------------------------------------------------*/ // ********************************************************************************************************************  int main () {   setup ();     DrawFatouImage (data, IterMax);	// first find Fatou   SaveArray2PGMFile (data, iWidth, IterMax, 0, "Fatou, name = iWidth_IterMax_n");      ComputeBoundaries(data,edge);   SaveArray2PGMFile (edge, iWidth, IterMax, 1, "Boundaries of Fatou; name = iWidth_IterMax_n");       CopyBoundaries(edge,data);   SaveArray2PGMFile (data, iWidth, IterMax, 2, "Fatou with boundaries; name = iWidth_IterMax_n");       //MarkAttractors(data);   MarkTraps(data);   SaveArray2PGMFile (data, iWidth, IterMax, 4, "Fatou with boundaries and traps; name = iWidth_IterMax_n");    	   DrawForwardOrbit(zc1p, 10000, data);   SaveArray2PGMFile (data, iWidth, IterMax, 5, "Fatou with boundaries and traps, critical orbit; name = iWidth_IterMax_n");     	    end ();    return 0; }