numerical problem with bailout test

[Adam Majewski]

I have two numbers z and zf.

In first case ( easy ) :
double zf = 1.0; // fixed point for fc(z) = z*z
double z = zf + pow(2.0,-53); // point of exterior of Julia set
After comparing its memory representations I see that these numbers  z and zf are the same for computer because z<nextafter (zf , 2.0);
So z does not escapes to infinity under f(z) = z*z but falls into fixed point zf.

In second case :
double zf = 0.5; // parabolic fixed point
double z = zf + pow(2.0,-27); // point of exterior of Julia set
When I iterate z using fc(z) = z*z + 0.25 then z is not escaping from zf.
Why it is so ?
Here z>nextafter (zf , 2.0);

What happens in second case ?
Can you check second case using your software ?

Theory :
http://arxiv.org/abs/math/0505036

TIA

[Adam Majewski]

I do not know what is the name of this numerical error ( underflow ?) but changing double to long double lets go futher :

n=   1 distance = 5.000000e-01  = 0.5000000000 LI =          3 log2(LI) =   2 time =  0 seconds
n=   2 distance = 2.500000e-01  = 0.2500000000 LI =          5 log2(LI) =   2 time =  0 seconds
n=   3 distance = 1.250000e-01  = 0.1250000000 LI =         10 log2(LI) =   3 time =  0 seconds
n=   4 distance = 6.250000e-02  = 0.0625000000 LI =         19 log2(LI) =   4 time =  0 seconds
n=   5 distance = 3.125000e-02  = 0.0312500000 LI =         35 log2(LI) =   5 time =  0 seconds
n=   6 distance = 1.562500e-02  = 0.0156250000 LI =         68 log2(LI) =   6 time =  0 seconds
n=   7 distance = 7.812500e-03  = 0.0078125000 LI =        133 log2(LI) =   7 time =  0 seconds
n=   8 distance = 3.906250e-03  = 0.0039062500 LI =        261 log2(LI) =   8 time =  0 seconds
n=   9 distance = 1.953125e-03  = 0.0019531250 LI =        518 log2(LI) =   9 time =  0 seconds
n=  10 distance = 9.765625e-04  = 0.0009765625 LI =       1031 log2(LI) =  10 time =  0 seconds
n=  11 distance = 4.882812e-04  = 0.0004882812 LI =       2055 log2(LI) =  11 time =  0 seconds
n=  12 distance = 2.441406e-04  = 0.0002441406 LI =       4104 log2(LI) =  12 time =  0 seconds
n=  13 distance = 1.220703e-04  = 0.0001220703 LI =       8201 log2(LI) =  13 time =  0 seconds
n=  14 distance = 6.103516e-05  = 0.0000610352 LI =      16394 log2(LI) =  14 time =  0 seconds
n=  15 distance = 3.051758e-05  = 0.0000305176 LI =      32778 log2(LI) =  15 time =  0 seconds
n=  16 distance = 1.525879e-05  = 0.0000152588 LI =      65547 log2(LI) =  16 time =  0 seconds
n=  17 distance = 7.629395e-06  = 0.0000076294 LI =     131084 log2(LI) =  17 time =  0 seconds
n=  18 distance = 3.814697e-06  = 0.0000038147 LI =     262156 log2(LI) =  18 time =  0 seconds
n=  19 distance = 1.907349e-06  = 0.0000019073 LI =     524301 log2(LI) =  19 time =  0 seconds
n=  20 distance = 9.536743e-07  = 0.0000009537 LI =    1048590 log2(LI) =  20 time =  0 seconds
n=  21 distance = 4.768372e-07  = 0.0000004768 LI =    2097166 log2(LI) =  21 time =  0 seconds
n=  22 distance = 2.384186e-07  = 0.0000002384 LI =    4194319 log2(LI) =  22 time =  0 seconds
n=  23 distance = 1.192093e-07  = 0.0000001192 LI =    8388624 log2(LI) =  23 time =  1 seconds
n=  24 distance = 5.960464e-08  = 0.0000000596 LI =   16777233 log2(LI) =  24 time =  1 seconds
n=  25 distance = 2.980232e-08  = 0.0000000298 LI =   33554450 log2(LI) =  25 time =  2 seconds
n=  26 distance = 1.490116e-08  = 0.0000000149 LI =   67108884 log2(LI) =  26 time =  3 seconds
n=  27 distance = 7.450581e-09  = 0.0000000075 LI =  134217761 log2(LI) =  27 time =  8 seconds
n=  28 distance = 3.725290e-09  = 0.0000000037 LI =  268435874 log2(LI) =  28 time = 16 seconds
n=  29 distance = 1.862645e-09  = 0.0000000019 LI =  536883630 log2(LI) =  29 time = 30 seconds
n=  30 distance = 9.313226e-10  = 0.0000000009 LI = 1074147225 log2(LI) =  30 time = 62 seconds
....

