primes - prime - sieve eratosthenes implementation

Tamiz de explicacion atkin

Question

Estoy haciendo un proyecto en este momento y necesito un método eficiente para calcular los números primos. He usado el tamiz de Eratóstenes , pero he estado buscando y he encontrado que el tamiz de Atkin es un método más eficiente. Me ha resultado difícil encontrar una explicación (que haya podido entender!) De este método. ¿Como funciona? El código de ejemplo (preferiblemente en C o python) sería brillante.

Edición: gracias por su ayuda, lo único que todavía no entiendo es a qué se refieren las variables x e y en el pseudo código. ¿Podría alguien por favor arrojar algo de luz sobre esto para mí?

Answer 1

Edición: lo único que todavía no entiendo es a qué se refieren las variables x e y en el pseudo código. ¿Podría alguien por favor arrojar algo de luz sobre esto para mí?

Para obtener una explicación básica del uso común de las variables ''x'' e ''y'' en el pseudocódigo de la página de Wikipedia (u otras implementaciones mejores del Tamiz de Atkin), es posible que encuentre útil mi respuesta a otra pregunta relacionada .

Answer 2

Aquí hay una implementación en c ++ de un tamiz de atkins que podría ayudarte ...

#include <iostream>
#include <cmath>
#include <fstream>
#include<stdio.h>
#include<conio.h>
using namespace std;

#define  limit  1000000

int root = (int)ceil(sqrt(limit));
bool sieve[limit];
int primes[(limit/2)+1];

int main (int argc, char* argv[])
{
   //Create the various different variables required
   FILE *fp=fopen("primes.txt","w");
   int insert = 2;
   primes[0] = 2;
   primes[1] = 3;
   for (int z = 0; z < limit; z++) sieve[z] = false; //Not all compilers have false as the       default boolean value
   for (int x = 1; x <= root; x++)
   {
        for (int y = 1; y <= root; y++)
        {
             //Main part of Sieve of Atkin
             int n = (4*x*x)+(y*y);
             if (n <= limit && (n % 12 == 1 || n % 12 == 5)) sieve[n] ^= true;
             n = (3*x*x)+(y*y);
             if (n <= limit && n % 12 == 7) sieve[n] ^= true;
             n = (3*x*x)-(y*y);
             if (x > y && n <= limit && n % 12 == 11) sieve[n] ^= true;
        }
   }
        //Mark all multiples of squares as non-prime
   for (int r = 5; r <= root; r++) if (sieve[r]) for (int i = r*r; i < limit; i += r*r) sieve[i] = false;
   //Add into prime array
   for (int a = 5; a < limit; a++)
   {
            if (sieve[a])
            {
                  primes[insert] = a;
                  insert++;
            }
   }
   //The following code just writes the array to a file
   for(int i=0;i<1000;i++){
             fprintf(fp,"%d",primes[i]);
             fprintf(fp,", ");
   }

   return 0;
 }

Answer 3

La página wiki es siempre un buen lugar para comenzar, ya que explica el algoritmo en su totalidad y proporciona un pseudocódigo comentado. (NB: hay muchos detalles, y dado que el sitio web de wiki es confiable, no lo citaré aquí).

Para referencias en los idiomas específicos que mencionaste:

• Implementación de C (optimizado)
• Implementación de Python

Espero que ayude.

Answer 4

El artículo de Wikipedia tiene una explicación:

• "El algoritmo ignora por completo cualquier número divisible por dos, tres o cinco. Todos los números con un resto de módulo sesenta son divisibles por dos y no son primos. Todos los números con el resto de módulo divisible por tres también son divisibles por tres y no primo. Todos los números con módulo sesenta resto divisible por cinco son divisibles por cinco y no primos. Todos estos restos se ignoran ".

Empecemos por los famosos.

primes = sieve [2..] = sieve (2:[3..])   
  -- primes are sieve of list of 2,3,4... , i.e. 2 prepended to 3,4,5...
sieve (x:xs) = x : sieve [y | y <- xs, rem y x /= 0]   -- set notation
  -- sieve of list of (x prepended to xs) is x prepended to the sieve of 
  --                  list of `y`s where y is drawn from xs and y % x /= 0

