DBL_EPSILON i Calculated Machine epsilon

Showing 1-10 of 10 messages
DBL_EPSILON i Calculated Machine epsilon Adam Majewski 4/15/12 6:13 AM
Witam.

Napisałem program w c ( kod poniżej). Oblicza ile potrzeba kroków (
iteracji tu LastIteration) aby punkt z1 uciekł z koła o promieniu ER,
czyli sprawdza abs(z1)<ER




Program działa, ale dla n=27 wpada w nieskończoną pętlę.
Przypuszczam, że wtedy program nie rozróżnia z1 i zp :
bo wtedy distance(z1,zp)<Calculated  Machine epsilon



gdzie :
distance(z1,zp) = 1.490116e-08
Calculated  Machine epsilon = 1.19209E-07 ( używając programu z wikipedii )
http://en.wikipedia.org/wiki/Machine_epsilon
punkt jest liczbą zespoloną)
punkt zp jest punktem stałym = nie ucieka z koła o promieniu ER


Jak to się ma do :

DBL_EPSILON           = 2.220446e-16  =
0.00000000000000022204460492503130808473 ( z float.h )
?



Pozdrawiam

Adam












============== Uruchomienie==========================

gcc -lm -Wall i.c
./a.out


=============wyniki ====================


n=  1 distance(z1,zp) = 5.000000e-01  LastIteration =        3
n=  2 distance(z1,zp) = 2.500000e-01  LastIteration =        5
n=  3 distance(z1,zp) = 1.250000e-01  LastIteration =       10
n=  4 distance(z1,zp) = 6.250000e-02  LastIteration =       19
n=  5 distance(z1,zp) = 3.125000e-02  LastIteration =       35
n=  6 distance(z1,zp) = 1.562500e-02  LastIteration =       68
n=  7 distance(z1,zp) = 7.812500e-03  LastIteration =      133
n=  8 distance(z1,zp) = 3.906250e-03  LastIteration =      261
n=  9 distance(z1,zp) = 1.953125e-03  LastIteration =      518
n= 10 distance(z1,zp) = 9.765625e-04  LastIteration =     1031
n= 11 distance(z1,zp) = 4.882812e-04  LastIteration =     2055
n= 12 distance(z1,zp) = 2.441406e-04  LastIteration =     4104
n= 13 distance(z1,zp) = 1.220703e-04  LastIteration =     8201
n= 14 distance(z1,zp) = 6.103516e-05  LastIteration =    16394
n= 15 distance(z1,zp) = 3.051758e-05  LastIteration =    32778
n= 16 distance(z1,zp) = 1.525879e-05  LastIteration =    65547
n= 17 distance(z1,zp) = 7.629395e-06  LastIteration =   131084
n= 18 distance(z1,zp) = 3.814697e-06  LastIteration =   262156
n= 19 distance(z1,zp) = 1.907349e-06  LastIteration =   524301
n= 20 distance(z1,zp) = 9.536743e-07  LastIteration =  1048590
n= 21 distance(z1,zp) = 4.768372e-07  LastIteration =  2097167
n= 22 distance(z1,zp) = 2.384186e-07  LastIteration =  4194322
n= 23 distance(z1,zp) = 1.192093e-07  LastIteration =  8388675
n= 24 distance(z1,zp) = 5.960464e-08  LastIteration = 16778821
n= 25 distance(z1,zp) = 2.980232e-08  LastIteration = 33604781
n= 26 distance(z1,zp) = 1.490116e-08  LastIteration = 68563774


=================== kod ==================


#include <stdio.h>
#include <math.h> // pow
#include <float.h> // DBL_EPSILON



double GiveLastIteration(double Zx, double Zy, double Cx, double Cy, int
IterationMax, 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;
    //printf("i = %f Zx = %e = %40.38f \n", i,Zx, Zx);
   }
   return i;
  }


int main()

{
   double distance;

    int n;
    // parabolic fixed point zp = zpx = zpy*I = 0.5
    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;


    //n = 1000;
    ER2= ER*ER;



    for (n =1; n<101; n++)
    {
      distance = pow(2.0,-n);
      printf("n= %d distance(z1,zp) = %e  = %40.38f ",n, distance,
distance);
      z1x = zpx + distance;
      LastIteration = GiveLastIteration(z1x,z1y, cx,cy, n+n,ER2 );
      printf("LastIteration = %.0f \n", LastIteration);
    }

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


   return 0;
}
Re: DBL_EPSILON i Calculated Machine epsilon Adam Majewski 4/16/12 9:26 AM
Zapomniałem dodać link do pracy skąd zaczerpnąłem pomysł na taki program  :
http://arxiv.org/abs/math/0505036


hth


Adam
Re: DBL_EPSILON i Calculated Machine epsilon Adam Majewski 4/22/12 12:55 AM

