Sieb des Atkin



  • Hallo,

    Ich versuche mich an dem Sieb des Atkin, aber irgendwie bekomme ich von 2..100 eine Zahl zu viel, die 49.

    std::vector<std::size_t> sieve_of_atkin(std::size_t max){
        assert(max >= 2);
    
        std::vector<std::size_t> remainders1{1, 13, 17, 29, 37, 41, 49, 53},
                                 remainders2{7, 19, 31, 43},
                                 remainders3{11, 23, 47, 59};
    
        std::vector<char> primes(max - 1);
    
        for(std::size_t n{2}; n <= max; ++n){
            if(!primes[n - 2]){
                if(std::find(remainders1.begin(), remainders1.end(), n % 60) != remainders1.end()){
                    for(std::size_t y{}; y <= std::sqrt(max); ++y){
                        for(std::size_t x{}; x <= std::sqrt(max); ++x){
                            if(4 * std::pow(x, 2) + std::pow(y, 2) == n)
                                primes[n - 2] = true;
                        }
                    }
                }
    
                else if(std::find(remainders2.begin(), remainders2.end(), n % 60) != remainders2.end()){
                    for(std::size_t y{}; y <= std::sqrt(max); ++y){
                        for(std::size_t x{}; x <= std::sqrt(max); ++x){
                            if(3 * std::pow(x, 2) + std::pow(y, 2) == n)
                                primes[n - 2] = true;
                        }
                    }
                }
    
                else if(std::find(remainders3.begin(), remainders3.end(), n % 60) != remainders3.end()){
                    for(std::size_t y{}; y <= std::sqrt(max); ++y){
                        for(std::size_t x{}; x <= std::sqrt(max); ++x){
                            if(x > y && 3 * std::pow(x, 2) - std::pow(y, 2) == n)
                                primes[n - 2] = true;
                        }
                    }
                }
            }
        }
    
        std::vector<std::size_t> results{2, 3, 5};
    
        for(std::size_t n{}; n != primes.size(); ++n){
            if(primes[n])
                results.push_back(n + 2);
        }
    
        return results;
    }
    

    Ich hab es nach der Anleitung von Wikipedia versucht. Die beschreiben was von x und y wobei ich mir nichts anderes vorstellen kann, als das, wie ich es angewendet habe. Kann mir jemand sagen was ich falsch mache?



  • @eigenartig

    Im Artikel ist die Rede von invertieren, du benutzt = true.



  • Außerdem eine ziemlich unperformante Implementierung:

    • du läßt jedesmal std::sqrt(max) berechnen (anstatt einmalig vor der Hauptschleife)
    • und auch n % 60 solltest du nur einmalig berechnen lassen
    • statt std::pow() solltest du einfach x * x (bzw. y * y) benutzen (keine Konvertierung zu double!)

    PS: Im englischen Wiki dazu (Sieve of Atkin) ist Pseudocode dafür beschrieben.

    PPS: Wenn du die 3 remainders-Vektoren selbst wieder in einen std::vector stecken würdest, dann bräuchtest du nicht 3x denselben Code zu schreiben (zusätzlich mit einem Vektor von Funktionen (z.B. als Lambda-Ausdrücke) für die Formeln).



  • 
    std::vector<std::size_t> sieve_of_atkin(std::size_t max){
        assert(max >= 2);
    
        std::vector<std::size_t> remainders1{1, 13, 17, 29, 37, 41, 49, 53},
                                 remainders2{7, 19, 31, 43},
                                 remainders3{11, 23, 47, 59};
    
        std::vector<char> primes(max - 1);
    
        for(std::size_t y{}; y <= std::sqrt(max); ++y){
            for(std::size_t x{}; x <= std::sqrt(max); ++x){
                for(std::size_t n{2}; n <= max; ++n){
                    if(primes[n - 2]){
                        if(std::find(remainders1.begin(), remainders1.end(), n % 60) != remainders1.end()){
                            if(4 * x*x + y*y == n)
                                primes[n - 2] = !primes[n - 2];
                        }
    
                        else if(std::find(remainders2.begin(), remainders2.end(), n % 60) != remainders2.end()){
                            if(3 * x*x + y*y == n)
                                primes[n - 2] = !primes[n - 2];
                        }
    
                        else if(std::find(remainders3.begin(), remainders3.end(), n % 60) != remainders3.end()){
                            if(x > y && 3 * x*x - y*y == n)
                                primes[n - 2] = !primes[n - 2];
                        }
                    }
                }
            }
        }
    
        std::vector<std::size_t> results{2, 3, 5};
    
        for(std::size_t n{}; n != primes.size(); ++n){
            if(primes[n]){
                for(std::size_t n2{(n+2)*(n+2)}; n2 <= max; n2 += n + 2)
                    primes[n2 - 2] = false;
    
                results.push_back(n + 2);
            }
        }
    
        return results;
    }
    

    Bekomme immer noch das gleiche Ergebnis und ist sehr langsam.



  • @eigenartig sagte in Sieb des Atkin:

    sehr langsam.

    Du sparst dir 2 std::size_t und baust dafür ständig Subtraktionen ein?
    Immer wieder modulo von selben Werten?
    Das zusammenbauen der results sieht merkwürdig aus.

    template<typename T>
    bool is_odd(T value)
    {
    	return value & 1;
    }
    
    std::deque<std::size_t> atkin(std::size_t max)
    {
    	std::deque<std::size_t> result{ 2, 3, 5 };
    	std::deque<bool> sieve(max, false);
    	
    	for (std::size_t n{ 7 }; n < sieve.size(); ++n) {
    		auto rem{ n % 60 };
    		if (!is_odd(rem))
    			continue;
    
    		if (rem == 1 || rem == 13 || rem == 17 || rem == 29 || rem == 37 || rem == 41 || rem == 49 || rem == 53) {
    			// 4x² + y² = n
    			for (std::size_t y{ 0 }; y < n; ++y)
    				for (std::size_t x{ 0 }; x < n; ++x)
    					if (4 * (x * x) + y * y == n)
    						sieve[n] = !sieve[n];			
    		}
    		else if (rem == 7 || rem == 19 || rem == 31 || rem == 43) {
    			// 3x² + y² = n
    			for (std::size_t y{ 0 }; y < n; ++y)
    				for (std::size_t x{ 0 }; x < n; ++x)
    					if (3 * (x * x) + y * y == n)
    						sieve[n] = !sieve[n];
    		}
    		else if (rem == 11 || rem == 23 || rem == 47 || rem == 59) {
    			// 3x² − y² = n | x > y
    			for (std::size_t y{ 0 }; y < n; ++y)
    				for (std::size_t x{ y + 1 }; x < n; ++x)
    					if (3 * (x * x) - y * y == n)
    						sieve[n] = !sieve[n];
    		}
    	}
    
    	for (std::size_t n{}; n < sieve.size(); ++n) {
    		if (!sieve[n])
    			continue;
    		
    		result.push_back(n);
    
    		for (auto not_prime{ n * n }; not_prime <= max; not_prime += n)
    			sieve[not_prime] = false;		
    	}
    	
    	return result;
    }
    


  • @Swordfish, ok danke.

    Dann ist es aber immer noch langsamer (bei 10000 viel langsamer) als bei meinem Sieb des Eratosthenes.



  • @eigenartig ich dachte deins rechnet falsch?



  • @Swordfish sagte in Sieb des Atkin:

    @eigenartig ich dachte deins rechnet falsch?

    Wie meinst du das? Mit deinem gehts auch nicht schneller.



  • Hab' mich verlesen. Nenn mal Code mit dem Du gemessen hast und die Wertebereiche.



  • @Swordfish, nein nicht seitdem ich die Multiplen am Ende rausnehm.



  • @Swordfish sagte in Sieb des Atkin:

    Hab' mich verlesen. Nenn mal Code mit dem Du gemessen hast und die Wertebereiche.

    Dauert ewig lang bei mir. Ideone ist da etwas schneller aber gibt mir keinen Output.

    Und bei mir lokal ist erst nach 1 - 2 Minuten fertig gewesen. Mein Sieb des Eratosthenes macht das im Millisekundenbereich.



  • Beachte den letzten Satz dazu in der englischen Wiki:

    Therefore, this special version is likely more of value as an intellectual exercise than a practical prime sieve with reduced real time expended for a given large practical sieving range.

    Die Schleifen für x und y müßten also noch sehr hinsichtlich der Durchläufe optimiert werden (wie z.B. im Wiki-Pseudocode, daß teilweise nur ungerade Werte durchlaufen werden), aber auch besonders die obere Schranke n ist viel zu ineffizient (da sie zumindestens für die ersten beiden Formeln - da nur positive Terme - deutlich nach unten korrigiert werden kann).


Log in to reply