Adam Majewski
|
|
« 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/0505036TIA
|
|
« Last Edit: April 20, 2012, 08:58:41 PM by Adam Majewski, Reason: spelling »
|
Logged
|
|
|
|
Adam Majewski
|
|
« 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: (){ 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
|
|
« 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 // 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
|
|
« 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
|
|
|
|
Adam Majewski
|
|
« 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. 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 ). /* 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
|
|
« 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
|
|
« Reply #8 on: April 24, 2012, 07:08:24 PM » |
|
here is beter version , but still not working. Help is wellcome /*
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
|
|
« 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
|
|
« 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 #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
|
|
« 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
|
|
|
|
|
|