Veamos cómo procede para unos primeros pasos:

primes = sieve [2..] = sieve (2:[3..]) 
                     = 2 : sieve p2     -- list starting w/ 2, the rest is (sieve p2)
p2 = [y | y <- [3..], rem y 2 /= 0]     -- for y from 3 step 1: if y%2 /= 0: yield y

p2 no contiene múltiplos de 2 . IOW, solo contendrá números coprime con 2 ; ningún número en p2 tiene 2 como su factor. Para encontrar p2 , no necesitamos probar dividir por 2 cada número en [3..] (es decir, 3 y más, 3,4,5,6,7, ... ), porque podemos enumerar todos los múltiplos de 2 por adelantado:

rem y 2 /= 0  ===  not (ordElem y [2,4..])     -- "y is not one of 2,4,6,8,10,..."

ordElem es como elem ( es decir , la prueba del elemento ), simplemente asume que su argumento de lista es una lista ordenada y creciente de números, por lo que puede detectar la no presencia de forma segura y la presencia:

ordElem y xs = take 1 (dropWhile (< y) xs) == [y]   -- = elem y (takeWhile (<= y) xs) 

El elemento ordinario no asume nada, por lo que debe inspeccionar cada elemento de su argumento de lista, por lo que no puede manejar listas infinitas. ordElem puede. Por lo que entonces,

p2 = [y | y <- [3..], not (ordElem y [2,4..])]   -- abstract this as a function `diff`:
   = diff [3..] [2,4..]               -- diff ys xs = [y | y <- ys, not (ordElem y xs)]
   -- 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
   -- . 4 . 6 . 8 . 10  . 12  . 14  . 16  . 18  . 20  . 22  .
   = diff [3..] (map (2*) [2..])           --  y > 2, so [4,6..] is enough
   = diff [2*k+x | k <- [0..], x <- [3,4]] -- "for k from 0 step 1: for x in [3,4]:
          [2*k+x | k <- [0..], x <- [  4]] --                            yield (2*k+x)"
   = [     2*k+x | k <- [0..], x <- [3  ]]                 -- 2 = 1*2 = 2*1
   = [3,5..]                                               -- that''s 3,5,7,9,11,...

p2 es solo una lista de todas las probabilidades por encima de 2 . Bueno, eso lo sabíamos. ¿Qué pasa con el siguiente paso?

sieve p2 = sieve [3,5..] = sieve (3:[5,7..]) 
                         = 3 : sieve p3
p3 = [y | y <- [5,7..], rem y 3 /= 0]
   = [y | y <- [5,7..], not (ordElem y [3,6..])]           -- 3,6,9,12,...
   = diff [5,7..] [6,9..]         -- but, we''ve already removed the multiples of 2, (!)
   -- 5 . 7 . 9 . 11 . 13 . 15 . 17 . 19 . 21 . 23 . 25 . 27 .
   -- . 6 . . 9 . . 12  . . 15 . . 18  . . 21 . . 24  . . 27 .
   = diff [5,7..] (map (3*) [3,5..])                       -- so, [9,15..] is enough
   = diff [2*k+x | k <- [0..], x <- [5]] (map (3*)
          [2*k+x | k <- [0..], x <- [3]] )
   = diff [6*k+x | k <- [0..], x <- [5,7,9]]               -- 6 = 2*3 = 3*2
          [6*k+x | k <- [0..], x <- [    9]] 
   = [     6*k+x | k <- [0..], x <- [5,7  ]]               -- 5,7,11,13,17,19,...

Y la siguiente,

sieve p3 = sieve (5 : [6*k+x | k <- [0..], x <- [7,11]])
         = 5 : sieve p5
p5 = [y | y <-        [6*k+x | k <- [0..], x <- [7,11]], rem y 5 /= 0]
   = diff [ 6*k+x | k <- [0..], x <- [7,11]]   (map (5*)
          [ 6*k+x | k <- [0..], x <- [5, 7]]   )                    -- no mults of 2 or 3!
   = diff [30*k+x | k <- [0..], x <- [7,11,13,17,19,23,25,29,31,35]] -- 30 = 6*5 = 5*6
          [30*k+x | k <- [0..], x <- [                 25,      35]]
   = [     30*k+x | k <- [0..], x <- [7,11,13,17,19,23,   29,31   ]]