So problem was cased be number precision.

[knighty]

I gess it's called a cancellation error. With double floats we have (1+2^(-53))-1=0 and (1+2^(-52))-1=2^(-52). So in the first case z==zf==1 and it will obviousely not escape for c=0.

With z=1+2^(-52) it will escape in 52 iterations. In general (and using arbitrary precision floats) with z=1+2^(-n), it needs n iteration to escape which is O(n).

in the second case with z=0.5+2^(-27) it needs approximately 2^27 iterations to escape. With z=0.5+2^(-n) it needs approximately 2^n iterations to escape which is O(2^n).

here is the (evaldraw) scripts I used:

Code:

(){
   cls(0);
   printf("c=0:\n");
   cx=0;cy=0;
   k=52;
   x=1;
   y=x+2^-k;
   z=(y-x)*2^k;
   printnum(z);
   printnum(iter(x,0,cx,cy,10^6));
   printnum(iter(y,0,cx,cy,10^6));
   printf("\nc=0.25:\n");
   cx=0.25;cy=0;
   k=19;
   x=0.5;
   y=x+2^-k;
   z=(y-x)*2^k;
   printnum(z);
   printnum(iter(x,0,cx,cy,10^6));
   printnum(iter(y,0,cx,cy,10^6));
}
iter(x,y,cx,cy,nmax){
   r2=x*x+y*y;
   for(i=0;i<nmax && r2<4;i++){
      xx=x*x-y*y+cx;
      y=2*x*y+cy;x=xx;
      r2=x*x+y*y;
   }
   i
}

[Adam Majewski]

Thx for answer.
Below is my C code using long double.
I thing about using gmp or gsl.
Help is wellcome 

Code:

// c code
#include <stdio.h>
#include <math.h> // pow
#include <float.h> // DBL_EPSILON
#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 distance; // distance(z1,zp)
  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;
  //distance = Zx-0.5; // distance(z1,zp)
 // printf("i= %.0f log2(i)= %.0f distance= %.38e = %.38f\n", i,log2(i), Zx, Zx );
 // printf("i= %.0Lf Zx = %.50Le\n", i,Zx);
  }
  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;
 
  //n=27;

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

 printf("DBL_EPSILON   = %e  = %40.38f ",DBL_EPSILON , DBL_EPSILON );
 

  return 0;
}

[Adam Majewski]

Is it possible to compute n for which  I should increase precision ?

[knighty]

Help is wellcome 

Regarding what exactly? 

Is it possible to compute n for which  I should increase precision ?

I'm not an expert but in this case it seems that that n is the number of bits in the mantissa +1.

[Adam Majewski]

Regarding what exactly? 
I'm not an expert but in this case it seems that that n is the number of bits in the mantissa +1.

Thx for answer.
Help is wellcome in ( for example) :
Converting C code to arbitrary precision ( for example MPRF like in http://jwm-art.net/mdz/ , or other libraries ).

Code:

/*
It does not work now !!!!!!!!!!!!!!!!


gcc  -lmpfr -lgmp -Wall -O2 pmpfr.c

mpfr without use of explicit definition of complex numbers


*/
#include <stdio.h>
// #include <math.h>
#include <gmp.h>
#include <mpfr.h>
//#include "<types.h>"
#include <time.h>

   // MPFR general settings
  mpfr_prec_t precision = 200; // 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),



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;
    mpfr_t zx2; // zx^2
    mpfr_t zy2; // zy^2

   
    mpfr_mul(   zy2,   zy,  zy,    GMP_RNDN); /* zy2 = zy * zy; */
    mpfr_mul(   zx2,   zx,  zx,    GMP_RNDN); /* zx2 = zx * zx; */
   
    while (mpfr_greater_p(ER2,temp))
    {
        /* zy = 2.0*zx*zy + cy; */
        mpfr_mul(   temp,     zx,    zy,    GMP_RNDN); // temp = zx*zy
        mpfr_mul_si(temp,     temp,     2.0,      GMP_RNDN); // temp = 2 * temp = 2*zx*zy
        mpfr_add(  zy,    temp,     cy,   GMP_RNDN);
        /* zx = zx^2 - zy^2 + cx; */
        mpfr_sub(   temp,     zx2,   zy2,   GMP_RNDN); // temp = zx^2 - zy^2
        mpfr_add(   zx,    temp,     cx,   GMP_RNDN); //  temp = temp + cx
        //
        mpfr_mul(   zy2,   zy,  zy,    GMP_RNDN); /* zy2 = zy * zy; */
        mpfr_mul(   zx2,   zx,  zx,    GMP_RNDN); /* zx2 = zx * zx; */
        i+=1.0;
    }
    return i;
}





