/*

  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

  
  =======================
  # gnuplot "i.plt"
set terminal svg enhanced background rgb 'white'
set xlabel "re(z)"
set ylabel "DLD"
set title "Relation between z and DLD in the interior of Julia set for c = -1"
set output "interior.svg"
plot "i.txt" with lines

  ----------------------
  d0 - db  = 5.0000000000000000 - 4.5389870050569598 = 0.4610129949430402
 allways free memory (deallocate )  to avoid memory leaks 
Numerical approximation of Julia set for fc(z)= z^2 + c 
parameter c = ( -1.0000000000000000 ; 0.0000000000000000 ) 
Image Width = 4.000000 in world coordinate
PixelWidth = 0.004004 
Maximal number of iterations = iterMax = 1000 
ratio of image  = 1.000000 ; it should be 1.000 ...
gcc version: 7.5.0
  
  
*/

#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


static const double ZxMin = -1.8;	//-0.05;
static const double ZxMax = 1.8;	//0.75;
static const double ZyMin = -1.8;	//-0.1;
static const double ZyMax = 1.8;	//0.7;
static double PixelWidth;	// =(ZxMax-ZxMin)/ixMax;
static double PixelHeight;	// =(ZyMax-ZyMin)/iyMax;
static double ratio;


// complex numbers of parametr plane 
//https://fractalforums.org/code-snippets-fragments/74/lagrangian-descriptors-fragment-code/3612/msg22426#msg22426
double complex c = -1.05204872;	// parameter of function fc(z)=z^2 + c


// attracting period 2 cycle 
double complex z21 = 0.049822582293598;
double complex z22 = -1.049822582293598;
double complex z1a = -0.641073494565534;	// alfa
complex double z1b = 1.641073494565534;	// beta
double complex zb = 0.203208556149733 + 0.377005347593583 * I;	// point on the boubndary of main componnet

/*
ER = pow(10,ERe);
   AR = pow(10,-ARe);
 */
int ARe = 3;			// increase ARe until black ( unknown) points disapear 
int ERe = 3;
double ER;			//= 1e60;
double AR;			//= 1e-16; // bigger values do not works




int IterMax = 1000;

// DLD
const int N = 20;		// fixed number : maximal number of iterations
double p = 0.01444322;		//


// DLD colors
double me = 1.0;
double mi = 0.9;
double d21;			// = lagrangian(z21, c, N, p);
double d22;			//  = lagrangian(z22, c, N, p);
double db;			// = lagrangian(z1a, c, N, p);        
double dd;			// = d1a-d21;



/* colors = shades of gray from 0 to 255 */
unsigned char iColorOfExterior = 150;
unsigned char iColorOfInterior = 50;
unsigned char iColorOfBoundary = 0;
unsigned char iColorOfUnknown = 255;

// pixel counters
unsigned long long int uUnknown = 0;
unsigned long long int uInterior = 0;
unsigned long long int uExterior = 0;




/* ------------------------------------------ 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;




}




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


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







// ***************************************************************************************************************************
// ************************** DLD/J*****************************************
// ****************************************************************************************************************************



/* partial pnorm 
   input: z , zn = f(z), p
   output ppn
   
   
*/
double
ppnorm (complex double z, complex double zn, double p)
{

  double s[2][3];		// array for 2 points on the Riemann sphere
  int j;
  double d;			// denominator 
  double x;
  double y;

  double ds;
  double ppn = 0.0;

  // map from complex plane to riemann sphere
  // z
  x = creal (z);
  y = cimag (z);
  d = x * x + y * y + 1.0;

  s[0][0] = (2.0 * x) / d;
  s[0][1] = (2.0 * y) / d;
  s[0][2] = (d - 2.0) / d;	// (x^2 + y^2 - 1)/d

  // zn
  x = creal (zn);
  y = cimag (zn);
  d = x * x + y * y + 1.0;
  s[1][0] = (2.0 * x) / d;
  s[1][1] = (2.0 * y) / d;
  s[1][2] = (d - 2.0) / d;	// (x^2 + y^2 - 1)/d

  // sum 
  for (j = 0; j < 3; ++j)
    {
      ds = fabs (s[1][j] - s[0][j]);
      //  normal:  neither zero, subnormal, infinite, nor NaN
      //if (fpclassify (ds) !=FP_INFINITE)
      //if (isnormal(ds)) 
      // it is solved by if (cabs(z) > 1e60 ) break; procedure in parent function 
      ppn += pow (ds, p);	// |ds|^p
      //      else {ppn = 10000.0; printf("ds = infty\t");} // 

    }


  return ppn;







}