http://www.fractalforums.com/programming/numerical-problem-with-bailout-test/
Re: DBL_EPSILON i Calculated Machine epsilon Maciej Witkowiak 4/28/12 5:39 AM
Adam Majewski napisał:

> Witam.
>
> Napisałem program w c ( kod poniżej). Oblicza ile potrzeba kroków (
> iteracji tu LastIteration) aby punkt z1 uciekł z koła o promieniu ER,
> czyli sprawdza abs(z1)<ER

(...)

> Jak to się ma do :
>
> DBL_EPSILON           = 2.220446e-16  =
> 0.00000000000000022204460492503130808473 ( z float.h )
> ?

Przy wywołaniu:
gcc -lm ...
oraz
gcc -lm -mfpmath=sse ...
dostałem identyczny wynik, jak Ty.

Ale po kompilacji przez:
gcc -lm -fpmath=387 ...
pojawiły się nowe wiersze wyniku:

n= 27 distance(z1,zp) = 7.450581e-09  = 0.00000000745058059692382812500000000000 LastIteration = 134217761
n= 28 distance(z1,zp) = 3.725290e-09  = 0.00000000372529029846191406250000000000 LastIteration = 268435874
n= 29 distance(z1,zp) = 1.862645e-09  = 0.00000000186264514923095703125000000000 LastIteration = 536883630
n= 30 distance(z1,zp) = 9.313226e-10  = 0.00000000093132257461547851562500000000 LastIteration = 1074147225
n= 31 distance(z1,zp) = 4.656613e-10  = 0.00000000046566128730773925781250000000 LastIteration = 2160053374
n= 32 distance(z1,zp) = 2.328306e-10  = 0.00000000023283064365386962890625000000 LastIteration = 4245834220

Tak na pierwszy rzut oka - może lepiej przepisać ten program tak, żeby używał
GMP?

I.

--
Najlepsza sygnatura to brak sygnatury.
http://bossstation.dnsalias.org/
Re: DBL_EPSILON i Calculated Machine epsilon Adam Majewski 5/13/12 4:06 AM
Dzięki za odpowiedź. Nie testowałem ( czytaj nie znałem ) tych opcji.

http://gcc.gnu.org/onlinedocs/gcc-4.0.0/gcc/i386-and-x86_002d64-Options.html
Ja mam AMD Athlon 64 i ubuntu 64 bitowe więc prawdopodobnie kompilator wybiera
je automatycznie.
Dziwne (?) że dla 387 daje więcej wyników.(?)

Co do GMP to zobacz tu :
http://en.wikibooks.org/wiki/Fractals/Mathematics/Numerical
Jest tam program używający MPFR ( i GMP) pozwala na większą pracyzję obliczeń ale ceną jest czas. Wg moich szacunków testowanie 1 punktu dla n = 50 zajmie rok. Popraw mnie jeśli się pomyliłem.
Wydaje mi się że wiem skąd jest ten problem :

precyzja każdego typu jest ograniczona do ( w uproszczeniu) m liczby cyfr dziesiętnych. Jeśli dodam małą liczbę distance do dużej zf  to jeszcze mieszczę się w zakresie przechowywanych cyfr dziesiętnych. Ale w trakcie działań następuje potęgowanie z2 = z * z
i wtedy cześć związana z distance przesuwa się w prawo i wychodzi poza zakres przechowywanych liczb. Jest obcięta czyli z= zf . Powoduje to powstanie nieskończonej pętli.
Podwyżąszenie precyzji odsuwa problem ale go nie likwiduje.

Pomoc mile widziana.

Pozdrawiam
Re: DBL_EPSILON i Calculated Machine epsilon Adam Majewski 5/13/12 4:23 AM
On Sunday, April 15, 2012 3:13:51 PM UTC+2, Adam Majewski wrote:
POzostaje czekać na GMP dla GPU (:-)
Re: DBL_EPSILON i Calculated Machine epsilon Edek Pienkowski 5/15/12 11:10 AM
Dnia Sun, 13 May 2012 04:23:41 -0700, Adam Majewski napisal:

>
> POzostaje czekać na GMP dla GPU (

Czyli naprawdę jeszcze nie ma?

Edek
Re: DBL_EPSILON i Calculated Machine epsilon Adam Majewski 5/15/12 2:23 PM
Popraw mnie jeśli się mylę : oficjalna strona nie ma
Myślisz o :
http://code.google.com/p/gpuprec/
?
??
Re: DBL_EPSILON i Calculated Machine epsilon Adam Majewski 5/15/12 2:27 PM
lub
http://www.hpcs.cs.tsukuba.ac.jp/~nakayama/cump/
Re: DBL_EPSILON i Calculated Machine epsilon Edek Pienkowski 5/15/12 3:38 PM
Dnia Tue, 15 May 2012 14:27:56 -0700, Adam Majewski napisal:
No, liczyłem na to że coś znajdziesz, jutro sam będę miał gotowe ;p

Edek