int main (void)
{
   
 
  // declares variables

   
   double LastIteration;
   time_t start,end;
   double dif;
   unsigned int n;
   mpfr_t zpx, zpy;  // parabolic fixed point zp = zpx = zpy*I = 0.5 for c=1/4
   mpfr_t distance;
   mpfr_t base, exponent;
   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_t t;

  mpfr_set_default_prec(precision);
  // initializes MPFR variables
  mpfr_inits(ER, ER2,zpx, zpy, distance, base, exponent, t, (mpfr_ptr) 0 );

  // Assignment Functions
  mpfr_set_ld (zpx, 0.5, RoundMode);
  mpfr_set_ld (zpy, 0.0, RoundMode);
  mpfr_set_ld (base, 2.0, RoundMode);
  mpfr_set_ld (ER, 2.0, RoundMode);
  mpfr_set_ld (t, -1.0, RoundMode);
   
  // computations






 
 
   mpfr_mul(ER2,  ER,     ER,      GMP_RNDN); //

   

  for (n = 1; n <= 100; n++)
    { //mpfr_mul(n,n,  t, RoundMode); // !!!!!!!!!!!!!!!!!!!!!
     //mpfr_set_ld(exponent, ,  RoundMode);
   
     mpfr_pow_si(distance,base, -n, RoundMode);
     //LastIteration = GiveLastIteration_mpfr(zx, zy, cx, cy, ER2);
     mpfr_printf ("n = %d distance = %.17e \n", n,  distance);
     
    }

 

  // free the space used by the MPFR variables
  mpfr_clears (ER, ER2, zpx, zpy, distance,base, exponent,  (mpfr_ptr) 0);
  mpfr_free_cache ();  /* free the cache for constants like pi */
 

return 0;
}

[jwm-art]

In the function GiveLastIteration_mpfr, you are using the mpfr variables zx2, zy2, and temp, without initializing them beforehand.

[Adam Majewski]

here is beter version , but still not working. Help is wellcome 

Code:

