News: Check out the originating "3d Mandelbulb" thread here  ## The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!

 Pages:    Go Down       Author Topic: numerical problem with bailout test  (Read 2137 times) Description: strange result of numerical computations 0 Members and 1 Guest are viewing this topic.
Fractal Lover  Posts: 221  « 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
Fractal Lover  Posts: 221  « 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
Fractal Lover  Posts: 221  « Reply #3 on: April 21, 2012, 09:43:16 PM »

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;
} Logged
Fractal Lover  Posts: 221  « 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 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. Logged
Fractal Lover  Posts: 221  « Reply #6 on: April 23, 2012, 09:18:16 PM »

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.

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
/* 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   « 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
Fractal Lover  Posts: 221  « Reply #8 on: April 24, 2012, 07:08:24 PM »

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
/* 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);
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
Fractal Lover  Posts: 221  « Reply #9 on: April 24, 2012, 07:14:43 PM »

line 109 should be :

Now it works Logged
Fractal Lover  Posts: 221  « 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
/* 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);
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
Fractal Lover  Posts: 221  « 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
Fractal Lover  Posts: 221  « Reply #12 on: May 12, 2012, 05:15:50 PM »

See :

http://en.wikibooks.org/wiki/Fractals/Mathematics/Numerical
for results Logged
 Pages:    Go Down
 Related Topics Subject Started by Replies Views Last post  Numerical Dream Images Showcase (Rate My Fractal) lyron 0 746 April 04, 2011, 10:11:32 AM by lyron  [link to] Stand-alone code for numerical computing Programming ker2x 0 1110 June 10, 2011, 02:31:39 PM by ker2x  numerical monsters Images Showcase (Rate My Fractal) Eric B 0 403 February 07, 2013, 05:46:12 AM by Eric B  Numerical precision of formula values feature request cyseal 0 450 December 04, 2015, 09:14:49 AM by cyseal  Test re: Tiling problem Fragmentarium Gallery Sabine 0 242 September 04, 2016, 09:56:23 AM by Sabine