// DLD = Discret Lagrangian Descriptior
double
lagrangian (complex double z0, complex double c, int iMax, double p)
{

  int i;			// number of iteration
  double d = 0.0;		// DLD = sum
  double ppn;			// partial pnorm
  complex double z = z0;
  complex double zn;		// next z

  for (i = 0; i < iMax; ++i)
    {




      zn = z * z + c;		// complex iteration
      ppn = ppnorm (z, zn, p);
      d += ppn;			// sum
      //
      z = zn;

      //if (! isnormal(d)) { return 0.0; } // not works
      if (cabs (z) > 1e6)
	break;			// exterior : big values produces artifacts on the image  



    }





  //if (d<0.0) {// interior
  // d(z1a) - d(z21) = -0.0804163521959989        
  //      d = - d;
  //      d = (db - d) /dd ; // normalize, see test_interior
  //d = d*d;
  //if (d>1.0) {printf("d int > 1.0\n");
  ///     }
  //      else {

  d = d / ((double) i);		// averaging not summation
  //d = d*me;} // exterior

  return d;




}





unsigned char
ComputeColor_DLD (complex double z, int FatouType)
{


  //double cabsz;
  int iColor;
  double d;

  if (FatouType == 1)
    {				// interior
      d = lagrangian (z, c, N, p);
      // modify gradient position

      //{d = d - (int)d;} // only fractional part
      d = d * d * mi;
      //if ( d< 1.0 ) d = 0.0;

    }				//  
  else
    {
      d = lagrangian (z, c, 10 * N, p);
    }

  iColor = (int) (d * 255) % 255;	// nMax or lower walues in denominator



  return (unsigned char) iColor;


}