/*

gcc  -lm -lmpfr -lgmp -Wall -O2 pmpfr.c

mpfr without use of explicit definition of complex numbers


*/
#include <stdio.h>
#include <math.h> // log2
#include <gmp.h>
#include <mpfr.h>
//#include "<types.h>"
#include <time.h>

   // MPFR general settings
  mpfr_prec_t precision = 200; // 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),



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;
    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,    GMP_RNDN); /* zy2 = zy * zy; */
    mpfr_mul(   zx2,   zx,  zx,    GMP_RNDN); /* zx2 = zx * zx; */
    mpfr_add(   temp,   zx2,  zy2,   GMP_RNDN); //  temp = zx2 + zy2
   
    while (mpfr_greater_p(ER2, temp)) //  while  (Zx2+Zy2<ER2
    {
        /* zy = 2.0*zx*zy + cy; */
        mpfr_mul(   temp,     zx,    zy,    GMP_RNDN); // temp = zx*zy
        mpfr_mul_ui(temp,     temp,     2.0,      GMP_RNDN); // temp = 2 * temp = 2*zx*zy
        mpfr_add(  zy,    temp,     cy,   GMP_RNDN);
        /* zx = zx^2 - zy^2 + cx; */
        mpfr_sub(   temp,     zx2,   zy2,   GMP_RNDN); // temp = zx^2 - zy^2
        mpfr_add(   zx,    temp,     cx,   GMP_RNDN); //  temp = temp + cx
        //
        mpfr_mul(   zy2,   zy,  zy,    GMP_RNDN); /* zy2 = zy * zy; */
        mpfr_mul(   zx2,   zx,  zx,    GMP_RNDN); /* zx2 = zx * zx; */
        mpfr_add(   temp,   zx2,  zy2,   GMP_RNDN); //  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, exponent;
   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_t t;

  mpfr_set_default_prec(precision);
  // initializes MPFR variables with default_prec
  mpfr_inits(ER, ER2, zpx, zpy, zx, zy, distance, cx, cy, base, exponent, t, (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);
  mpfr_set_ld (t, -1.0, RoundMode);
   
  // computations
  mpfr_mul(ER2,  ER,     ER,      GMP_RNDN); //
 
  for (n = 1; n <= 100; n++)
    {
     time (&start);
     mpfr_pow_si(distance,base, - n, RoundMode);
     mpfr_add(zx, zx, 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);
     //mpfr_printf ("n = %3d distance = %.10Re ; time = %2.0f seconds \n", n,  distance, dif);
    }

 

  // free the space used by the MPFR variables
  mpfr_clears (ER, ER2, zpx, zpy, zx, zy, distance, cx, cy, base, exponent,  (mpfr_ptr) 0);
  mpfr_free_cache ();  /* free the cache for constants like pi */
 

return 0;
}

[Adam Majewski]

line 109 should be :
mpfr_add(zx, zpx, distance, RoundMode);

Now it works

[Adam Majewski]

Here are results using MPFR-3.0.0-p8 with GMP-4.3.2 with precision = 80
It is about 20 times longer than gcc and long double.

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 =  1 seconds
n =  23 distance = 1.1920928955e-07 LI =    8388624 log2(LI) =  23; time =  3 seconds
n =  24 distance = 5.9604644775e-08 LI =   16777232 log2(LI) =  24; time =  4 seconds
n =  25 distance = 2.9802322388e-08 LI =   33554449 log2(LI) =  25; time = 10 seconds
n =  26 distance = 1.4901161194e-08 LI =   67108882 log2(LI) =  26; time = 20 seconds
n =  27 distance = 7.4505805969e-09 LI =  134217747 log2(LI) =  27; time = 39 seconds
n =  28 distance = 3.7252902985e-09 LI =  268435476 log2(LI) =  28; time = 78 seconds
n =  29 distance = 1.8626451492e-09 LI =  536870933 log2(LI) =  29; time = 152 seconds
n =  30 distance = 9.3132257462e-10 LI = 1073741846 log2(LI) =  30; time = 298 seconds
n =  31 distance = 4.6566128731e-10 LI = 2147483671 log2(LI) =  31; time = 596 seconds
n =  32 distance = 2.3283064365e-10 LI = 4294967320 log2(LI) =  32; time = 1206 seconds

To compile :
gcc  -lm -lmpfr -lgmp -Wall -O2 pmpfr.c

Code:

#include <stdio.h>
#include <math.h> // log2
#include <gmp.h>
#include <mpfr.h>
#include <time.h>

   // MPFR general settings
  mpfr_prec_t precision = 80; // 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),



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;
    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,    GMP_RNDN); /* zy2 = zy * zy; */
    mpfr_mul(   zx2,   zx,  zx,    GMP_RNDN); /* zx2 = zx * zx; */
    mpfr_add(   temp,   zx2,  zy2,   GMP_RNDN); //  temp = zx2 + zy2
   
    while (mpfr_greater_p(ER2, temp)) //  while  (Zx2+Zy2<ER2
    {
        /* zy = 2.0*zx*zy + cy; */
        mpfr_mul(   temp,     zx,    zy,    GMP_RNDN); // temp = zx*zy
        mpfr_mul_ui(temp,     temp,     2.0,      GMP_RNDN); // temp = 2 * temp = 2*zx*zy
        mpfr_add(  zy,    temp,     cy,   GMP_RNDN);
        /* zx = zx^2 - zy^2 + cx; */
        mpfr_sub(   temp,     zx2,   zy2,   GMP_RNDN); // temp = zx^2 - zy^2
        mpfr_add(   zx,    temp,     cx,   GMP_RNDN); //  temp = temp + cx
        //
        mpfr_mul(   zy2,   zy,  zy,    GMP_RNDN); /* zy2 = zy * zy; */
        mpfr_mul(   zx2,   zx,  zx,    GMP_RNDN); /* zx2 = zx * zx; */
        mpfr_add(   temp,   zx2,  zy2,   GMP_RNDN); //  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, exponent;
   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_t t;

  mpfr_set_default_prec(precision);
  // initializes MPFR variables with default_prec
  mpfr_inits(ER, ER2, zpx, zpy, zx, zy, distance, cx, cy, base, exponent, t, (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);
  mpfr_set_ld (t, -1.0, RoundMode);
   
  // computations
  mpfr_mul(ER2,  ER,     ER,      GMP_RNDN); //
 
  for (n = 1; n <= 100; 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);
    }
  printf ("Using MPFR-%s with GMP-%s with precision = %u bits \n", mpfr_version, gmp_version, (unsigned int) precision);
 

  // free the space used by the MPFR variables
  mpfr_clears (ER, ER2, zpx, zpy, zx, zy, distance, cx, cy, base, exponent,  (mpfr_ptr) 0);
  mpfr_free_cache ();  /* free the cache for constants like pi */
 

return 0;
}

[Adam Majewski]

It seems that precision ( the number of bits used to represent the significand of a floating-point number )
is :

precision > 2* nMax

[Adam Majewski]

See :

http://en.wikibooks.org/wiki/Fractals/Mathematics/Numerical
for results