Esta es la secuencia con la que trabaja el tamiz de Atkin. No hay múltiplos de 2, 3 o 5 en él. Atkin y Bernstein trabajan el módulo 60 , es decir, duplican el rango:

p60 = [ 60*k+x | k <- [0..], x <- [1, 7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59]]

A continuación, utilizan algunos teoremas sofisticados para conocer algunas propiedades de cada uno de esos números y tratar cada uno en consecuencia. Todo se repite (a la la "rueda") con un período de 60 .

• "Todos los números n con módulo-sesenta resto 1, 13, 17, 29, 37, 41, 49 o 53 (...) son primos si y solo si el número de soluciones a 4x^2 + y^2 = n es impar y el número es squarefree ".

Qué significa eso? Si sabemos que la cantidad de soluciones a esa ecuación es impar para tal número, entonces es primo si está libre de cuadrados. Eso significa que no tiene factores repetidos (como 49, 121, etc.).

Atkin y Bernstein lo utilizan para reducir el número de operaciones en general: para cada cebado (de 7 en adelante), enumeramos (y marcamos para remover) los múltiplos de su cuadrado , por lo que están mucho más separados que en el tamiz de Eratóstenes. así que hay menos de ellos en un rango dado.

El resto de las reglas son:

• "Todos los números n con modulo sesenta resto 7, 19, 31 o 43 (...) son primos si y solo si el número de soluciones a 3x^2 + y^2 = n es impar y el número es squarefree. "

• "Todos los números n con modulo sesenta resto 11, 23, 47 o 59 (...) son primos si y solo si el número de soluciones a 3x^2 − y^2 = n es impar y el número es squarefree. "

Esto se ocupa de todos los 16 números de núcleo en la definición de p60 .

vea también: ¿Cuál es el algoritmo más rápido para encontrar números primos?

La complejidad temporal del tamiz de Eratóstenes en la producción de primos hasta N es O (N log log N) , mientras que la del tamiz de Atkin es O (N) (dejando a un lado el log log N adicional log log N optimizaciones que no se escalan bien ). La sabiduría aceptada, sin embargo, es que los factores constantes en el tamiz de Atkin son mucho más altos y, por lo tanto, solo pueden rendir muy por encima de los números de 32 bits ( N> 2 ³² ) , en todo caso .

Answer 5

// Title : Seive of Atkin ( Prime number Generator) 

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

int main()
{
    ios_base::sync_with_stdio(false);
    long long int n;
    cout<<"Enter the value of n : ";
    cin>>n;
    vector<bool> is_prime(n+1);
    for(long long int i = 5; i <= n; i++)
    {
        is_prime[i] = false;
    }
    long long int lim = ceil(sqrt(n));

    for(long long int x = 1; x <= lim; x++)
    {
        for(long long int y = 1; y <= lim; y++)
        {
            long long int num = (4*x*x+y*y);
            if(num <= n && (num % 12 == 1 || num%12 == 5))
            {
                is_prime[num] = true;
            }

            num = (3*x*x + y*y);
            if(num <= n && (num % 12 == 7))
            {
                is_prime[num] = true;
            }

            if(x > y)
            {
                num = (3*x*x - y*y);
                if(num <= n && (num % 12 == 11))
                {
                    is_prime[num] = true;
                }
            }
        }
    }
    // Eliminating the composite by seiveing
    for(long long int i = 5; i <= lim; i++)
    {
        if(is_prime[i])
            for(long long int j = i*i; j <= n; j += i)
            {
                is_prime[j] = false;
            }
    }
    // Here we will start printing of prime number
   if(n > 2)
   {
       cout<<"2/t"<<"3/t";
   }
   for(long long int i = 5; i <= n; i++)
   {
         if(is_prime[i])
         {
             cout<<i<<"/t";
         }
    }
    cout<<"/n";
    return 0;
}