// plots raster point (ix,iy) 
int
DrawDLDPoint (unsigned char A[], int ix, int iy)
{
  int i;			/* index of 1D array */
  unsigned char iColor;
  complex double z;
  int FatouType;


  i = Give_i (ix, iy);		/* compute index of 1D array from indices of 2D array */
  z = GiveZ (ix, iy);
  iColor = A[i];		// read color = read the information about Fatou component type ( interior/exterior)
  if (iColor == iColorOfInterior)
    {
      FatouType = 1;
    }				// tru = interior
  else
    {
      FatouType = 0;
    }

  iColor = ComputeColor_DLD (z, FatouType);	// compute new color 
  A[i] = iColor;		// save new colr to the array         

  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int
DrawDLDImage (unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 

  printf ("compute DLD image \n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  for (iy = iyMin; iy <= iyMax; ++iy)
    {
      printf (" %d from %d \r", iy, iyMax);	//info 
      for (ix = ixMin; ix <= ixMax; ++ix)
	DrawDLDPoint (A, ix, iy);	//  
    }

  return 0;
}



//=========================================




unsigned char
ComputeColor_Fatou (complex double z, int IterMax)
{

  int i;			// number of iteration
  for (i = 0; i < IterMax; ++i)
    {




      z = z * z + c;		// complex iteration

      //if (! isnormal(d)) { return 0.0; } // not works
      if (cabs (z) > 2.0)
	{
	  uExterior += 1;
	  return iColorOfExterior;
	}			// exterior : big values produces NAN error in ppnorm computing 


      if (cabs (z - z21) < AR || cabs (z - z22) < AR)
	{			// interior
	  uInterior += 1;
	  return iColorOfInterior;


	}


    }

  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 

  printf ("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)
    {
      printf (" %d from %d \r", iy, iyMax);	//info 
      for (ix = ixMin; ix <= ixMax; ++ix)
	DrawFatouPoint (A, ix, iy, IterMax);	//  
    }

  return 0;
}


//=========







// uses global var :  ...
// scanning complex plane 
int
CheckFatouImage (unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 
  unsigned long long int u_Unknown = 0;
  unsigned long long int u_Interior = 0;
  unsigned long long int u_Exterior = 0;

  printf ("check Fatou array \n");
  // for all pixels of image 
  //#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  for (iy = iyMin; iy <= iyMax; ++iy)
    {
      //printf (" %d from %d \r", iy, iyMax); //info 
      for (ix = ixMin; ix <= ixMax; ++ix)

	{
	  unsigned char color = A[Give_i (ix, iy)];
	  if (color == iColorOfInterior)
	    {
	      u_Interior += 1;
	    }
	  else
	    {
	      if (color == iColorOfExterior)
		{
		  u_Exterior += 1;
		}
	      {
		if (color == iColorOfUnknown)
		  {
		    u_Unknown += 1;
		  }
	      }
	    }
	  //printf("error\n");

	}



    }



  printf ("pixel counters\n");
  printf ("uUnknown = %llu\n", u_Unknown);
  printf ("uExterior = %llu\n", u_Exterior);
  printf ("uInterior = %llu\n", u_Interior);
  printf ("Sum of pixels  = %llu\n", u_Interior + u_Exterior + u_Unknown);
  printf ("all pixels of the array = iSize = %llu\n", iSize);
  printf ("Maximum value for an unsigned long long int = ULLONG_MAX = %llu\n",
	  ULLONG_MAX);

  return 0;
}



//-------------------------------------------------------

// test how values changes to tune color 
// from z21 to zb
int
test_interior ()
{

  complex double z = z21;
  complex double dz = cabs (z21 - zb) / 15.0;
  double r = cabs (zb);


  printf ("# re(z) \t d\n");	// gnuplot
  while (creal (z) < r)
    {				// from z21 to z1a

      double d = lagrangian (z, c, N, p);
      //int iColor = ComputeColorOfDLD(z);

      // printf(" z = %.16f d = %.16f color = %d \n",creal(z), d, iColor);
      printf (" %.16f %.16f \n", creal (z), d);	// gnuplot 
      z += dz;
    }

  //              
  d21 = lagrangian (z21, c, N, p);
  db = lagrangian (zb, c, N, p);
  d22 = lagrangian (z22, c, N, p);
  dd = db - d21;
  printf ("d(z21) = %.16f \n", d21);
  printf ("d(z22) = %.16f \n", d22);
  printf ("d(zb) = %.16f \n", db);
  printf ("d(zb) - d(z21) = %.16f\n", dd);;



  return 0;


}




// test how values changes to tune color 
int
test_exterior ()
{

  complex double z;
  complex double z0 = 1.62;
  complex double z1 = 3.0;
  complex double dz = 0.001;


  z = z0;
  printf ("# z d\n");		// gnuplot
  while (creal (z) < creal (z1))
    {

      double d = lagrangian (z, c, N, p);
      //int iColor = ComputeColorOfDLD(z);

      // printf(" z = %.16f d = %.16f color = %d \n",creal(z), d, iColor);
      printf (" %.16f %.16f \n", creal (z), d);	// gnuplot 
      z += dz;
    }

  //              
  double d0 = lagrangian (z0, c, N, p);
  double d1 = lagrangian (z1, c, N, p);
  double dd = d0 - d1;
  printf ("d0 - d1  = %.16f - %.16f = %.16f\n", d0, d1, dd);

  return 0;


}










// *******************************************************************************************
// ********************************** save A array to pgm file ****************************
// *********************************************************************************************

int
SaveArray2PGMFile (unsigned char A[], int a, int b, int c, double d,
		   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, "dd300%d_%d_%d_%.2f", a, b, c, d);	/*  */
  char *filename = strcat (name, ".pgm");
  char long_comment[200];
  sprintf (long_comment, "%s\tER = %e\tAR =%e", comment, ER, AR);





  // 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
  fwrite (A, iSize, 1, fp);	// write array with image data bytes to the file in one step 
  fclose (fp);

  // info 
  printf ("File %s saved ", filename);
  if (long_comment == NULL || strlen (long_comment) == 0)
    printf ("\n");
  else
    printf (". Comment = %s \n", long_comment);

  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 fc(z)= z^2 + c \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 ("DLD : N = Maximal number of iterations = iterMax = %d \n", N);
  printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  //




  return 0;
}






// *****************************************************************************
//;;;;;;;;;;;;;;;;;;;;;;  setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************

int
setup ()
{

  printf ("setup start\n");






  /* 2D array ranges */

  iWidth = iHeight;
  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].

  /* 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 = pow (10, ERe);
  AR = pow (10, -ARe);



  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc (iSize * sizeof (unsigned char));


  if (data == NULL)
    {
      fprintf (stderr, " Could not allocate memory");
      return 1;
    }





  test_interior ();


  printf (" end of setup \n");

  return 0;

}				// ;;;;;;;;;;;;;;;;;;;;;;;;; end of the setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




int
end ()
{


  printf (" allways free memory (deallocate )  to avoid memory leaks \n");	// https://en.wikipedia.org/wiki/C_dynamic_memory_allocation
  free (data);


  PrintProgramInfo ();
  PrintCInfo ();
  return 0;

}

// ********************************************************************************************************************
/* -----------------------------------------  main   -------------------------------------------------------------*/
// ********************************************************************************************************************

int
main ()
{
  setup ();


  DrawFatouImage (data, IterMax);	// first find Fatou
  SaveArray2PGMFile (data, iWidth, IterMax, ERe, ARe,
		     "DLD/J , name = iWidth_IterMax_ER_AR");

  CheckFatouImage (data);	// 

  DrawDLDImage (data);
  SaveArray2PGMFile (data, iWidth, N, ERe, mi,
		     "DLD/J , name = iWidth_N_ER_mi");


  //test_exterior();

  end ();

  return 0;
}