Logo by stereoman - Contribute your own Logo!

END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

it was a great time but no longer maintainable by c.Kleinhuis contact him for any data retrieval,
thanks and see you perhaps in 10 years again

this forum will stay online for reference
News: Support us via Flattr FLATTR Link
 
*
Welcome, Guest. Please login or register. January 16, 2018, 08:11:21 AM


Login with username, password and session length


The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!


Pages: [1]   Go Down
  Print  
Share this topic on DiggShare this topic on FacebookShare this topic on GoogleShare this topic on RedditShare this topic on StumbleUponShare this topic on Twitter
Author Topic: numerical problem with bailout test  (Read 1779 times)
Description: strange result of numerical computations
0 Members and 1 Guest are viewing this topic.
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« on: April 20, 2012, 08:23:08 PM »

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



« Last Edit: April 20, 2012, 08:58:41 PM by Adam Majewski, Reason: spelling » Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #1 on: April 21, 2012, 01:54:59 PM »

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.

Logged
knighty
Fractal Iambus
***
Posts: 819


« Reply #2 on: April 21, 2012, 04:59:12 PM »

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
}
Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #3 on: April 21, 2012, 09:43:16 PM »

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

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

Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #4 on: April 22, 2012, 09:14:38 AM »

Is it possible to compute n for which  I should increase precision ?
« Last Edit: April 22, 2012, 09:56:25 AM by Adam Majewski » Logged
knighty
Fractal Iambus
***
Posts: 819


« Reply #5 on: April 22, 2012, 10:54:58 AM »

Help is wellcome  smiley
Regarding what exactly?  smiley

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.
Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #6 on: April 23, 2012, 09:18:16 PM »

Regarding what exactly?  smiley
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;
}



Logged
jwm-art
Iterator
*
Posts: 171



WWW
« Reply #7 on: April 24, 2012, 12:25:40 AM »

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

Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #8 on: April 24, 2012, 07:08:24 PM »

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

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


Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #9 on: April 24, 2012, 07:14:43 PM »

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

Now it works
Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #10 on: April 24, 2012, 09:04:04 PM »

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

Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #11 on: April 25, 2012, 06:21:24 PM »

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

precision > 2* nMax
Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #12 on: May 12, 2012, 05:15:50 PM »

See :

http://en.wikibooks.org/wiki/Fractals/Mathematics/Numerical
for results
Logged
Pages: [1]   Go Down
  Print  
 
Jump to:  

Related Topics
Subject Started by Replies Views Last post
Numerical Dream Images Showcase (Rate My Fractal) lyron 0 635 Last post April 04, 2011, 10:11:32 AM
by lyron
[link to] Stand-alone code for numerical computing Programming ker2x 0 881 Last post June 10, 2011, 02:31:39 PM
by ker2x
numerical monsters Images Showcase (Rate My Fractal) Eric B 0 305 Last post February 07, 2013, 05:46:12 AM
by Eric B
Numerical precision of formula values feature request cyseal 0 267 Last post December 04, 2015, 09:14:49 AM
by cyseal
Test re: Tiling problem Fragmentarium Gallery Sabine 0 79 Last post September 04, 2016, 09:56:23 AM
by Sabine

Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines

Valid XHTML 1.0! Valid CSS! Dilber MC Theme by HarzeM
Page created in 0.512 seconds with 27 queries. (Pretty URLs adds 0.019s, 2q)