h1. Ciurul lui Erathostene
(Creat de '_Cosmin_':user/Cosmin la data de _2004-11-24_ categoria _Teoria numerelor_, autor(i) _Cosmin_)
*Continut scurt:*
Articolul de fata incearca o implementarea mai eficienta a acestui algoritm clasic. Se poate optimiza pentru a folosi doar O(sqrt(n)) memorie, varianta prezentata aici folosind O(n / log n) memorie, unde log n e numarul de biti al unui cuvant.
*Continut lung:*
Ciurul lui Erathostene e un algoritm clasic care se invata la scoala impreuna cu conceptul de numere prime inca din clasa a 6. Acest algoritm determina toate numerele prime mai mici decat un numar dat ca parametru. Ideea lui de abordare a acestei probleme poate fi modificata pentru a rezolva si alte probleme precum problema Fractii din arhiva infoarena, problema Riemann vs Mertens din arhiva uva ([1]http://acm.uva.es/p/v107/10738.html), problema Divizibilitate a rundei 4 a concursului algoritmus ([2]http://algoritmus.org/probleme/Probleme_Runda04.php), sau problema Square Free a SRM-ului 190 de pe TOPCODER ([3]http://www.topcoder.com/stat?c=problem_statement&pm=2342&rd=4770, rezolvarea ei o gasiti aici : [4]http://www.topcoder.com/index?t=statistics&c=srm190_prob) . (linkurile ce contin litere mari nu merg pentru ca editorul le trasforma automat in litere mici si nu mai pot fi modificate, daca vreti sa urmati un link faceti un copy paste cu adresa in browserul vostru)
In articolul acesta nu ne vom concentra asupra acestor probleme ci asupra implementarii optimizate ale algoritmului original.
Am mai scris despre aceste optimizari intr-un articol din ginfo, dar codul de acolo nu era testat si nu merge :).
Ideea la acest algoritm e ca marcam intr-un sir fiecare multiplu al unui numar prim si numerele ramase nemarcate sunt numere prime, o descriere mai grafica gasiti la adresa: [5]http://mathworld.wolfram.com/SieveofEratosthenes.html.
Sa incercam o prima implementare a acestui algoritm (implementarile vor folosi limbajul java, dar sunt foarte usor transformabile in C/C++).
//class PrimeNumbersSieve1
final int MAXSIZE = 1000001;
char[] p = new char[MAXSIZE];
//p[i] == 0 if i is prime
public int getTheNumber(int n) {
int i, j, nr = 0;
for (i = 2; i <= n; ++i) {
if (p[i] == 0) {
nr++;
for (j = i + i; j <= n; j += i) {
p[j] = 1;
}
}
}
return nr;
}
O prima idee de optimizare ar fi sa nu mai luam in calcul numerele pare pentru ca stim ca singurul numar prim par e 2. Deci sa vedem noua varianta a programului:
//class PrimeNumbersSieve2
final int MAXSIZE = 1000001;
char[] p = new char[MAXSIZE];
//p[i] == 0 if i is prime
public int getTheNumber(int n) {
int i, j, nr = 1;
for (i = 3; i <= n; i += 2) {
if (p[i] == 0) {
nr++;
for (j = i + i + i; j <= n; j += i << 1) {
p[j] = 1;
}
}
}
return nr;
}
Putem incerca o optimizare de memorie, pentru ca nu mai avem nevoie de elementele cu index par din sirul p. Acum semnificatia lui p[i] s-a schimbat p[i] fiind 0 daca 2*i+1 e numar prim si 1 daca nu.
// class PrimeNumbersSieve3
final int MAXSIZE = 1000000/2+1;
char[] p = new char[MAXSIZE];
//p[i] == 0 if 2*i + 1 is prime
public int getTheNumber(int n) {
int i, j, nr = 1;
for (i = 1; (i << 1) + 1 <= n; i += 1) {
if (p[i] == 0) {
nr++;
for (j = i + i + i + 1;
(j << 1) + 1 <= n;
j += (i << 1) + 1) {
p[j] = 1;
}
}
}
return nr;
}
Urmatoarea optimizare va fi marcarea multiplilor numarului prim i de la i*i nu de la 2*i cum am facut in prima varianta sau de la 3*i cum am facut in a 2-a. Aceasta optimizare este evidenta: orice numar prim compus multiplu de i mai mic decat i*i are un factor prim mai mic decat i, si acel factor l-a marcat mai devreme, deci nu are rost sa il marcam si la pasul i. Ideea aceasta este exact ideea ce se foloseste la rezolvarea problemei Numere Prime din arhiva infoarena. Sa vedem acum codul sursa:
//class PrimeNumbersSieve4
final int MAXSIZE = 1000000/2+1;
char[] p = new char[MAXSIZE];
//p[i] == 0 if 2*i + 1 is prime
public int getTheNumber(int n) {
int i, j, nr = 1;
for (i = 1;
((i * i) << 1) + (i << 1) <= n;
i += 1) {
if (p[i] == 0) {
for (j = ((i * i) << 1) + (i << 1);
(j << 1) + 1 <= n;
j += (i << 1) + 1) {
p[j] = 1;
}
}
}
for (i=1; 2 * i + 1 <= n; ++i)
if (p[i] == 0) nr++;
return nr;
}
Codul sursa arata putin urat pentru ca nu lucram direct cu i ci cu 2*i+1, am mai facut optimizarea ce apare si in mathworld, nu parcurgem numerele pana la n pentru marcarea multiplilor ci pana la sqrt(n) lucru care e evident dupa cele explicate mai sus.
Ultima imbunatatire care o vom aduce este aceea de a folosi mai putina memorie. Cum pentru fiecare numar e necesara doar o informatie booleana, aceasta o putem tine intr-un bit, nu este necesar un char intreg. Sa vedem cum arata codul:
//class PrimeNumbersSieve5
final int MAXSIZE = 100000000/2/8+1;
char[] p = new char[MAXSIZE];
//p[i] == 0 if 2*i + 1 is prime
public int getTheNumber(int n) {
int i, j, nr = 1;
for (i = 1;
((i * i) << 1) + (i << 1) <= n;
i += 1) {
if ((p[i >> 3] & (1 << (i & 7))) == 0) {
for (j = ((i * i) << 1) + (i << 1);
(j << 1) + 1 <= n;
j += (i << 1) + 1) {
p[j >> 3] |= (1 << (j & 7));
}
}
}
for (i = 1; 2 * i + 1 <= n; ++i)
if ((p[i >> 3] & (1 << (i & 7))) == 0)
nr++;
return nr;
}
Codul p[i >> 3] & (1 << (i & 7)) == 0 testeaza daca al i-lea bit din sirul de biti e 0 (deci daca 2*i+1 e prim). i >> 3 e echivalent cu i / 8 deci se gaseste pt bitul al i-lea in ce char e el stocat din sirul nostru de charuri. i & 7 e echivalent cu i % 8 ca sa aflam pe ce bit al charului e stocat numarul prim i.
Codul p[j >> 3] |= (1 << (j & 7)) seteaza bitul j % 8 al charului j / 8 in 1, pt ca sa stim ca numarul 2 * j + 1 nu e prim.
Ultima varianta arata destul de urat fata de prima, sa vedem daca s-a meritat efortul. Urmatoarele date sunt obtinute pe un procesor Athlon XP la 1.8 Ghz:
10 runs of BruteForcePrimes().method1(1000000) lasted: 24.956 seconds
10 runs of PrimeNumbersSieve1().getTheNumber(1000000) lasted: 3.024 seconds
10 runs of PrimeNumbersSieve2().getTheNumber(1000000) lasted: 1.842 seconds
10 runs of PrimeNumbersSieve3().getTheNumber(1000000) lasted: 1.422 seconds
10 runs of PrimeNumbersSieve4().getTheNumber(1000000) lasted: 0.971 seconds
10 runs of PrimeNumbersSieve5().getTheNumber(1000000) lasted: 0.22 seconds
Codul folosit pentru testare e:
public static void main(String[] arg) {
double tick = System.currentTimeMillis();
for (int i = 0; i < 10; ++i) {
System.out.println("Run " + i
+ " result: "+
new PrimeNumbersSieve1().
getTheNumber(1000000));
}
System.out.println(" 10 runs of "
+ "PrimeNumbersSieve1()."
+ "getTheNumber(1000000) "
+ "lasted: "
+ (System.currentTimeMillis()
- tick)
* 1e-3 + " seconds");
}
S-a rulat fiecare cod de 10 ori pentru a da ocazia optimizarilor Just In Time din java sa isi faca treaba.
Se pare ca optimizarile codului initial au dus la o imbunatatire a vitezei cu un factor de 10.
Problemele sugerate la inceputul articolului sunt destul de frumoase si ar trebui rezolvate. Pt fiecare dintre ele puteti sa va testati rezolvarea: la algoritmus sunt disponibile fisierele de intrare si iesire, infoarena si uva au evaluatoare automate si pentru problema de pe TOPCODER puteti intra in applet-ul lor si sa incercati sa o rezolvati in practice roomul asociat SRM-ului 190 dupa ce o rezolvati folositi optiunea Run System Test din Practice Room Options.
References
Visible links
1. http://acm.uva.es/p/v107/10738.html
2. http://algoritmus.org/probleme/probleme_runda04.php
3. http://www.topcoder.com/stat?c=problem_statement&pm=2342&rd=4770
4. http://www.topcoder.com/index?t=statistics&c=srm190_prob
5. http://mathworld.wolfram.com/sieveoferatosthenes.html
h1. Ciurul lui Eratostene
(Categoria _Matematica_, Autor _Cosmin Negruseri_)
Articolul de fata incearca o implementarea mai eficienta a acestui algoritm clasic. Se poate optimiza pentru a folosi doar $O(sqrt(n))$ memorie, varianta prezentata aici folosind $O(n / log n)$ memorie, unde log n e numarul de biti al unui cuvant.
Ciurul lui Eratostene e un algoritm clasic care se invata la scoala impreuna cu conceptul de numere prime inca din clasa a 6. Acest algoritm determina toate numerele prime mai mici decat un numar dat ca parametru. Ideea lui de abordare a acestei probleme poate fi modificata pentru a rezolva si alte probleme precum problema "Fractii":problema/Fractii din arhiva infoarena, problema "Riemann vs Mertens":http://acm.uva.es/p/v107/10738.html din arhiva "uva":http://acm.uva.es/, problema "Divizibilitate":http://algoritmus.org/probleme/Probleme_Runda04.php a rundei 4 a concursului algoritmus, sau problema "Square Free":http://www.topcoder.com/stat?c=problem_statement&pm=2342&rd=4770 a SRM-ului 190 de pe TOPCODER, rezolvarea ei o gasiti "aici":http://www.topcoder.com/index?t=statistics&c=srm190_prob.
In articolul acesta nu ne vom concentra asupra acestor probleme ci asupra implementarii optimizate ale algoritmului original.
Am mai scris despre aceste optimizari intr-un articol din ginfo, dar codul de acolo nu era testat si nu merge :).
Ideea la acest algoritm e ca marcam intr-un sir fiecare multiplu al unui numar prim si numerele ramase nemarcate sunt numere prime, o descriere mai grafica gasiti la "adresa":http://mathworld.wolfram.com/SieveofEratosthenes.html.
Sa incercam o prima implementare a acestui algoritm (implementarile vor folosi limbajul Java, dar sunt foarte usor transformabile in C/C++).
== code(java) |
//class PrimeNumbersSieve1
final int MAXSIZE = 1000001;
char[] p = new char[MAXSIZE];
//p[i] == 0 if i is prime
public int getTheNumber(int n) {
int i, j, nr = 0;
for (i = 2; i <= n; ++i) {
if (p[i] == 0) {
nr++;
for (j = i + i; j <= n; j += i) {
p[j] = 1;
}
}
}
return nr;
}
==
O prima idee de optimizare ar fi sa nu mai luam in calcul numerele pare pentru ca stim ca singurul numar prim par e {$2$}. Deci sa vedem noua varianta a programului:
== code(java) |
//class PrimeNumbersSieve2
final int MAXSIZE = 1000001;
char[] p = new char[MAXSIZE];
//p[i] == 0 if i is prime
public int getTheNumber(int n) {
int i, j, nr = 1;
for (i = 3; i <= n; i += 2) {
if (p[i] == 0) {
nr++;
for (j = i + i + i; j <= n; j += i << 1) {
p[j] = 1;
}
}
}
return nr;
}
==
Putem incerca o optimizare de memorie, pentru ca nu mai avem nevoie de elementele cu index par din sirul {$p$}. Acum semnificatia lui $p{~i~}$ s-a schimbat $p{~i~}$ fiind 0 daca $2*i+1$ e numar prim si $1$ daca nu.
== code(java) |
// class PrimeNumbersSieve3
final int MAXSIZE = 1000000/2+1;
char[] p = new char[MAXSIZE];
//p[i] == 0 if 2*i + 1 is prime
public int getTheNumber(int n) {
int i, j, nr = 1;
for (i = 1; (i << 1) + 1 <= n; i += 1) {
if (p[i] == 0) {
nr++;
for (j = i + i + i + 1; (j << 1) + 1 <= n; j += (i << 1) + 1) {
p[j] = 1;
}
}
}
return nr;
}
==
Urmatoarea optimizare va fi marcarea multiplilor numarului prim $i$ de la $i*i$ nu de la $2*i$ cum am facut in prima varianta sau de la $3*i$ cum am facut in a {$2$}-a. Aceasta optimizare este evidenta: orice numar prim compus multiplu de $i$ mai mic decat $i*i$ are un factor prim mai mic decat {$i$}, si acel factor l-a marcat mai devreme, deci nu are rost sa il marcam si la pasul {$i$}. Ideea aceasta este exact ideea ce se foloseste la rezolvarea problemei "Numere Prime":problema/prim din arhiva infoarena. Sa vedem acum codul sursa:
== code(java) |
//class PrimeNumbersSieve4
final int MAXSIZE = 1000000/2+1;
char[] p = new char[MAXSIZE];
//p[i] == 0 if 2*i + 1 is prime
public int getTheNumber(int n) {
int i, j, nr = 1;
for (i = 1; ((i * i) << 1) + (i << 1) <= n; i += 1) {
if (p[i] == 0) {
for (j = ((i * i) << 1) + (i << 1); (j << 1) + 1 <= n; j += (i << 1) + 1) {
p[j] = 1;
}
}
}
for (i=1; 2 * i + 1 <= n; ++i)
if (p[i] == 0) nr++;
return nr;
}
==
Codul sursa arata putin urat pentru ca nu lucram direct cu $i$ ci cu {$2*i+1$}, am mai facut optimizarea ce apare si in "mathworld":http://mathworld.wolfram.com/, nu parcurgem numerele pana la $n$ pentru marcarea multiplilor ci pana la $sqrt(n)$ lucru care e evident dupa cele explicate mai sus.
Ultima imbunatatire care o vom aduce este aceea de a folosi mai putina memorie. Cum pentru fiecare numar e necesara doar o informatie booleana, aceasta o putem tine intr-un bit, nu este necesar un char intreg. Sa vedem cum arata codul:
== code(java) |
//class PrimeNumbersSieve5
final int MAXSIZE = 100000000/2/8+1;
char[] p = new char[MAXSIZE];
//(p[i>>3]&(1<<(i&7))) == 0 if 2*i + 1 is prime
public int getTheNumber(int n) {
int i, j, nr = 1;
for (i = 1; ((i * i) << 1) + (i << 1) <= n; i += 1) {
if ((p[i >> 3] & (1 << (i & 7))) == 0) {
for (j = ((i * i) << 1) + (i << 1); (j << 1) + 1 <= n; j += (i << 1) + 1) {
p[j >> 3] |= (1 << (j & 7));
}
}
}
for (i = 1; 2 * i + 1 <= n; ++i)
if ((p[i >> 3] & (1 << (i & 7))) == 0)
nr++;
return nr;
}
==
Codul {@p[i >> 3] & (1 << (i & 7)) @}=={@ 0@} testeaza daca al {@i@}-lea bit din sirul de biti e $0$ (deci daca {@2*i+1@} e prim). {@i >> 3@} e echivalent cu {@i / 8@} deci se gaseste pt bitul al {@i@}-lea in ce char e el stocat din sirul nostru de charuri. {@i & 7@} e echivalent cu {@i % 8@} ca sa aflam pe ce bit al charului e stocat numarul prim {@i@}.
Codul {@p[j >> 3] |= (1 << (j & 7))@} seteaza bitul {@j % 8@} al charului {@j / 8@} in {@1@}, pentru ca sa stim ca numarul {@2 * j + 1@} nu e prim.
Ultima varianta arata destul de urat fata de prima, sa vedem daca s-a meritat efortul. Urmatoarele date sunt obtinute pe un procesor Athlon XP la 1.8 Ghz:
* 10 runs of {@BruteForcePrimes().method1(1000000)@} lasted: 24.956 seconds
* 10 runs of {@PrimeNumbersSieve1().getTheNumber(1000000)@} lasted: 3.024 seconds
* 10 runs of {@PrimeNumbersSieve2().getTheNumber(1000000)@} lasted: 1.842 seconds
* 10 runs of {@PrimeNumbersSieve3().getTheNumber(1000000)@} lasted: 1.422 seconds
* 10 runs of {@PrimeNumbersSieve4().getTheNumber(1000000)@} lasted: 0.971 seconds
* 10 runs of {@PrimeNumbersSieve5().getTheNumber(1000000)@} lasted: 0.22 seconds
Codul folosit pentru testare e:
== code(java) |
public static void main(String[] arg) {
double tick = System.currentTimeMillis();
for (int i = 0; i < 10; ++i) {
System.out.println("Run " + i
+ " result: "+
new PrimeNumbersSieve1().
getTheNumber(1000000));
}
System.out.println(" 10 runs of "
+ "PrimeNumbersSieve1()."
+ "getTheNumber(1000000) "
+ "lasted: "
+ (System.currentTimeMillis()
- tick)
* 1e-3 + " seconds");
}
==
S-a rulat fiecare cod de $10$ ori pentru a da ocazia optimizarilor Just In Time din java sa isi faca treaba.
Se pare ca optimizarile codului initial au dus la o imbunatatire a vitezei cu un factor de {$10$}.
Problemele sugerate la inceputul articolului sunt destul de frumoase si ar trebui rezolvate. Pentru fiecare dintre ele puteti sa va testati rezolvarea: la algoritmus sunt disponibile fisierele de intrare si iesire, infoarena si uva au evaluatoare automate si pentru problema de pe "TOPCODER":http://www.topcoder.com/tc puteti intra in applet-ul lor si sa incercati sa o rezolvati in practice roomul asociat SRM-ului 190 dupa ce o rezolvati folositi optiunea Run System Test din Practice Room Options.
