Fractals/Mathematics/Numerical

Convergence

" ... methods of acceleration of convergence. Suppose you have a slowly convergent series, and want to know its sum (numerically). Just by summing x_1 + x_2 + ... + x_{1000} + ... x_{1000000} + ... + x_{1000000000} + ... you will get the required accuracy after 100 years. If you fit your x_n to c_2/n^2 + c_3/n^3 + (a few more terms), you will get the same accuracy of the sum in 1 second." Andrey (theoretical high energy physicist) ^[1]

Above problem you can find in parabolic Julia sets.

Precision

How to deal with lack of floating point precision:^[2]

implement higher precision arithmetic than your hardware natively supports
- software emulation ( emulating a double with two floats, fixed points numbers , ...)
use algorithms that are more numerically stable (

Escaping test

This test was introduced by by John Milnor ^[3]. See also analysis by Mark Braverman ^[4] and roundoff error by Robert P. Munafo^[5]

Cases

Nonparabolic case

Lets take simple hyperbolic case where parameter c is :

Here repelling fixed point z_f is :

Parabolic case

Lets take simple parabolic case where parameter c is : ^[6]

Here parabolic fixed point z_f is :

Test

Lets take point z of exterior of Julia set but lying near fixed point :

where n is a positive integer

Check how many iterates i needs point z to reach target set ( = escape) :

where ER is Escape Radius

Show relationship between :

n
Last Iteration
type of numbers used for computations ( float, double, long double, extended, arbitrary precision )

Programs

See FractalForum for evaldraw script^[7]

C source code for CPU one thread - hyperbolic case and float numbers . Click on the right to view

/*
 gcc -lm -Wall p.c
./a.out
 */
 
#include <stdio.h>
#include <math.h> // pow
#include <time.h>
 
 
 
 
 
float GiveLastIteration(float Zx, float Zy, float Cx, float Cy, int ER2)
 {
  float Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
 
  float i=0.0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  while  (Zx2+Zy2<ER2)  /* ER2=ER*ER */
  {
   Zy=2*Zx*Zy + Cy;
   Zx=Zx2-Zy2 + Cx;
   Zx2=Zx*Zx;
   Zy2=Zy*Zy;
   i+=1.0;
 
  }
  return i;
 }
 
 
int main()
 
{
  float distance;
 
   int n;
   //  fixed point zf = zfx = zfy*I = 1.0 for c=0
   float zfx = 1.0;
   float zfy = 0.0; 
// z1 = z1x + z1y*I = point of exterior of Julia set but near parabolic fixed point zp
   float z1x;
   float z1y = zfy;
   float cx= 0.0;
   float cy= 0.0;
   // Escape Radius ; it defines target set  = { z: abs(z)>ER}
   // all points z in the target set are escaping to infinity
   float ER=2.0;
   float ER2;
 
   float LastIteration;
 
 
   time_t start,end;
   float dif;
 
 
 
 
    ER2= ER*ER;
 
 
   printf ("Using c with float and Escape Radius = %f \n",  ER);
   for (n =1; n<101; n++)
   {
     time (&start);
     distance = pow(2.0,-n);
     printf("n= %3d distance = %e  = %.20f ",n, distance, distance);
     z1x = zfx + distance;
     LastIteration = GiveLastIteration(z1x,z1y, cx,cy,ER2 );
     time (&end);
     dif = difftime (end,start);
     printf("LI = %3.0f time = %2.0lf seconds\n", LastIteration, dif);
   }
 
 
 
 
  return 0;
}

C source code for CPU one thread - hyperbolic case and double numbers . Click on the right to view

/*
 gcc -lm -Wall p.c
./a.out
 */
#include <stdio.h>
#include <math.h> // pow
#include <time.h>
 
double GiveLastIteration(double Zx, double Zy, double Cx, double Cy, int ER2)
 {
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  double i=0.0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  while  (Zx2+Zy2<ER2)  /* ER2=ER*ER */
  {
   Zy=2*Zx*Zy + Cy;
   Zx=Zx2-Zy2 + Cx;
   Zx2=Zx*Zx;
   Zy2=Zy*Zy;
   i+=1.0;
 
  }
  return i;
 }
 
 
int main()
 
{
  double distance;
 
   int n;
   //  fixed point zf = zfx = zfy*I = 1.0 for c=0
   double zfx = 1.0;
   double zfy = 0.0; 
// z1 = z1x + z1y*I = point of exterior of Julia set but near parabolic fixed point zp
   double z1x;
   double z1y = zfy;
   double cx= 0.0;
   double cy= 0.0;
   // Escape Radius ; it defines target set  = { z: abs(z)>ER}
   // all points z in the target set are escaping to infinity
   double ER=2.0;
   double ER2;
 
   double LastIteration;
 
 
   time_t start,end;
   double dif;
 
 
 
 
    ER2= ER*ER;
 
 
   printf ("Using c with double and Escape Radius = %f \n",  ER);
   for (n =1; n<101; n++)
  {
     time (&start);
     distance = pow(2.0,-n);
     printf("n= %3d distance = %e  = %.20f ",n, distance, distance);
     z1x = zfx + distance;
 
 
     LastIteration = GiveLastIteration(z1x,z1y, cx,cy,ER2 );
     time (&end);
     dif = difftime (end,start);
     printf("LI = %3.0f time = %2.0lf seconds\n", LastIteration, dif);
   }
 
 return 0;
}

C source code for CPU one thread - hyperbolic case and long double numbers . Click on the right to view

/*
 gcc -lm -Wall p.c
./a.out
 */
 
 
 
#include <stdio.h>
#include <math.h> // pow
#include <time.h>
 
 
long double GiveLastIteration(long double Zx0, long double Zy0, long double Cx, long double Cy, int ER2)
 {
  long double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  long double Zx = Zx0;
  long double Zy = Zy0;
  long double i=0.0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  while  (Zx2+Zy2<ER2)  /* ER2=ER*ER */
  {
   Zy=2*Zx*Zy + Cy;
   Zx=Zx2-Zy2 + Cx;
   Zx2=Zx*Zx;
   Zy2=Zy*Zy;
   i+=1.0;
 
  }
  return i;
 }
 
 
int main()
 
{
  long double distance;
 
   int n;
   //  fixed point zp = zpx = zpy*I = 1.0 for c=0
   long double zpx = 1.0;
   long double zpy = 0.0; 
// z1 = z1x + z1y*I = point of exterior of Julia set but near parabolic fixed point zp
   long double z1x;
   long double z1y = zpy;
   long double cx= 0.0;
   long double cy= 0.0;
   // Escape Radius ; it defines target set  = { z: abs(z)>ER}
   // all points z in the target set are escaping to infinity
   long double ER=2.0;
   long double ER2;
 
   long double LastIteration;
 
 
   time_t start,end;
   long double dif;
 
 
 
    ER2= ER*ER;
 
 
 printf ("Using c with long double and Escape Radius = %Lf \n",  ER);
   for (n =1; n<101; n++)
   {
     time (&start);
     distance = pow(2.0,-n);
     z1x = zpx + distance;
     LastIteration = GiveLastIteration(z1x,z1y, cx,cy,ER2 );
     time (&end);
     dif = difftime (end,start);
     printf("n= %3d distance = %Le  = %.10Lf LI = %10.0Lf log2(LI) = %3.0Lf time = %2.0Lf seconds\a\n",n, distance, distance, LastIteration, log2l(LastIteration), dif);
   }
 
 
 
 
  return 0;
}

C source code for CPU one thread - parabolic case and float numbers . Click on the right to view

/*
 gcc -lm -Wall p.c
./a.out
 */
#include <stdio.h>
#include <math.h> // pow
#include <time.h>
 
 
float GiveLastIteration(float Zx0, float Zy0, float Cx, float Cy, int ER2)
 {
  float Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  float Zx = Zx0;
  float Zy = Zy0;
  float i=0.0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  while  (Zx2+Zy2<ER2)  /* ER2=ER*ER */
  {
   Zy=2*Zx*Zy + Cy;
   Zx=Zx2-Zy2 + Cx;
   Zx2=Zx*Zx;
   Zy2=Zy*Zy;
   i+=1.0;
  }
  return i;
 }
 
 
int main()
 
{
  float distance;
 
   int n;
   // parabolic fixed point zp = zpx = zpy*I = 0.5 for c=1/4
   float zpx = 0.5;
   float zpy = 0.0; 
// z1 = z1x + z1y*I = point of exterior of Julia set but near parabolic fixed point zp
   float z1x;
   float z1y = zpy;
   float cx= 0.25;
   float cy= 0.0;
   // Escape Radius ; it defines target set  = { z: abs(z)>ER}
   // all points z in the target set are escaping to infinity
   float ER=2.0;
   float ER2;
 
   float LastIteration;
 
 
   time_t start,end;
   float dif;
 
 
 
    ER2= ER*ER;
 
 
   printf ("Using c with float and Escape Radius = %f \n",  ER);
   for (n =1; n<101; n++)
   {
     time (&start);
     distance = pow(2.0,-n);
     z1x = zpx + distance;
     LastIteration = GiveLastIteration(z1x,z1y, cx,cy,ER2 );
     time (&end);
     dif = difftime (end,start);
     printf("n= %3d distance = %e  = %.10f LI = %10.0f log2(LI) = %3.0f time = %2.0lf seconds\n",n, distance, distance, LastIteration, log2(LastIteration), dif);
   }
 
 
 
 
  return 0;
}

C source code for CPU one thread - parabolic case and double numbers . Click on the right to view

/*
 gcc -lm -Wall p.c
./a.out
 */
 
 
#include <stdio.h>
#include <math.h> // pow
#include <time.h>
 
 
double GiveLastIteration(double Zx0, double Zy0, double Cx, double Cy, int ER2)
 {
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  double Zx = Zx0;
  double Zy = Zy0;
  double i=0.0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  while  (Zx2+Zy2<ER2)  /* ER2=ER*ER */
  {
   Zy=2*Zx*Zy + Cy;
   Zx=Zx2-Zy2 + Cx;
   Zx2=Zx*Zx;
   Zy2=Zy*Zy;
   i+=1.0;
  }
  return i;
 }
 
 
 
double max(double a , double b)
{
  if (a>b )
    return a;
    else return b;
 
}
 
 
int main()
 
{  double distance;
   int n;
   // parabolic fixed point zp = zpx = zpy*I = 0.5 for c=1/4
   double zpx = 0.5;
   double zpy = 0.0; 
// z1 = z1x + z1y*I = point of exterior of Julia set but near parabolic fixed point zp
   double z1x;
   double z1y = zpy;
   double cx= 0.25;
   double cy= 0.0;
   // Escape Radius ; it defines target set  = { z: abs(z)>ER}
   // all points z in the target set are escaping to infinity
   double ER=2.0;
   double ER2;
 
   double LastIteration;
 
 
   time_t start,end;
   double dif;
 
 
 
 
 
    ER2= ER*ER;
 
  for (n =1; n<101; n++)
   {
     time (&start);
     distance = pow(2.0,-n);
     z1x = zpx + distance;
     LastIteration = GiveLastIteration(z1x,z1y, cx,cy,ER2 );
     time (&end);
     dif = difftime (end,start);
     printf("n= %3d distance = %e  = %.10f LI = %10.0f log2(LI) = %3.0f time = %2.0lf seconds\n",    
             n, distance, distance, LastIteration, log2(LastIteration), dif);
    }
  return 0;
}

C source code for CPU one thread - parabolic case and long double numbers . Click on the right to view

/*
 gcc -lm -Wall p.c
./a.out
 */
 
 
 
#include <stdio.h>
#include <math.h> // pow
#include <time.h>
 
 
long double GiveLastIteration(long double Zx0, long double Zy0, long double Cx, long double Cy, int ER2)
 {
  long double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  long double Zx = Zx0;
  long double Zy = Zy0;
 
  long double i=0.0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  while  (Zx2+Zy2<ER2)  /* ER2=ER*ER */
  {
   Zy=2*Zx*Zy + Cy;
   Zx=Zx2-Zy2 + Cx;
   Zx2=Zx*Zx;
   Zy2=Zy*Zy;
   i+=1.0;
  }
  return i;
 }
 
 
int main()
 
{
  long double distance;
 
   int n;
   // parabolic fixed point zp = zpx = zpy*I = 0.5 for c=1/4
   long double zpx = 0.5;
   long double zpy = 0.0; 
// z1 = z1x + z1y*I = point of exterior of Julia set but near parabolic fixed point zp
   long double z1x;
   long double z1y = zpy;
   long double cx= 0.25;
   long double cy= 0.0;
   // Escape Radius ; it defines target set  = { z: abs(z)>ER}
   // all points z in the target set are escaping to infinity
   long double ER=2.0;
   long double ER2;
 
   long double LastIteration;
 
 
   time_t start,end;
   long double dif;
 
    ER2= ER*ER;
 
  for (n =1; n<101; n++)
   {
     time (&start);
     distance = pow(2.0,-n);
     z1x = zpx + distance;
     LastIteration = GiveLastIteration(z1x,z1y, cx,cy,ER2 );
     time (&end);
     dif = difftime (end,start);
     printf("n= %3d distance = %Le  = %.10Lf LI = %10.0Lf log2(LI) = %3.0Lf time = %2.0Lf seconds\a\n",n, distance, distance, LastIteration, log2l(LastIteration), dif);
   }
 
 return 0;
}

C source code for CPU one thread - parabolic case and GSL library . Click on the right to view

/*
http://www.gnu.org/software/gsl/
 
 gcc -L/usr/lib -lgsl -lgslcblas -lm -Wall pgsl.c 
./a.out
 
or 
 
GSL_IEEE_MODE="extended-precision" ./a.out
 
*/
 
#include <stdio.h>
#include <math.h> // pow
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <time.h>
 
 
long double GiveLastIteration(gsl_complex z, gsl_complex c, int ER2)
 {
 
  long double i=0.0;
 
  while  (gsl_complex_abs2 (z)< ER2)  /* ER2=ER*ER */
  {
   z = gsl_complex_add(c, gsl_complex_mul(z,z));
   i+=1.0;
 
  }
  return i;
 }
 
 
int main()
 
{
  //long double distance;
 
   int n;
   // parabolic fixed point zp = zpx = zpy*I = 0.5 for c=1/4
 
  gsl_complex zp = gsl_complex_rect(0.5, 0.0);
   //long double zpx = 0.5;
   //long double zpy = 0.0; 
// z1 = z1x + z1y*I = point of exterior of Julia set but near parabolic fixed point zp
   gsl_complex z1;
   //long double z1x;
   //long double z1y = zpy;
   //
   gsl_complex c = gsl_complex_rect(0.25, 0.0);
   //long double cx= 0.25;
   //long double cy= 0.0;
   // Escape Radius ; it defines target set  = { z: abs(z)>ER}
   // all points z in the target set are escaping to infinity
   long double ER=2.0;
   long double ER2;
 
   long double LastIteration;
 
 
   time_t start,end;
   long double dif;
 
 
 
    ER2= ER*ER;
 
  //n=27;
 
   for (n =1; n<101; n++)
   {
     time (&start);
     z1 = gsl_complex_rect(GSL_REAL(zp) + pow(2.0,-n), 0.0);
     LastIteration = GiveLastIteration(z1, c,ER2 );
     time (&end);
     dif = difftime (end,start);
     printf("n= %3d LI = %10.0Lf log2(LI) = %3.0Lf time = %2.0Lf seconds\a\n",n, LastIteration, log2l(LastIteration), dif);
   }
 
 
 
 
  return 0;
}

C source code for CPU one thread - parabolic case and mpfr library . Click on the right to view

/*
 gcc  -lm -lmpfr -lgmp -Wall -O2 pmpfr.c
./a.out
 
*/
 
#include <stdio.h>
#include <math.h> // log2
#include <gmp.h>
#include <mpfr.h>
#include <time.h>
 
   // MPFR general settings
  mpfr_prec_t precision = 12; // the number of bits used to represent the significand of a floating-point number
  mpfr_rnd_t RoundMode = MPFR_RNDD; // the way to round the result of a floating-point operation, round toward minus infinity (roundTowardNegative in IEEE 754-2008),
 
 
 
 
// z(i+1) = fc(zi) = zi^2 + c
// iteration of z under complex quadratic polynomial
double GiveLastIteration_mpfr(
                     mpfr_t zx,
                     mpfr_t zy,
                     mpfr_t cx,
                     mpfr_t cy,    
                     mpfr_t ER2)
{
 
    double i=0.0;
    mpfr_t temp; // temporary variable 
    mpfr_t zx2; // zx^2
    mpfr_t zy2; // zy^2
 
    // initializes MPFR variables with default_prec
   mpfr_inits(temp, zx2, zy2, (mpfr_ptr) 0 ); 
 
    mpfr_mul(   zy2,   zy,  zy,    RoundMode); /* zy2 = zy * zy; */
    mpfr_mul(   zx2,   zx,  zx,    RoundMode); /* zx2 = zx * zx; */
    mpfr_add(   temp,   zx2,  zy2,   RoundMode); //  temp = zx2 + zy2
 
    while (mpfr_greater_p(ER2, temp)) //  while  (Zx2+Zy2<ER2
    {
        /* zy = 2.0*zx*zy + cy; */
        mpfr_mul(   temp,     zx,    zy,    RoundMode); // temp = zx*zy
        mpfr_mul_ui(temp,     temp,     2.0,      RoundMode); // temp = 2 * temp = 2*zx*zy
        mpfr_add(  zy,    temp,     cy,   RoundMode);
        /* zx = zx^2 - zy^2 + cx; */
        mpfr_sub(   temp,     zx2,   zy2,   RoundMode); // temp = zx^2 - zy^2
        mpfr_add(   zx,    temp,     cx,   RoundMode); //  temp = temp + cx
        //
        mpfr_mul(   zy2,   zy,  zy,    RoundMode); /* zy2 = zy * zy; */
        mpfr_mul(   zx2,   zx,  zx,    RoundMode); /* zx2 = zx * zx; */
        mpfr_add(   temp,   zx2,  zy2,   RoundMode); //  temp = zx2 + zy2
        i+=1.0;
    }
 
 
     mpfr_clears (temp, zx2, zy2, (mpfr_ptr) 0 );
    return i;
}
 
 
 
 
 
int main (void)
{
 
 
  // declares variables
 
 
   double LastIteration;
   time_t start,end;
   double dif;
   int n;
   mpfr_t zpx, zpy;  // parabolic fixed point zp = zpx + zpy*I = 0.5 for c=1/4
   mpfr_t zx, zy;  //  point z = zx + zy*I 
   mpfr_t distance; // between zp and z
   mpfr_t cx, cy;  //  point c = cx + cy*
 
   mpfr_t base;
   mpfr_t ER ; // Escape Radius ; it defines target set  = { z: abs(z)>ER}; all points z in the target set are escaping to infinity
   mpfr_t ER2; // ER2= ER * ER
 
 
  mpfr_set_default_prec(precision);
  // initializes MPFR variables with default_prec
  mpfr_inits(ER, ER2, zpx, zpy, zx, zy, distance, cx, cy, base, (mpfr_ptr) 0 ); 
 
  // Assignment Functions
  mpfr_set_ld (zpx, 0.5, RoundMode);
  mpfr_set_ld (zpy, 0.0, RoundMode);
  mpfr_set_ld (zy, 0.0, RoundMode);
  mpfr_set_ld (cx, 0.25, RoundMode);
  mpfr_set_ld (cy, 0.0, RoundMode);
  mpfr_set_ld (base, 2.0, RoundMode);
  mpfr_set_ld (ER, 2.0, RoundMode);
 
 
  // computations
  mpfr_mul(ER2,  ER,     ER,      RoundMode); // 
  mpfr_printf ("Using MPFR-%s with GMP-%s with precision = %u bits and Escape Radius = %Rf \n", mpfr_version, gmp_version, (unsigned int) precision, ER);
  for (n = 1; n <= 60; n++)
    {
     time (&start);
     mpfr_pow_si(distance,base, - n, RoundMode);
     mpfr_add(zx, zpx, distance, RoundMode);
     LastIteration = GiveLastIteration_mpfr(zx, zy, cx, cy, ER2);
     time (&end);
     dif = difftime(end,start);
     mpfr_printf ("n = %3d distance = %.10Re LI = %10.0f log2(LI) = %3.0f; time = %2.0f seconds \n", n,  distance, LastIteration, log2(LastIteration), dif);
    }
 
 
 
  // free the space used by the MPFR variables
  mpfr_clears (ER, ER2, zpx, zpy, zx, zy, distance, cx, cy, base, (mpfr_ptr) 0);
  mpfr_free_cache ();  /* free the cache for constants like pi */
 
 
return 0;
}

Results

Maximal n for which program does not fall into fixed point (columns = number types (precision of significand); rows = case)
	float (24)	double (53)	long double (64)	MPFR (80)	MPFR (100)
hyperbolic	23	52	63	79	99
parabolic	12	26	32

The results for standard C types ( float, double, long double) and MPFR precision are the same for the same precision

Relation between : Last iteration, n ( in columns) and case ( in rows)
	1	2	3	4	5	24	53	64	80	100
hyperbolic	1	2	3	4	5	24	53	64	80	100
parabolic	3	5	10	19	35	16 778 821

Relation between number of iterations and time of computation in hyperbolic case :

Using MPFR-3.0.0-p8 with GMP-4.3.2 with precision = 128 bits and Escape Radius = 2.000000 
n =   1 distance = 5.0000000000e-01 LI =          1 log2(LI) =   0; time =  0 seconds 
n =   2 distance = 2.5000000000e-01 LI =          2 log2(LI) =   1; time =  0 seconds 
n =   3 distance = 1.2500000000e-01 LI =          3 log2(LI) =   2; time =  0 seconds 
n =   4 distance = 6.2500000000e-02 LI =          4 log2(LI) =   2; time =  0 seconds 
n =   5 distance = 3.1250000000e-02 LI =          5 log2(LI) =   2; time =  0 seconds 
n =   6 distance = 1.5625000000e-02 LI =          6 log2(LI) =   3; time =  0 seconds 
n =   7 distance = 7.8125000000e-03 LI =          7 log2(LI) =   3; time =  0 seconds 
n =   8 distance = 3.9062500000e-03 LI =          8 log2(LI) =   3; time =  0 seconds 
n =   9 distance = 1.9531250000e-03 LI =          9 log2(LI) =   3; time =  0 seconds 
n =  10 distance = 9.7656250000e-04 LI =         10 log2(LI) =   3; time =  0 seconds 
n =  11 distance = 4.8828125000e-04 LI =         11 log2(LI) =   3; time =  0 seconds 
n =  12 distance = 2.4414062500e-04 LI =         12 log2(LI) =   4; time =  0 seconds 
n =  13 distance = 1.2207031250e-04 LI =         13 log2(LI) =   4; time =  0 seconds 
n =  14 distance = 6.1035156250e-05 LI =         14 log2(LI) =   4; time =  0 seconds 
n =  15 distance = 3.0517578125e-05 LI =         15 log2(LI) =   4; time =  0 seconds 
n =  16 distance = 1.5258789062e-05 LI =         16 log2(LI) =   4; time =  0 seconds 
n =  17 distance = 7.6293945312e-06 LI =         17 log2(LI) =   4; time =  0 seconds 
n =  18 distance = 3.8146972656e-06 LI =         18 log2(LI) =   4; time =  0 seconds 
n =  19 distance = 1.9073486328e-06 LI =         19 log2(LI) =   4; time =  0 seconds 
n =  20 distance = 9.5367431641e-07 LI =         20 log2(LI) =   4; time =  0 seconds 
n =  21 distance = 4.7683715820e-07 LI =         21 log2(LI) =   4; time =  0 seconds 
n =  22 distance = 2.3841857910e-07 LI =         22 log2(LI) =   4; time =  0 seconds 
n =  23 distance = 1.1920928955e-07 LI =         23 log2(LI) =   5; time =  0 seconds 
n =  24 distance = 5.9604644775e-08 LI =         24 log2(LI) =   5; time =  0 seconds 
n =  25 distance = 2.9802322388e-08 LI =         25 log2(LI) =   5; time =  0 seconds 
n =  26 distance = 1.4901161194e-08 LI =         26 log2(LI) =   5; time =  0 seconds 
n =  27 distance = 7.4505805969e-09 LI =         27 log2(LI) =   5; time =  0 seconds 
n =  28 distance = 3.7252902985e-09 LI =         28 log2(LI) =   5; time =  0 seconds 
n =  29 distance = 1.8626451492e-09 LI =         29 log2(LI) =   5; time =  0 seconds 
n =  30 distance = 9.3132257462e-10 LI =         30 log2(LI) =   5; time =  0 seconds 
n =  31 distance = 4.6566128731e-10 LI =         31 log2(LI) =   5; time =  0 seconds 
n =  32 distance = 2.3283064365e-10 LI =         32 log2(LI) =   5; time =  0 seconds 
n =  33 distance = 1.1641532183e-10 LI =         33 log2(LI) =   5; time =  0 seconds 
n =  34 distance = 5.8207660913e-11 LI =         34 log2(LI) =   5; time =  0 seconds 
n =  35 distance = 2.9103830457e-11 LI =         35 log2(LI) =   5; time =  0 seconds 
n =  36 distance = 1.4551915228e-11 LI =         36 log2(LI) =   5; time =  0 seconds 
n =  37 distance = 7.2759576142e-12 LI =         37 log2(LI) =   5; time =  0 seconds 
n =  38 distance = 3.6379788071e-12 LI =         38 log2(LI) =   5; time =  0 seconds 
n =  39 distance = 1.8189894035e-12 LI =         39 log2(LI) =   5; time =  0 seconds 
n =  40 distance = 9.0949470177e-13 LI =         40 log2(LI) =   5; time =  0 seconds 
n =  41 distance = 4.5474735089e-13 LI =         41 log2(LI) =   5; time =  0 seconds 
n =  42 distance = 2.2737367544e-13 LI =         42 log2(LI) =   5; time =  0 seconds 
n =  43 distance = 1.1368683772e-13 LI =         43 log2(LI) =   5; time =  0 seconds 
n =  44 distance = 5.6843418861e-14 LI =         44 log2(LI) =   5; time =  0 seconds 
n =  45 distance = 2.8421709430e-14 LI =         45 log2(LI) =   5; time =  0 seconds 
n =  46 distance = 1.4210854715e-14 LI =         46 log2(LI) =   6; time =  0 seconds 
n =  47 distance = 7.1054273576e-15 LI =         47 log2(LI) =   6; time =  0 seconds 
n =  48 distance = 3.5527136788e-15 LI =         48 log2(LI) =   6; time =  0 seconds 
n =  49 distance = 1.7763568394e-15 LI =         49 log2(LI) =   6; time =  0 seconds 
n =  50 distance = 8.8817841970e-16 LI =         50 log2(LI) =   6; time =  0 seconds 
n =  51 distance = 4.4408920985e-16 LI =         51 log2(LI) =   6; time =  0 seconds 
n =  52 distance = 2.2204460493e-16 LI =         52 log2(LI) =   6; time =  0 seconds 
n =  53 distance = 1.1102230246e-16 LI =         53 log2(LI) =   6; time =  0 seconds 
n =  54 distance = 5.5511151231e-17 LI =         54 log2(LI) =   6; time =  0 seconds 
n =  55 distance = 2.7755575616e-17 LI =         55 log2(LI) =   6; time =  0 seconds 
n =  56 distance = 1.3877787808e-17 LI =         56 log2(LI) =   6; time =  0 seconds 
n =  57 distance = 6.9388939039e-18 LI =         57 log2(LI) =   6; time =  0 seconds 
n =  58 distance = 3.4694469520e-18 LI =         58 log2(LI) =   6; time =  0 seconds 
n =  59 distance = 1.7347234760e-18 LI =         59 log2(LI) =   6; time =  0 seconds 
n =  60 distance = 8.6736173799e-19 LI =         60 log2(LI) =   6; time =  0 seconds 
n =  61 distance = 4.3368086899e-19 LI =         61 log2(LI) =   6; time =  0 seconds 
n =  62 distance = 2.1684043450e-19 LI =         62 log2(LI) =   6; time =  0 seconds 
n =  63 distance = 1.0842021725e-19 LI =         63 log2(LI) =   6; time =  0 seconds 
n =  64 distance = 5.4210108624e-20 LI =         64 log2(LI) =   6; time =  0 seconds 
n =  65 distance = 2.7105054312e-20 LI =         65 log2(LI) =   6; time =  0 seconds 
n =  66 distance = 1.3552527156e-20 LI =         66 log2(LI) =   6; time =  0 seconds 
n =  67 distance = 6.7762635780e-21 LI =         67 log2(LI) =   6; time =  0 seconds 
n =  68 distance = 3.3881317890e-21 LI =         68 log2(LI) =   6; time =  0 seconds 
n =  69 distance = 1.6940658945e-21 LI =         69 log2(LI) =   6; time =  0 seconds 
n =  70 distance = 8.4703294725e-22 LI =         70 log2(LI) =   6; time =  0 seconds 
n =  71 distance = 4.2351647363e-22 LI =         71 log2(LI) =   6; time =  0 seconds 
n =  72 distance = 2.1175823681e-22 LI =         72 log2(LI) =   6; time =  0 seconds 
n =  73 distance = 1.0587911841e-22 LI =         73 log2(LI) =   6; time =  0 seconds 
n =  74 distance = 5.2939559203e-23 LI =         74 log2(LI) =   6; time =  0 seconds 
n =  75 distance = 2.6469779602e-23 LI =         75 log2(LI) =   6; time =  0 seconds 
n =  76 distance = 1.3234889801e-23 LI =         76 log2(LI) =   6; time =  0 seconds 
n =  77 distance = 6.6174449004e-24 LI =         77 log2(LI) =   6; time =  0 seconds 
n =  78 distance = 3.3087224502e-24 LI =         78 log2(LI) =   6; time =  0 seconds 
n =  79 distance = 1.6543612251e-24 LI =         79 log2(LI) =   6; time =  0 seconds 
n =  80 distance = 8.2718061255e-25 LI =         80 log2(LI) =   6; time =  0 seconds 
n =  81 distance = 4.1359030628e-25 LI =         81 log2(LI) =   6; time =  0 seconds 
n =  82 distance = 2.0679515314e-25 LI =         82 log2(LI) =   6; time =  0 seconds 
n =  83 distance = 1.0339757657e-25 LI =         83 log2(LI) =   6; time =  0 seconds 
n =  84 distance = 5.1698788285e-26 LI =         84 log2(LI) =   6; time =  0 seconds 
n =  85 distance = 2.5849394142e-26 LI =         85 log2(LI) =   6; time =  0 seconds 
n =  86 distance = 1.2924697071e-26 LI =         86 log2(LI) =   6; time =  0 seconds 
n =  87 distance = 6.4623485356e-27 LI =         87 log2(LI) =   6; time =  0 seconds 
n =  88 distance = 3.2311742678e-27 LI =         88 log2(LI) =   6; time =  0 seconds 
n =  89 distance = 1.6155871339e-27 LI =         89 log2(LI) =   6; time =  0 seconds 
n =  90 distance = 8.0779356695e-28 LI =         90 log2(LI) =   6; time =  0 seconds 
n =  91 distance = 4.0389678347e-28 LI =         91 log2(LI) =   7; time =  0 seconds 
n =  92 distance = 2.0194839174e-28 LI =         92 log2(LI) =   7; time =  0 seconds 
n =  93 distance = 1.0097419587e-28 LI =         93 log2(LI) =   7; time =  0 seconds 
n =  94 distance = 5.0487097934e-29 LI =         94 log2(LI) =   7; time =  0 seconds 
n =  95 distance = 2.5243548967e-29 LI =         95 log2(LI) =   7; time =  0 seconds 
n =  96 distance = 1.2621774484e-29 LI =         96 log2(LI) =   7; time =  0 seconds 
n =  97 distance = 6.3108872418e-30 LI =         97 log2(LI) =   7; time =  0 seconds 
n =  98 distance = 3.1554436209e-30 LI =         98 log2(LI) =   7; time =  0 seconds 
n =  99 distance = 1.5777218104e-30 LI =         99 log2(LI) =   7; time =  0 seconds 
n = 100 distance = 7.8886090522e-31 LI =        100 log2(LI) =   7; time =  0 seconds 
n = 101 distance = 3.9443045261e-31 LI =        101 log2(LI) =   7; time =  0 seconds 
n = 102 distance = 1.9721522631e-31 LI =        102 log2(LI) =   7; time =  0 seconds 
n = 103 distance = 9.8607613153e-32 LI =        103 log2(LI) =   7; time =  0 seconds 
n = 104 distance = 4.9303806576e-32 LI =        104 log2(LI) =   7; time =  0 seconds 
n = 105 distance = 2.4651903288e-32 LI =        105 log2(LI) =   7; time =  0 seconds 
n = 106 distance = 1.2325951644e-32 LI =        106 log2(LI) =   7; time =  0 seconds 
n = 107 distance = 6.1629758220e-33 LI =        107 log2(LI) =   7; time =  0 seconds 
n = 108 distance = 3.0814879110e-33 LI =        108 log2(LI) =   7; time =  0 seconds 
n = 109 distance = 1.5407439555e-33 LI =        109 log2(LI) =   7; time =  0 seconds 
n = 110 distance = 7.7037197775e-34 LI =        110 log2(LI) =   7; time =  0 seconds 
n = 111 distance = 3.8518598888e-34 LI =        111 log2(LI) =   7; time =  0 seconds 
n = 112 distance = 1.9259299444e-34 LI =        112 log2(LI) =   7; time =  0 seconds 
n = 113 distance = 9.6296497219e-35 LI =        113 log2(LI) =   7; time =  0 seconds 
n = 114 distance = 4.8148248610e-35 LI =        114 log2(LI) =   7; time =  0 seconds 
n = 115 distance = 2.4074124305e-35 LI =        115 log2(LI) =   7; time =  0 seconds 
n = 116 distance = 1.2037062152e-35 LI =        116 log2(LI) =   7; time =  0 seconds 
n = 117 distance = 6.0185310762e-36 LI =        117 log2(LI) =   7; time =  0 seconds 
n = 118 distance = 3.0092655381e-36 LI =        118 log2(LI) =   7; time =  0 seconds 
n = 119 distance = 1.5046327691e-36 LI =        119 log2(LI) =   7; time =  0 seconds 
n = 120 distance = 7.5231638453e-37 LI =        120 log2(LI) =   7; time =  0 seconds 
n = 121 distance = 3.7615819226e-37 LI =        121 log2(LI) =   7; time =  0 seconds 
n = 122 distance = 1.8807909613e-37 LI =        122 log2(LI) =   7; time =  0 seconds 
n = 123 distance = 9.4039548066e-38 LI =        123 log2(LI) =   7; time =  0 seconds 
n = 124 distance = 4.7019774033e-38 LI =        124 log2(LI) =   7; time =  0 seconds 
n = 125 distance = 2.3509887016e-38 LI =        125 log2(LI) =   7; time =  0 seconds 
n = 126 distance = 1.1754943508e-38 LI =        126 log2(LI) =   7; time =  0 seconds 
n = 127 distance = 5.8774717541e-39 LI =        127 log2(LI) =   7; time =  0 seconds

Parabolic case :

Using MPFR-3.0.0-p8 with GMP-4.3.2 with precision = 100 bits and Escape Radius = 2.000000 
n =   1 distance = 5.0000000000e-01 LI =           3 log2(LI) =   2; time =     0 seconds 
n =   2 distance = 2.5000000000e-01 LI =           5 log2(LI) =   2; time =     0 seconds 
n =   3 distance = 1.2500000000e-01 LI =          10 log2(LI) =   3; time =     0 seconds 
n =   4 distance = 6.2500000000e-02 LI =          19 log2(LI) =   4; time =     0 seconds 
n =   5 distance = 3.1250000000e-02 LI =          35 log2(LI) =   5; time =     0 seconds 
n =   6 distance = 1.5625000000e-02 LI =          68 log2(LI) =   6; time =     0 seconds 
n =   7 distance = 7.8125000000e-03 LI =         133 log2(LI) =   7; time =     0 seconds 
n =   8 distance = 3.9062500000e-03 LI =         261 log2(LI) =   8; time =     0 seconds 
n =   9 distance = 1.9531250000e-03 LI =         518 log2(LI) =   9; time =     0 seconds 
n =  10 distance = 9.7656250000e-04 LI =        1031 log2(LI) =  10; time =     0 seconds 
n =  11 distance = 4.8828125000e-04 LI =        2055 log2(LI) =  11; time =     0 seconds 
n =  12 distance = 2.4414062500e-04 LI =        4104 log2(LI) =  12; time =     0 seconds 
n =  13 distance = 1.2207031250e-04 LI =        8201 log2(LI) =  13; time =     0 seconds 
n =  14 distance = 6.1035156250e-05 LI =       16394 log2(LI) =  14; time =     0 seconds 
n =  15 distance = 3.0517578125e-05 LI =       32778 log2(LI) =  15; time =     0 seconds 
n =  16 distance = 1.5258789062e-05 LI =       65547 log2(LI) =  16; time =     0 seconds 
n =  17 distance = 7.6293945312e-06 LI =      131084 log2(LI) =  17; time =     0 seconds 
n =  18 distance = 3.8146972656e-06 LI =      262156 log2(LI) =  18; time =     0 seconds 
n =  19 distance = 1.9073486328e-06 LI =      524301 log2(LI) =  19; time =     0 seconds 
n =  20 distance = 9.5367431641e-07 LI =     1048590 log2(LI) =  20; time =     1 seconds 
n =  21 distance = 4.7683715820e-07 LI =     2097166 log2(LI) =  21; time =     0 seconds 
n =  22 distance = 2.3841857910e-07 LI =     4194319 log2(LI) =  22; time =     2 seconds 
n =  23 distance = 1.1920928955e-07 LI =     8388624 log2(LI) =  23; time =     2 seconds 
n =  24 distance = 5.9604644775e-08 LI =    16777232 log2(LI) =  24; time =     6 seconds 
n =  25 distance = 2.9802322388e-08 LI =    33554449 log2(LI) =  25; time =    11 seconds 
n =  26 distance = 1.4901161194e-08 LI =    67108882 log2(LI) =  26; time =    21 seconds 
n =  27 distance = 7.4505805969e-09 LI =   134217747 log2(LI) =  27; time =    42 seconds 
n =  28 distance = 3.7252902985e-09 LI =   268435475 log2(LI) =  28; time =    87 seconds 
n =  29 distance = 1.8626451492e-09 LI =   536870932 log2(LI) =  29; time =   175 seconds 
n =  30 distance = 9.3132257462e-10 LI =  1073741845 log2(LI) =  30; time =   351 seconds 
n =  31 distance = 4.6566128731e-10 LI =  2147483669 log2(LI) =  31; time =   698 seconds 
n =  32 distance = 2.3283064365e-10 LI =  4294967318 log2(LI) =  32; time =  1386 seconds 
n =  33 distance = 1.1641532183e-10 LI =  8589934615 log2(LI) =  33; time =  2714 seconds 
n =  34 distance = 5.8207660913e-11 LI = 17179869207 log2(LI) =  34; time =  5595 seconds 
n =  35 distance = 2.9103830457e-11 LI = 34359738392 log2(LI) =  35; time = 11175 seconds 
n =  36 distance = 1.4551915228e-11 LI = 68719476762 log2(LI) =  36; time = 22081 seconds

Analysis

Maximal n in hyperbolic case $n = n_h$ is allmost the same as the precision of the significand ^[8]

In parabolic case maximal $n = n_p$ is

Last Iteration ( escape time = iteration fro which abs(zn) > ER ) is : in hyperbolic case equal to n :

in parabolic case equal to 2^n :

Time of computations is proportional to number of iterations. In hyperbolic case is is short. In parabolic case grows quickly as number of iterations.

Checking one point if escapes in parabolic case :

for n = 34 take about one hour ( 5 595 seconds )
for n = 40 take about one day
for n = 45 take about one month
for n = 50 take about one year

Q&A

Why programs fails ?

Cancellation of significant digits^[9] and loss of significance ( LOS).^[10]

The program fails because of limited precision of used number types. Addition of big (zp) and small number (distance) gives number which has more decimal digits then can be saved ( floating point type has only 7 decimal digits). Some of the most right digits are cancelled and iteration goes into an infinite loop.

For example : when using floating point in parabolic case lets take

float cx = 0.25;
float Zpx = 0.5;
float Zx ;
float distance;
float Zx2;
float n = 13;

distance = pow(2.0,-n); // = 1.2207031250e-04 = 0,00012207

It is greater then machine epsilon^[11] :

distance > FLT_EPSILON // = pow(2, -24) = 5.96e-08 = 0,00000006

so this addition still works :

Zx = Zpx + distance; // adding big and small number gives 0,50012207

After multiplication it gives :

Zx2 = Zx*Zx; // = 0,250122

next step is addition. Because floating point format saves only 7 decimal digits it is truncated to :

Zx = Zx2 + cx; // = 0,500122 = Zp + (distance/2)

Here relative error is to big and

d2= 0.0000000149 // distance*distance

is smaller then FLT_EPSILON/2.0 = 0.0000000596;

Solution : increase precision !

#include <stdio.h>
#include <math.h> // pow
#include <float.h> // FLT_EPSILON
#include <time.h>
#include <fenv.h> //fegetround()
 
 
 
 
 
int main()
 
{
 
 
 
 
   float cx= 0.25;
   // Escape Radius ; it defines target set  = { z: abs(z)>ER}
   // all points z in the target set are escaping to infinity
   float ER=2.0;
   float ER2;
 
 
 
 
   time_t start,end;
   float dif;
 
 
 
    ER2= ER*ER;
 
 
 
 
   float Zpx = 0.5; 
   float Zx ;// bad value = 0.5002; good value = 0.5004
   float Zx2; /* Zx2=Zx*Zx */
   float i=0.0;
   float d ; // distance between Zpx=1/2 and zx
   float d2; // d2=d*d;
   int n = 13;
 
   d = pow(2.0,-n);
   Zx = Zpx + d;
   d2=d*d;
 
 
    time (&start);
    Zx2=Zpx*Zpx + 2.0*d*Zpx +d2;
    printf ("Using c with float and Escape Radius = %f \n",  ER);
    printf("Round mode is    = %d \n",fegetround());
    printf("i= %3.0f; Zx = %f;  Zx2 = %10.8f ;  d = %f ;  d2 = %.10f\n", i,Zx, Zx2, d , d2);
    if (d2<(FLT_EPSILON/2.0)) 
       {printf("error : relative error to big and d2= %.10f is smaller then FLT_EPSILON/2.0 = %.10f; increase precision ! \a\n", 
                d2, FLT_EPSILON/2.0);    
         return 1;}
 
    while  (Zx2<ER2)  /* ER2=ER*ER */
    {
 
     Zx=Zx2 + cx;
     d = Zx-Zpx;
     d2=d*d; 
     Zx2= 0.25+ d +d2; // zx2 = zx*zx = (zp+d)*(zp+d)= zp2 +2*d*zp +d2 = 2.25+d +d2
     i+=1.0;
     //printf("i= %3.0f; Zx = %f;  Zx2 = %10.8f ;  d = %f ;  d2 = %f \n", i,Zx, Zx2, d,d2);
    }
 
   time (&end);
   dif = difftime (end,start);
   printf("n = %d; distance = %3f; LI = %10.0f log2(LI) = %3.0f time = %2.0lf seconds\n",n, d, i, log2(i), dif);
 
 
 
 
  return 0;
}