V
Firstrun schrieb:
Ich hab das Problem gelösst für ein Feld.Da beim Feld die Zahl der Primazahlen a priori fest stehen muss [N_MAX] möchtet ich dieses durch eine Liste lösen.
nee.
nimm dennoch ein array. sagen wir mal der größe 64k, oder wie groß dein cache auch immer sein mag. bei mir isses zwischen 32k und 64k optimal. die arraygröße nenne ich mal B.
umd die primzahlen von 2 bis N_MAX auszueben, siebe stückweise, zuerst von 0 bis B-1, dann von B bis 2B-1, dann von 2B bis 3B-1...
um die zahlen von 2B bis 3B-1 zu sieben, brauchste nur innerhalb des (B bits großen) arrays die vielfachen der primzahlen zwischen 2 und wurzel(3B) zu streichen.
ne implemetierung hab ich auch, ist sich aber c++, denn c++ ist wegen des abstraction bonus bei sowas schneller als c oder asm.
#ifndef VHLIB_PRIMEGENERATOR_H
#define VHLIB_PRIMEGENERATOR_H
#include <iostream>
#include "BitField.h"
#include "Prime.h"
/*Mit dem Sieb des Eratostenes werden Primzahlen gesucht. Weil ich nicht den
Speicher für 2^32 Bits habe, geschieht das Abschnittweise. Dabei ist zu beachten,
daß für einen neuen Abschnitt von z1 bis z2 alle Teiler bis sqrt(z2) ausgestrichen
werden müssen. Um das zu beschleunigen, gibt es noch ein kleines Sieb, das zur
Generierung der Streich-Zahlen zuständig ist. Das wiederum wird nicht mehr von
einem Sieb versorgt, sondern nur von der Funktion nextPrime().
*/
namespace PrimeGeneratorPrivate
{
const int BIGSIZE=128*256;
const int SMALLSIZE=256;
class PrimeGeneratorLevel0
{
private:
u32 pos;
public:
u32 findFirst()
{
pos=2;
return pos;
};
u32 findNext()
{
pos=nextPrime(pos);
return pos;
};
};
class PrimeGeneratorLevel1
{
private:
enum{SIZE=SMALLSIZE};
BitField field;
u32 pos;
u32 start;
public:
PrimeGeneratorLevel1()
:field(SIZE)
{
};
u32 findFirst()
{
for(u32 t=3;t*t<=2*SIZE;t+=2)
{
for(u32 s=t+t/2;s<SIZE;s+=t)
field.set(s);
};
pos=0;
start=0;
return 2;
};
void calculateNextPage()
{
do
{
start+=2*SIZE;
field.clearAll();
PrimeGeneratorLevel0 p;
p.findFirst();
for(u32 t=p.findNext();t*t<=start+2*SIZE;t=p.findNext())
{
u32 s=(t-start%t)%t;s+=(s%2==0?t:0);s/=2;
while(s<SIZE)
{
field.set(s);
s+=t;
};
};
pos=field.findFirstFalse();
}while(pos==-1);
};
u32 findNext()
{
pos=field.findNextFalse(pos);
if(pos==-1)
{
calculateNextPage();
};
return 2*pos+start+1;
};
};
class PrimeGeneratorLevel2
{
private:
enum{SIZE=BIGSIZE};
BitField field;
u32 pos;
u32 start;
public:
PrimeGeneratorLevel2()
:field(SIZE)
{
start=-1;
};
u32 findFirst()
{
if(start!=0)
{
for(u32 t=3;t*t<=2*SIZE;t+=2)
{
for(u32 s=t+t/2;s<SIZE;s+=t)
field.set(s);
};
start=0;
};
pos=0;
return 2;
};
void calculateNextPage()
{
do
{
u32 old=start;
start+=2*SIZE;
if(start<old)
{
pos=-1;
};
field.clearAll();
PrimeGeneratorLevel1 p;
p.findFirst();
for(u32 t=p.findNext();t*t<=start+2*SIZE;t=p.findNext())
{
u32 s=(t-start%t)%t;s+=(s%2==0?t:0);s/=2;
while(s<SIZE)
{
field.set(s);
s+=t;
};
};
pos=field.findFirstFalse();
}while(pos==-1);
};
u32 findNext()
{
pos=field.findNextFalse(pos);
if(pos==-1)
{
calculateNextPage();
if(pos==-1) return pos;
};
return 2*pos+start+1;
};
};
static bool test()
{
std::cout<<"testing class PrimeGenerator:\n";
PrimeGeneratorLevel2 p;
PrimeGeneratorLevel0 ps;
u32 i=p.findFirst();
u32 is=ps.findFirst();
u16 c=0;
while(i<1000000000)
{
if(!++c) std::cout<<i32(i)<<'\r'<<std::flush;
if(i!=is)
{
std::cout<<"Fehler bei "<<i32(i)<<'\n';
return false;
};
i=p.findNext();
is=ps.findNext();
};
std::cout<<"OK\n";
return true;
};
};//namespace PrimeGeneratorPrivate
typedef PrimeGeneratorPrivate::PrimeGeneratorLevel2 PrimeGenerator;
/*
class PrimeTwinGenerator
{//nicht getestet!!!
private:
PrimeGenerator pg;
u32 pos;
public:
u32 findFirst()
{
pg.findFirst();
pg.findNext();
pos=pg.findNext();
return 3;
};
u32 findNext()
{
u32 o;
do
{
o=pos;
pos=pg.findNext();
}while(pos-o!=2);
return o;
};
};*/
#endif
benötigt wird dafür ein schnelles bitfeld
#ifndef VHLIB_BITFIELD_H
#define VHLIB_BITFIELD_H
#include "msvc_bits.h"
#include <algorithm>
#include <list>
class BitField
{
private:
u32 *data;
size_t size;
BitField(const BitField &);
BitField &operator=(const BitField &);
public:
BitField(size_t _size)
:size(_size)
{
size_t wordCount=(size+31)/32;
data=new u32[wordCount+2];
data[wordCount]=0;//such-ende-marke
data[wordCount+1]=-1;//such-ende-marke
clearAll();
};
friend void swap(BitField& a,BitField& b)
{
std::swap(a.data,b.data);
std::swap(a.size,b.size);
}
~BitField()
{
delete[] data;
};
size_t getSize() const
{
return size;
}
void clearAll()
{
size_t wordCount=(size+31)/32;
for(size_t i=0;i<wordCount;++i)
data[i]=0;
};
size_t countBitsTrue()
{
size_t wordCount=(size+31)/32;
size_t result=0;
for(size_t i=0;i<wordCount;++i)
result+=::countBitsTrue(data[i]);
return result;
};
void set(size_t pos)
{
size_t wordPos=pos/32;
size_t bitPos=pos%32;
setBitTrue(&data[wordPos],bitPos);
};
bool const operator[](size_t pos) const
{//const, um nicht aus Versehen zu setzen.
//Zum Setzen gibt's set und clear, das ist effizienter
size_t wordPos=pos/32;
size_t bitPos=pos%32;
return getBit(data[wordPos],bitPos);
}
void clear(size_t pos)
{
size_t wordPos=pos/32;
size_t bitPos=pos%32;
setBitFalse(&data[wordPos],bitPos);
};
size_t findFirstFalse()
{
return findNextFalse(-1);
};
size_t findNextFalse(size_t pos)
{
++pos;
size_t wordPos=pos/32;
size_t bitPos=pos%32;
u32 word=data[wordPos];
if(bitPos!=0)
word=(word>>bitPos) | (-1<<(32-bitPos));
while(word==-1)
{
word=data[++wordPos];
bitPos=0;
};
size_t r=findFirstBitFalse(word)+bitPos+32*wordPos;
if(r>=size) r=-1;
return r;
};
};
#endif
und ne schneller nextprime()
//Offset zur nächstmöglichen Primzahlenposition
u32 diffPrime[30]={1,6,5,4,3,2,1,4,3,2,1,2,1,4,3,2,1,2,1,4,3,2,1,6,5,4,3,2,1,2};
//Selber Offset zyklisch in Z30 erspart: if(m>30)m-=30
//Zweck: Sprünge sind auf modernen Prozessoren zu teuer.
u32 diffPrimeM[30]={1,6,5,4,3,2,1,4,3,2,1,2,1,4,3,2,1,2,1,4,3,2,1,6,5,4,3,2,1,2-30};
//1.073
u32 nextPrime(u32 last)
{
if(last<11) return firstPrimeNext[last];
u32 p=last;
u32 mp=p%30;
do
{
p+=diffPrime[mp];
mp+=diffPrimeM[mp];
}while(!isPrimeCheckFactorsAbove5(p));
return p;
};