El tamiz de Atkin

He estado tratando de aprender algoritmos para generar números primos y me encontré con Sieve of Atkin en Wikipedia. Entiendo casi todas las partes del algoritmo, excepto algunas. Aquí están las preguntas:

  1. ¿Cómo se forman las tres ecuaciones cuadráticas a continuación? 4x ^ 2 + y ^ 2, 3x ^ 2 + y ^ 2 y 3x ^ 2-y2
  2. El algoritmo en wikipedia habla del módulo sesenta pero no entiendo cómo / dónde se usa en el psudocódigo a continuación.
  3. ¿Cómo se encuentran estos recordatorios 1,5,7 y 11?

A continuación se muestra el pseudocódigo de Wikipedia para referencia:

// arbitrary search limit limit ← 1000000 // initialize the sieve for i in [5, limit]: is_prime(i) ← false // put in candidate primes: // integers which have an odd number of // representations by certain quadratic forms for (x, y) in [1, √limit] × [1, √limit]: n ← 4x²+y² if (n ≤ limit) and (n mod 12 = 1 or n mod 12 = 5): is_prime(n) ← ¬is_prime(n) n ← 3x²+y² if (n ≤ limit) and (n mod 12 = 7): is_prime(n) ← ¬is_prime(n) n ← 3x²-y² if (x > y) and (n ≤ limit) and (n mod 12 = 11): is_prime(n) ← ¬is_prime(n) // eliminate composites by sieving for n in [5, √limit]: if is_prime(n): // n is prime, omit multiples of its square; this is // sufficient because composites which managed to get // on the list cannot be square-free is_prime(k) ← false, k ∈ {n², 2n², 3n², ..., limit} print 2, 3 for n in [5, limit]: if is_prime(n): print n 

El pseudo código de Sieve of Atkin del artículo de Wikipedia que ha citado contiene las respuestas a sus preguntas o la discusión sobre el artículo para el que Will Ness ha proporcionado un enlace, aunque es posible que no pueda juntar la información. Las respuestas cortas son las siguientes:

  1. Las tres ecuaciones provienen de la prueba matemática de Atkin de que todos los primos ocurrirán como las soluciones de una o más de esas tres ecuaciones con las condiciones de módulo apropiadas para todos los valores válidos de las variables ‘x’ e ‘y’, y de hecho él probó que los números primos verdaderos serán aquellos números que tengan un número impar de soluciones para esas tres ecuaciones (se deja como Verdadero cuando se alterna un número impar de veces desde las condiciones iniciales de Falso) con las condiciones del módulo excluyendo aquellos números que son divisibles por cuadrados de primos encontrados menores o iguales a la raíz cuadrada del número probado.

  2. El verdadero tamiz de Atkin se basa en un conjunto de condiciones de módulo 60; este pseudo código representa una simplificación donde hay menos rangos de condiciones para cada ecuación usando un conjunto de condiciones de módulo 12 (5 × 12 = 60) – sin embargo, esto resulta en un 20% de trabajo adicional debido a la redundancia introducida en este nuevo conjunto de condiciones. También es la razón por la cual este pseudo código simplificado necesita comenzar su barrido de cebado a 5 en lugar de 7 normal y hacer las eliminaciones libres de los cuadrados de los cebadores básicos comenzando en un primo base de 5 a un costo adicional en tiempo de ejecución como los factores de 5 no se manejan de otra manera correctamente. La razón de esta simplificación es quizás sacrificar algunos gastos adicionales de código para facilitar la complejidad del código al tener que verificar los valores, aunque esto se puede hacer con una sola tabla buscar mientras que el extra más del 30% en trabajo debido al uso del módulo 12 no se puede reducir

  3. Los “mods” (deberían ser rests deletreados) los encuentran los operadores ‘mod’ en el pseudo código que citó, que representa el operador de módulo en cualquier lenguaje de computadora que se pueda usar, a menudo representado por el símbolo “%” en lenguajes de computadora tales como C, Java, C #, F #, etcétera. Este operador produce el rest entero después de una división entera del dividendo (el primero de los números, antes del ‘mod’) por el divisor (el segundo de los números, después del símbolo ‘mod’). La razón por la que los rests son solo cuatro en lugar de los 16 utilizados en el algoritmo completo del módulo 60 se debe a las simplificaciones de la reducción a un algoritmo del módulo 12. Notará que con las condiciones del módulo 12, la condición “4x” pasa 25 que normalmente sería eliminada por las condiciones del módulo 60, por lo tanto, muchos compuestos que contienen 25 como factor deben ser eliminados por el primo adicional de verificación libre de 5 cuadrados . Del mismo modo, 55 pasa el control “3x +” y 35 pasa el control “3x-” donde no lo harían para el algoritmo completo del módulo 60.

Como se discutió en la sección “Conversación” del artículo de Wikipedia y se insinuó más arriba, este pseudo código nunca es mucho más rápido que las implementaciones básicas de Sieve of Eratóstenes (SoE), solo una que usa el mismo grado de factorización de rueda debido a sus ineficiencias : las variables ‘x’ e ‘y’ no necesitan extenderse en todo el rango hasta la raíz cuadrada del rango tamizado en muchos casos (mencionado en el artículo), el uso apropiado de la rueda módulo 60 recupera las redundancias en el uso de la simplificación del módulo 12 (como mencioné anteriormente), y existen patrones de celosía del módulo en las soluciones a las ecuaciones cuadráticas tales que las condiciones que usan las operaciones del módulo (computacionalmente lentas) no necesitan ser probadas mediante el uso de un algoritmo que automáticamente omite aquellas soluciones que no satisfarían esas condiciones de módulo de acuerdo con los patrones de celosía (solo se menciona muy oscuramente en el documento completo de Atkin y Bernstein).

Para responder a una pregunta que no hizo, debería tener: “¿Por qué usar el tamiz de Atkin?” .

La razón principal por la que se implementa el Sieve of Atkin (SoA) en lugar de Sieve of Eratosthenes (SoE) es que es el “Conocimiento de Internet” que el SOA es más rápido. Hay dos razones para esta creencia, de la siguiente manera:

  1. Se supone que el SoA es más rápido porque la complejidad computacional asintótica es menor para él que para el SoE por un factor de log (log n) donde n es el rango de primos tamizados. Hablando en términos prácticos, yendo de un rango de dos a la potencia 32 (cuatro mil millones más) a dos a la potencia 64 (alrededor de 2 seguido de 19 ceros), este factor es seis sobre cinco es igual a 1.2. Eso significa que el verdadero SoE debería tomar 1.2 veces más de lo esperado por solo una relación lineal cuando se tamiza al rango numérico de 64 bits en comparación con el rango numérico de 32 bits, mientras que el SoA tendrá una relación lineal si todos fueran ideales . Puede apreciar que, en primer lugar, este es un factor muy pequeño para una gran diferencia en rangos de números y, segundo, que este beneficio solo se cumple si los dos algoritmos tienen el mismo rendimiento o casi el mismo en un rango razonable de primos.

    Lo que no se entiende claramente por ese “conocimiento de Internet” es que estas cifras solo se aplican cuando se compara la relación de rendimiento en un rango dado en comparación con otro rango dado para el mismo algoritmo , no como una comparación entre diferentes algoritmos. Por lo tanto, es inútil como una prueba de que el SoA será más rápido que el SoE ya que el SoA podría comenzar con una sobrecarga mayor para un rango de tamiz determinado de un algoritmo SoE particular como en el siguiente ejemplo optimizado de SoE.

  2. Se cree que el SoA es más rápido debido a la comparación computacional hecha y publicada por Atkin y Bernstein según su artículo enlazado en el artículo de Wikipedia. Si bien el trabajo es preciso, solo se aplica a las limitaciones artificiales que hicieron en el algoritmo SoE que compararon: dado que el algoritmo SoA se basa en la factorización de módulo 60 que representa dos rotaciones de rueda de factorización de 2,3,5, también limitaron el SoE algoritmo para esa misma factorización de rueda. Al hacer esto, el SoE realiza aproximadamente 424,000 operaciones de eliminación de número compuesto en el rango de mil millones probado, mientras que un SoA completamente optimizado como probado tiene aproximadamente 326,000 operaciones combinadas de desactivación alternante y cuadrada, que toman aproximadamente el mismo tiempo que cada operación de eliminación de números compuestos en el SoE debido a que está escrito en el mismo estilo. Esto significa que el SoA es aproximadamente un 30% más rápido que el SoE para este conjunto particular de condiciones de factorización de ruedas , y eso es exactamente lo que mostraron las pruebas de comparación de Atkins y Bernstein.

    Sin embargo, el SoA no responde a niveles adicionales de factorización de las ruedas ya que el nivel 2,3,5 está “integrado” al algoritmo, mientras que el SoE sí responde a niveles adicionales: usando una rueda 2,3,5,7 la factorización da como resultado aproximadamente el mismo número de operaciones, lo que significa aproximadamente el mismo rendimiento para cada una. Se puede usar incluso un nivel de factorización de ruedas parciales mayor que el nivel 2,3,5,7 para obtener el número de operaciones para SoE aproximadamente un 16,7% menos que SoA, para un mejor rendimiento proporcional. Las optimizaciones requeridas para implementar estos niveles adicionales de factorización de ruedas son en realidad más simples que la complejidad del código para implementar el SoA original optimizado. La huella de memoria para implementar estas implementaciones comparables segmentadas por páginas es casi la misma, ya que es el tamaño de los búferes de página más el conjunto de primos base.

    Además, ambos se beneficiarían de estar escritos en un estilo de “búsqueda del estado de la máquina” que ayudaría a una mejor optimización utilizando código en línea y desenrollado extremo del bucle, pero el SoE es más apropiado que el SoA debido a que es un algoritmo más simple.

Por lo tanto, para tamices oscila hasta aproximadamente el rango numérico de 32 bits, el SoE optimizado al máximo es aproximadamente un 20% más rápido (o más con factorización de rueda adicional) que el SoA; sin embargo, el SoA tiene esta ventaja de complejidad computacional asintótica, por lo que habrá un punto en el que se pondrá al día. Ese punto estará aproximadamente en el rango donde la razón de log (log n) a log (log (2 ^ 32)) o 5 es 1.2 o un rango de aproximadamente 2 veces diez a la potencia diecinueve, un número excesivamente grande. Si la ejecución de SoA optimizada en una computadora de escritorio moderna tomara aproximadamente un tercio de segundo para calcular los números primos en el rango numérico de 32 bits, y si la implementación fuera ideal, como al tomar un aumento de tiempo lineal con un rango creciente, tomaría aproximadamente 45 años para calcular los números primos a este rango. Sin embargo, el análisis de los números primos en estos rangos altos se hace a menudo en trozos pequeños, para los cuales el uso del SoA sería teóricamente práctico en comparación con el SoE para tamices de gran número, pero por un factor muy pequeño.

EDIT_ADD: en realidad, ni la página segmentada SoE ni la SoA continúan ejecutándose en tiempo lineal para estos enormes rangos en comparación con los rangos más bajos, ya que ambos tienen problemas con las operaciones de “alternancia” y “descarte” que pierden eficiencia debido a tener saltear grandes cantidades de páginas por salto; esto significa que ambos requerirán algoritmos modificados para manejar este “salto de página” y la pequeña ventaja del SoA puede borrarse por completo si hay pequeñas diferencias en la forma en que se hace esto. El SoA tiene muchos más términos en las “tablas de salto” que el SoE por aproximadamente la razón inversa entre el número de primos encontrados hasta la raíz cuadrada del rango en ese rango, y esto probablemente agregará una O (log n) término para ambos en el procesamiento, pero un aumento constante del factor constante para el SoA que tiene un mayor número de entradas en la “tabla de salto”. Este hecho adicional anulará completamente cualquier ventaja del SoA sobre el SoE, incluso para rangos extremadamente grandes. Además, el SoA tiene una sobrecarga constante de más bucles individuales requeridos para muchos más casos al implementar las condiciones para las tres ecuaciones cuadráticas separadas más el bucle “primo libre de cuadrados” en lugar de solo un bucle de selección de primos, por lo que nunca puede tener un valor tan bajo tiempo de cálculo promedio por operación como el SoE cuando está completamente optimizado. END_EDIT_ADD

EDIT_ADD2: En mi opinión, gran parte de la confusión sobre el Tamiz de Atkin se debe a las deficiencias del pseudo código del artículo de Wikipedia citado en la pregunta, por lo que se ha propuesto una versión algo mejor del pseudo código que aborda al menos algunas de las deficiencias que tienen que ver con el rango de las variables ‘x’ e ‘y’ y la confusión módulo 12 versus módulo 60 de la siguiente manera:

 limit ← 1000000000 // arbitrary search limit // Initialize the sieve for i in {7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59,61,...}: is_prime(i) ← false // Put in candidate primes: // integers which have an odd number of // representations by certain quadratic forms. while n ≤ limit, n ← 4x²+y² where x ∈ {1,2,...} and y ∈ {1,3,...} // odd y's if n mod 60 ∈ {1,13,17,29,37,41,49,53}: is_prime(n) ← ¬is_prime(n) while n ≤ limit, n ← 3x²+y² where x ∈ {1,3,...} and y ∈ {2,4,...} // only odd if n mod 60 ∈ {7,19,31,43}: // x's and even y's is_prime(n) ← ¬is_prime(n) while n ≤ limit, n ← 3x²-y² where x ∈ {2,3,...} and y ∈ {x-1,x-3,...,1} //all if n mod 60 ∈ {11,23,47,59}: // even/odd odd/even combos is_prime(n) ← ¬is_prime(n) // Eliminate composites by sieving, only for those occurrences on the wheel for n² ≤ limit where n ∈ {7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59,61,...}: if is_prime(n): // n is prime, omit multiples of its square; this is // sufficient because composites which managed to get // on the list cannot be square-free while c ≤ limit, c ← n² × k where k ∈ {1,7,11,13,17,19,23,29, 31,37,41,43,47,49,53,59,...}: is_prime(c) ← false output 2, 3, 5 for n ≤ limit when n ← 60 × k + x where k ∈ {0..} and x ∈ {7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59,61}: if is_prime(n): output n 

Lo anterior parece bastante simple y es un algoritmo bastante bueno, excepto que todavía no es más rápido que un Tamiz básico de Eratóstenes que usa la misma rueda de factorización 2; 3; 5 porque desperdicia casi la mitad de sus operaciones de conmutación de bucle interno que fallan. pruebas de módulo Demostrar:

A continuación está el patrón de repetición 4x² + y² de los módulos de calificación 60 en un rango de 15 valores de ‘x’ (cada valor) verticalmente y 15 valores impares de ‘y’ horizontalmente; ambos comenzando en uno. Tenga en cuenta que son simétricos con respecto a un eje vertical, pero solo 128 de 225 o 256 de 450 para un rango completo de 30 números de valores ‘x’ son operaciones válidas de alternancia:

  0 13 29 53 0 0 53 49 53 0 0 53 29 13 0 17 0 41 0 37 17 0 1 0 17 37 0 41 0 17 37 0 1 0 0 37 0 0 0 37 0 0 1 0 37 0 13 29 53 0 0 53 49 53 0 0 53 29 13 0 41 49 0 29 1 41 29 0 29 41 1 29 0 49 41 0 0 49 13 0 0 13 0 13 0 0 13 49 0 0 17 0 41 0 37 17 0 1 0 17 37 0 41 0 17 17 0 41 0 37 17 0 1 0 17 37 0 41 0 17 0 0 49 13 0 0 13 0 13 0 0 13 49 0 0 41 49 0 29 1 41 29 0 29 41 1 29 0 49 41 0 13 29 53 0 0 53 49 53 0 0 53 29 13 0 37 0 1 0 0 37 0 0 0 37 0 0 1 0 37 17 0 41 0 37 17 0 1 0 17 37 0 41 0 17 0 13 29 53 0 0 53 49 53 0 0 53 29 13 0 1 0 0 49 0 1 49 0 49 1 0 49 0 0 1 

A continuación se muestra el patrón de repetición 3x² + y² de los módulos de calificación 60 en un rango de 5 valores impares de ‘x’ verticalmente y 15 valores pares de ‘y’ horizontalmente. Observe que son simétricos con respecto a un eje horizontal, pero solo 48 de 75 o 144 de 450 para un rango completo de 30 números de valores ‘x’ son operaciones de alternar válidas, ya que todas las 144 de 450 operaciones inválidas con ‘x’ e impar ‘y’ ya ha sido eliminado:

  7 19 0 7 43 0 19 19 0 43 7 0 19 7 0 31 43 0 31 7 0 43 43 0 7 31 0 43 31 0 19 31 0 19 0 0 31 31 0 0 19 0 31 19 0 31 43 0 31 7 0 43 43 0 7 31 0 43 31 0 7 19 0 7 43 0 19 19 0 43 7 0 19 7 0 

A continuación se muestra el patrón repetitivo de 3x²-y² de los módulos de calificación 60 en un rango de 5 valores impares de ‘x’ verticalmente y 15 valores pares de ‘y’ horizontalmente. Observe que son simétricos con respecto a un eje horizontal, pero solo 48 de 75 o 224 de 450 para un rango completo de 30 números de valores ‘x’ son operaciones válidas de alternancia:

  59 47 0 59 23 0 47 47 0 23 59 0 47 59 0 23 11 0 23 47 0 11 11 0 47 23 0 11 23 0 11 59 0 11 0 0 59 59 0 0 11 0 59 11 0 23 11 0 23 47 0 11 11 0 47 23 0 11 23 0 59 47 0 59 23 0 47 47 0 23 59 0 47 59 0 

A continuación se muestra el patrón repetitivo de 3x²-y² de los módulos de calificación 60 en un rango de 5 valores pares de ‘x’ verticalmente y 15 valores impares de ‘y’ horizontalmente. Tenga en cuenta que son simétricos con respecto a un eje vertical, pero solo 48 de 75 o 224 de 450 para un rango completo de 30 valores ‘x’ son operaciones válidas de alternancia:

  11 0 47 23 0 11 23 0 23 11 0 23 47 0 11 47 0 23 59 0 47 59 0 59 47 0 59 23 0 47 47 0 23 59 0 47 59 0 59 47 0 59 23 0 47 11 0 47 23 0 11 23 0 23 11 0 23 47 0 11 59 0 0 11 0 59 11 0 11 59 0 11 0 0 59 

Como se puede determinar mediante la inspección de las tablas anteriores, para el pseudo código como arriba hay un rango repetitivo general de 30 valores de x (incluyendo ambos pares y probabilidades) que solo tiene 688 operaciones válidas de un total de 1125 operaciones combinadas; por lo tanto, desperdicia casi la mitad de su procesamiento al saltar saltos condicionalmente debido a que las operaciones de salto no productivas forman parte de los bucles más internos. Existen dos métodos posibles para eliminar la redundancia de “golpe” de forma ineficiente, como se muestra a continuación:

  1. El método descrito en el trabajo de Atkin y Bernstein, que utiliza los patrones reconocidos en los módulos combinados de ‘x’ e ‘y’ para procesar solo el módulo “aciertos” utilizando el hecho de que una vez que uno localiza un módulo dado en el patrón, luego hay una secuencia infinita de valores en ese módulo (igual a la posición del bit de rueda) donde cada patrón está separado por un desplazamiento fácilmente calculado (variable) que tiene la forma “posición actual más A multiplicado por el valor actual de ‘x’ más B “y” posición actual más C multiplicado por el valor actual de ‘y’ más D “donde A, B, C y D son constantes simples (significado simple pequeño fácilmente manipulable). Bernstein usó ese método con la ayuda adicional de muchas tablas de búsqueda (LUT).

  2. El método para crear Tablas de búsqueda de estado de patrón (LUT) generales (una para cada una de las cuatro tablas anteriores más una para el bucle libre de primo menor) indexadas por los valores actuales de módulo de ‘x’ combinados con el módulo de ‘y’ para encontrar los parámetros A, B, C y D para omitir, no para el siguiente patrón, sino solo para la siguiente posición en la secuencia horizontal. Es probable que sea la opción de mayor rendimiento, ya que permite un ciclo de reloj extremo por reducción de operación usando el forro interno del código de bucle desenrollado y producirá un código más eficiente para rangos grandes cuando la segmentación de página es más pequeña en promedio. . Esto probablemente hará que los ciclos de reloj por bucle sean parecidos a los de un tamiz altamente optimizado de Eratóstenes como en primesieve , pero es probable que no lleguen a ese nivel debido a tener que calcular los tamaños de pasos variables en lugar de poder usar desplazamientos primarios fijos. para SoE.

Por lo tanto, hay tres objectives de gobierno en la reducción del tiempo de ejecución de un tamiz principal, de la siguiente manera:

  1. Un tamiz exitoso reduce la cantidad de operaciones, que incluso el SoA “optimizado para golpear” falla en comparación con un SoE con factor de rueda alto en aproximadamente un 16,7% para rangos de unos pocos miles de millones.

  2. Un tamiz exitoso reduce los ciclos de reloj de la CPU por operación, lo que SoA falla en comparación con un SoE altamente optimizado como primesieve porque sus operaciones son más complejas debido a los incrementos variables, de nuevo probablemente en un 20% a 60%. El primógeno (SoA) de Atkin y Bernstein toma alrededor de 4.5 en comparación con aproximadamente 2.5 ciclos de reloj por operación para primesieve (SoE) para un rango de mil millones para cada uno, pero quizás el SoA podría acelerarse un tanto tomando prestadas algunas de las técnicas de optimización de primesieve.

  3. Un tamiz exitoso tiene una complejidad de código razonablemente baja para que pueda extenderse más fácilmente a rangos más grandes utilizando técnicas como el “tamizado de cubetas” y otras optimizaciones de segmentación de páginas. Por esto, el Tamiz de Atkin falla miserablemente, ya que se vuelve exponencialmente más complejo para la página segmentando grandes rangos numéricos. Es extremadamente complejo escribir un progtwig SoA que compita incluso con el primógeno (SoA) de Bernstein y mucho menos con el primesieve, mientras que es bastante fácil escribir código que se acerque al mismo rendimiento que el primesieve.

  4. Un tamiz exitoso tiene una complejidad empírica menor, que el SoA tiene sobre el SoE por un factor de (log (log n) donde n es el rango superior para ser tamizado, pero esa proporción extra pequeña no es suficiente para compensar para las dos razones de pérdidas combinadas superiores, ya que este factor se reduce a medida que aumenta el rango. END_EDIT_ADD2

Entonces, la respuesta a la pregunta “¿Por qué usar el tamiz de Atkin?” es “No hay ninguna razón para usarlo en absoluto si el SoE se implementa con optimizaciones máximas de factorización de rueda hasta que los números tamizados son extremadamente grandes (rango numérico de 64 bits y superior) y entonces la ventaja SoA es muy pequeña y quizás no realizable en todo depende de ajustes muy pequeños en optimizaciones relativas “. .

Como respuesta a otra pregunta similar de Sieve of Atkin, publiqué una versión segmentada de C de la página de SoE optimizada según lo anterior en: https://stackoverflow.com/a/20559687/549617 .

En mi respuesta original a esta pregunta , con el fin de ayudar a una mejor comprensión del algoritmo Sieve of Atkin extendí el pseudo código Wikipedia Sieve of Atkin (SoA) para corregir algunas de las deficiencias en el algoritmo simplificado “módulo 12” dado, que muchos encuentran confuso en comparación con el algoritmo completo “módulo 60”. Además, el algoritmo original de pseudo código de Wikipedia conduce a implementaciones ineficientes; específicamente, los rangos de las variables ‘x’ e ‘y’ son cada uno más anchos de lo necesario y la simplificación del módulo 12 da como resultado un 20% adicional de trabajo redundante. Sin embargo, como se discutió en esa respuesta, ese pseudo código aún no es eficiente ya que las pruebas de módulo todavía están en los bucles de solución cuadrática más interna y por lo tanto casi la mitad de esos bucles de solución son improductivos para aquellas soluciones que no pasan las pruebas de módulo; esto significa que un progtwig que implemente ese pseudo código probablemente sea casi el doble de lento de lo necesario incluso cuando no haya un cuello de botella de velocidad de acceso a la memoria.

El algoritmo utilizado por Atkin y Bernstein en su progtwig “primigenio” para evitar que la redundancia de bucle se pueda entender mejor mediante una progresión a través de las siguientes observaciones basadas en las tablas de “aciertos” en mi respuesta anterior:

  1. La simetría de las tablas “aciertos” depende de si ‘y’ es impar (simetría vertical) o incluso (simetría horizontal).

  2. Por lo tanto, si quisiéramos que el avance de los patrones avanzara en la dirección horizontal como se muestra allí, podemos invertir el orden de los bucles ‘x’ e ‘y’ para “3x +” y “3x-; impar x -> incluso y “casos que significan casos ‘x’ impares.

  3. Ahora es fácil eliminar los casos cero para las dos tablas volteadas “3x” que solo tienen cinco repeticiones horizontalmente por una condición en un bucle externo como aquellos casos donde “y mod 3 es igual a 0″ más los tres casos adicionales solo para el medio ” eje “columna donde” y mod 5 es igual a 0 “.

  4. Para el último caso “3x-; incluso x -> impar y”, es fácil eliminar la mayoría de los ceros filtrando aquellos casos en los que “(y + 2) mod 3 es igual a 0” más los dos casos adicionales solo para el última línea donde “(x mod 60) es igual a 59” donde “(y + 2) mod 5 es igual a 0” (eliminaciones duplicadas para la columna de eje vertical simétrica, que siempre se elimina). Aunque estas “pruebas y” parecerían tener lugar en un bucle interno, ya que solo hay cinco “visitas” por cada lado del eje de simetría vertical, estas pueden codificarse en línea por el bucle ‘x’.

  5. Las eliminaciones más difíciles son para la tabla “4x” donde el patrón es más complejo: hay cinco patrones horizontales distintos donde cada uno se puede reconocer porque “(x mod 60) para la primera columna” es distinto, con uno de los principales patrones de cero filas en la primera tabla de arriba que tienen un módulo de 5 y el otro un módulo de 25; ahora los cinco patrones diferentes se pueden determinar “delimitando” los casos como desplazamientos a cada lado del eje simétrico vertical para cada patrón.

  6. El espaciado a la siguiente fila es exactamente 4 × ((x + 1) ² – x²) o 8x + 4 para la primera tabla, (y + 2) ² – y² o 4y + 4 para las dos nuevas tablas (intercambiadas) nuevas , y 3 × ((x + 2) ² – x²) o 12x + 12 para la última tabla.

  7. Para todas las tablas verticalmente simétricas (ahora), el desplazamiento del eje central se puede determinar a partir del desplazamiento de la primera columna como (y + 16) ² – y² o 32y + 256 para las filas largas y 3 × ((x + 4) ² – x²) o 24x + 48 para las filas cortas.

  8. De forma similar, los desplazamientos simétricos en cualquier lado del eje vertical se pueden calcular fácilmente a partir del desplazamiento del eje vertical, por ejemplo, mediante el par (y + 2) ² – y² o 4 + 4x con (y – 2) ² – y² o 4 – 4y donde y es el desplazamiento del eje central para los casos de fila larga y más y menos dos posiciones; el caso de los casos de la fila corta “x” es similar, pero solo se multiplica por tres para la escala general de los valores “3x” ‘x’.

  9. Los filtros anteriores determinados por las condiciones de variables horizontales parecen romper nuestro objective de eliminar el código condicional en los bucles internos, pero esa regla no se aplica al determinar que el patrón no es el bucle interno sino que comienza una serie completa de bucles interiores nuevos en todos los patrones horizontalmente en la matriz para cada elemento válido del patrón. Esto funciona porque podemos demostrar matemáticamente que cada una de las mismas posiciones equivalentes en el siguiente patrón (que tiene el mismo módulo) se separan de la siguiente manera: para el patrón horizontal avanza 30 en los casos de filas largas, los patrones están separados por (x + 30) ² – x² o 60x + 900 entonces 60 × (x + 15); para el patrón horizontal avanza 10 en los casos de filas cortas, los patrones están separados por 3 × ((x + 10) ² – x²) así que 60 * (x + 5). Como ambos tienen un módulo 60 de cero, esta es la explicación exacta de por qué hay patrones repetitivos del mismo módulo en cada posición. Por lo tanto, una vez que encontramos el desplazamiento para cada posición de esquina superior izquierda, podemos determinar el módulo y desplazamiento para las posiciones de inicio para las posiciones de patrón diferentes de cero y simplemente enrollar un bucle interno usando las relaciones de patrón desde ese punto hasta el límite de la matriz de tamizado .

  10. Todos los incrementos de columna y fila anteriores se pueden realizar sin utilizar una operación “dividir” usando módulo adicional donde el algoritmo realiza un seguimiento de los pares de cocientes y el módulo de 60 con solo una operación de “control y ajuste de desbordamiento” en el módulo / cociente par después de cada operación, pero el código condicional para hacer el “verificar y ajustar” es probablemente más costoso desde el punto de vista de la forma escrita, para lo cual la mayoría de los comstackdores de optimización generarán una operación múltiple usando características de truncamiento de desbordamiento para reemplazar la operación de división computacionalmente costosa para la división por enteros pequeños, que generalmente toma menos tiempo que la combinación de operaciones más el código condicional requerido para el método “verificar y ajustar” o usar una operación de división directa.

  11. Los conjuntos de eliminaciones anteriores funcionan muy bien para una representación de bits empaquetados de los primos candidatos en que cada uno de los 16 valores de módulo válidos se puede representar como un bit en una palabra de 16 bits y el índice de palabra presenta los índices de cada módulo 60 rueda; por lo tanto, podemos avanzar en incrementos uniformes uniformemente divisibles por 60 simplemente avanzando un número entero de palabras de 16 bits.

Para usar como un algoritmo segmentado no de página como para el pseudo código aquí, es suficiente mover las comprobaciones de módulo fuera del bucle más interno sin eliminar todas las “fallas de golpe” para el bucle medio resultante, ya que aunque todavía no son lo máximo eficientes en que el ciclo medio todavía procesa “fallos de éxito”, que no agrega más de un uno por ciento al tiempo de ejecución; sin embargo, para las implementaciones paginadas, “hits perdidos” incluso en los bucles externos pueden agregar mucho a la sobrecarga computacional ya que estos bucles se ejecutan una vez por página sobre el rango de valores ‘x’ para esa página, así que hay un factor sobre el logaritmo de la raíz cuadrada del rango tamizado más de los bucles externos para el SoA en comparación con el Tamiz de Eratóstenes (SoE). Una versión de pseudo código más complejo que represente una versión (opcionalmente empaquetada en bits) no segmentada por páginas del método Atkin y Bernstein “optimizado para resultados” se puede escribir de la siguiente manera:

 // arbitrary search limit limit ← 1000000000 // the set of base wheel prime candidates s ∈ {7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59,61} // Initialize the sieve as an array of wheels with // enough wheels to include the representation for limit for n ← 60 * w + x where w ∈ {0,...(limit - 7) ÷ 60}, x in s: sieve(n) ← false // start with no primes. // Put in candidate primes: // integers which have an odd number of // representations by certain quadratic forms. while n ≤ limit when n ← 4x²+y² // all x's and only odd y's where x ∈ {1,...}, y ∈ {1,31,...}: // skip by pattern width in y for cb ≤ limit when midy ← y + i, cb ← n - y² + midy² where i ∈ {0,2,...28}: // middle level loop to "hit" test modulos cbd ← (cb - 7) ÷ 60, cm ← cb mod 60 if cm in {1,13,17,29,37,41,49,53}: while c ≤ limit when c ← 60 × (cbd - midy² + ((midy + 15 × j) × j)²) + cm where j {0,...}: sieve(c) ← ¬sieve(c) // toggles the affected bit while n ≤ limit when n ← 3x²+y² // only odd x's and even y's where x ∈ {1,3,...}, y ∈ {2,32,...}: // skip by pattern width in y for cb ≤ limit when midy ← y + i, cb ← n - y² + midy² where i ∈ {0,2,...28}: // middle level loop to "hit" test modulos cbd ← (cb - 7) ÷ 60, cm ← cb mod 60 if cm in {7,19,31,43}: while c ≤ limit when c ← 60 × (cbd - midy² + ((midy + 15 × j) × j)²) + cm where j {0,...}: sieve(c) ← ¬sieve(c) // toggles the affected bit while n ≤ limit when n ← 3x²-y² //even/odd & odd/even combos where x ∈ {2,3,...}, y ∈ {x-1,x-31,...,1}: // skip by pattern width in y for midy ≥ 1 and cb ≤ limit when midy ← y - i, cb ← n + y² - midy² where i ∈ {0,2,...28}: // middle level loop to "hit" test modulos cbd ← (cb - 7) ÷ 60, cm ← cb mod 60 if cm in {11,23,47,59}: while yn ≥ 1 and c ≤ limit when yn = midy - 30 * j and c ← 60 × (cbd + midy² - ((midy - 15 × j) × j)²) + cm where j {0,...}: sieve(c) ← ¬sieve(c) // toggles the affected bit // Eliminate prime squares by sieving, only for those occurrences on the wheel for sqr ≤ limit where sqr ← n², n in {7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59,61,...}: if is_prime(n): // n is prime, omit multiples of its square; this is // sufficient because composites which managed to get // on the list cannot be square-free if c0 ≤ limit where c0 ← sqr × (m mod 60) where m in s: c0d ← (c0 - 7) ÷ 60, cm ← (c0 - 7) mod 60 + 7 // means cm in s while c ≤ limit for c ← 60 × (c0d + sqr × j) + cm where j ∈ {0,...}: sieve(c) ← false // cull composites of this prime square output 2, 3, 5 for n ≤ limit when n ← 60 × k + x where k ∈ {0..}, x in s: if sieve(n): output n 

In particular, notice how the innermost quadratic loops are very simple but increment by a linearly increasing variable amount per loop, whereas the prime squares free elimination, which uses a type of Eratosthenes progression, increments by a fixed amount “sqr”. These simple loops will execute quite quickly on a modern CPU, taking about 4 CPU clock cycles per loop if memory access is fast enough as for small sieve ranges or as in using a page segmented version where the pages fit within the L2/L1 CPU caches that have memory access times of about this range.

When the array is large as in many Megabytes, then the average memory access time is likely to average something between 20 and over a 100 CPU clock cycles (depending on the speed of RAM, efficiency of memory access by the CPU, and efficiency of use of intermediate caches) and this will be the main limiting factor in execution speed rather than the tightness of the loops other than the number of toggle/cull operations that the algorithm requires. These particular algorithms, both SoA and SoE, put a particular strain on the memory interface for large buffers in that they skip through the buffer incrementing by whole wheel circumferences times a multiplier at a time and repeat for a different offsets rather than handling all the “hits” of the wheel in sequence; this is due to moving the “hit” tests into a middle loop, which saves “hit” testing but has this larger buffer “jump” as a consequence, with the average “jump” for the SoA higher than for the SoE. Thus, an estimate of the execution time for this algorithm over the range of numbers to a billion is about 60 CPU clock cycles per operation (which will be less by a factor for the SoE due to a smaller average “jump”) multiplied by the number of toggle/cull operations (about 260,000 for the Sieve of Atkin) divided by the CPU clock speed. Therefore, a 3.5 GHz CPU will execute this algorithm over a range to a billion in about four seconds and the main difference between execution speed for non-paged implementations of the SoA and the SoE in this style is the difference in memory access efficiency for large buffers, where a highly wheel factorized SoE is more efficient by about a factor of two.

Given that for large buffers the memory access time is the limiting speed factor, it becomes very important to find an algorithm with a minimum number of operations. The SoA as implemented by Atkin and Bernstein in their “primegen” sample program takes about a constant factor of “hit” operations as a ratio to the range sieved, and by implementing even a simple version of a program from my first answer we can determine that this factor is about 0.26 of the prime range sieved: meaning that there are about 260,000 operations for a sieve range of a billion.

For a SoE implemented with a high wheel factorization of a 210 element wheel with 48 prime candidate “hits” per rotation (small wheel for primes of 2;3;5;7) in combination with a pre-cull by the additional primes of 11;13;17;19, the number of cull operations for a sieve range n is approximately n times (ln((ln n) / 2) + M – 1/(ln n) – 1/2 – 1/3 – 1/5 – 1/7 – 1/11 – 1/13 – 1/17 – 1/19) * 48 / 210 where ‘*’ means multiplied by, ‘/’ means divided by, M is the Meissel-Mertens constant M ≈ 0.2614972128… correcting the number of operations arising from the (ln ln n) estimate, the reciprocal “ln n” term is the adjustment for the culling starting at the square of the base primes (less than the square root of the range n); for n of one billion, the factor compared to n is about 0.250485568 or about 250 million cull operations for a range to a billion.

This is about the same number of operations as for the SoA and contrasts with the simple 2;3;5 wheel used in the Atkin and Bernstein comparison test, which has a factor of 0.404805008 or about 400,000 operations for a sieve range of one billion. This better ratio for the more highly factorized SoE means that the constant factors for this quite realizable implementation are about the same as for the SoA except for any differences in constant factors due to loop complexity.

EDIT_ADD1: For comparison, the following is equivalent pseudo code for a highly factorized and pre-culled SoE written in the same style as the above SoA algorithm:

 // arbitrary search limit limit ← 1000000000 // the set of base wheel primes r ∈ {23,29,31,37,41,43,47,53, 59,61,67,71,73,79,83, //positions + 19 89,97,101,103,107,109,113,121,127, 131,137,139,143,149,151,157,163, 167,169,173,179,181,187,191,193, 197,199,209,211,221,223,227,229} // an array of length 11 times 13 times 17 times 19 = 46189 wheels initialized // so that it doesn't contain multiples of the large wheel primes for n where n ← 210 × w + x where w ∈ {0,...46189}, x in r: // already if (n mod cp) not equal to 0 where cp ∈ {11,13,17,19}: // no 2,3,5,7 mstr(n) ← true else mstr(n) ← false // factors // Initialize the sieve as an array of the smaller wheels with // enough wheels to include the representation for limit for n where n ← 210 × w + x, w ∈ {0,...(limit - 19) ÷ 210}, x in r: sieve(n) ← mstr(n mod (210 × 46189)) // init pre-culled primes. // Eliminate composites by sieving, only for those occurrences on the // wheel using wheel factorization version of the Sieve of Eratosthenes for n² ≤ limit when n ← 210 × k + x where k ∈ {0..}, x in r if sieve(n): // n is prime, cull its multiples s ← n² - n × (x - 23) // zero'th modulo cull start position while c0 ≤ limit when c0 ← s + n × m where m in r: c0d ← (c0 - 23) ÷ 210, cm ← (c0 - 23) mod 210 + 23 //means cm in r while c ≤ limit for c ← 210 × (c0d + n × j) + cm where j ∈ {0,...}: sieve(c) ← false // cull composites of this prime output 2, 3, 5, 7, 11, 13, 17, 19, for n ≤ limit when n ← 210 × k + x where k ∈ {0..}, x in r: if sieve(n): output n 

Notice how much less complex this code is in spite of it being more highly wheel factorized so that there is less of a constant overhead of operations than the SoA code while the innermost culling loop is actually slightly simpler and therefore (if anything) likely slightly faster. Also notice how many more outer loops must be run for the SoA code in there are a total of very approximately the same number of outer loops as the square root of the sieving range for the SoA where there are only about the sieve range divided by the natural logarithm of the square root of that range outer loops for the SoE – a factor of about ten times less for the SoE compared to the SoA for sieve ranges of one billion. This is the trade-off that the SoA makes: more outer loops for less operations per inner loop with bigger buffer “jumps” between operations; however, this is also what makes the SoA less efficient per operation, especially for larger sieve ranges. END_EDIT_ADD1

As to the code complexity for larger ranges, these are always implemented as page segmented sieves in order to avoid the necessity of having the whole sieve array in RAM memory (memory constraints would limit the range) and in order to take advantage of CPU caching locality to reduce average memory access times as discussed above. For the SoE, this can easily be implemented based on a saved copy of the a base primes representation where there are approximately the square root of the range divided by log of that square root number of that range; for the SoA one needs to refer to that same base primes representation for the prime squares free elimination but as well needs to use values of sequence offsets plus saved values of ‘x’ or ‘y’ (or save the current values of both ‘x’ and ‘y’) in order to resume at the next page offset with a total number of saved sequences as about the square root of the range (not divided by natural logarithm of root n), or alternatively new ‘x’/’y’ pair offsets can be calculated for the start of each new page at a considerable computational overhead. The first method adds greatly to the memory use overhead, especially for multi-threaded implementations where a copy needs to be kept for each thread; the second doesn’t use more memory but uses more computational power as compared to the SoE for these larger ranges.

Page segmented algorithms are also important in order to take full advantage of modern multi-core CPU’s which can reduce the clock time required for a given task by the ratio of the number of cores used, especially for page divided algorithms where each page processing can be assigned to a separate thread.

The way that the “primesieve” SoE implementation manages average operation times of about 2.5 CPU clock cycles as compared to the about 4.5 CPU clock cycles for Atkin and Bernstein’s “primegen” sample implementation of the SoA is by extreme in-lining of loop unrolling with page segmentation; this technique can also be applied to the SoA but this particular pseudo code is unsuitable for doing that and an alternate technique of “hit rate” optimization that handles all of the modulos in a single loop (per quadratic equation loop) needs to be used; however, that is a complex algorithm that would seem to be beyond the scope of this question. Suffice it to say that the jumps in complexity so far are nothing as compared to what is required to add this extra optimization to the SoA, and even with that the ratio of outer loops still makes the SoE much more efficient.

Thus, while it is fairly easy to implement a naive SoA algorithm as in the original Wikipedia pseudo code or even using the more complex SoA pseudo code as in these two answers, it is extremely complex to write something that would compete with even Bernstein’s “primegen” (SoA) let alone with the faster “primesieve” (SoE), whereas it is quite easy to write code that comes close to the same performance as “primegen” and just a little more work to write something that competes with “primesieve” using the SoE algorithm with the slight modifications to support page segmentation.

EDIT_ADD2: To put this theory into practical terms, I would normally write an example program to demonstrate these principles; however, in this case we have all we really need in the “primegen” source code , which is the reference implementation by the inventors of the SoA. On my machine (i7-2600K, 3.5 GHz) with the buffer size adjusted to 32768 from 8192 to reflect my L1 cache size, this “primespeed” sieves the primes to one billion in about 0.365 seconds. In fact, Atkin and Bernstein cheated a little in their comparison to their Eratosthenes sieve algorithm in that they only assigned about four Kilobyte buffers to it; adjusting that buffer size to about the same size (the “#define B32 1001” to “#define B32 8008”) results in an execution time of about 0.42 seconds for “eratspeed” to about one billion, leaving it still slightly slower. However, as discussed above, it is doing more work in that it is doing about 0.4048 billion cull operations whereas the SoA is only doing about 0.26 billion combined toggle/cull operations. Using the algorithm from the above SoE pseudo code to change this is have maximum wheel factorization means that the SoE would actually do slightly less operation than the SoA at about 0.2505 billion, meaning the the execution time would be reduced to about two thirds or more or about 0.265 to 0.3 seconds making it faster than “primegen”.

Meanwhile, if one compares the “primespeed” to the “eratspeed” code one sees that the paged SoE is much less complex than the SoA code, just as indicated by the above pseudo code.

Of course, one could counter that the SoA maintains the ratio of 0.26 times the sieving range operations to as high a sieving range as one wants whereas the SoE ratio climbs as per the equations given above – by about an extra 20% for a sieving range even to one hundred billion. However, given a slight tweak to both algorithms to increase buffer size and “sub slice” it to process it in L1 caches sized slices, the SoE will perform as predicted whereas the overhead of the SoA continues to rise due to the extra ratio of middle/inner loops for the SoA as compared to the SoE. While the SoE has an increase in culling operations by this about 20%, the SoA has an increase in quadratic processing loops of about this same 20% for a very small net gain if any. If these optimizations were made to the SoA, then it might just catch up to the maximally wheel factorized SoE at some very large sieving range such as ten to the nineteenth power or so, but we are talking hundreds of years of processing on a desktop computer even using multi-processing.

The SoA might be practical for sieving a sub-range (ie. a single page segment) from a very large range offset (say ten raised to the twenty-second power), but there is still that cost of extra memory or computation in dealing with all those values of ‘x’ and ‘y’ that are not necessary with the SoE. END_EDIT_ADD2

Just as for my first answer, the comparison says that SoA loses to the highly wheel factorized SoE on almost all fronts, as follows:

  1. The SoA code is much more complex and it is harder to write and maintain, thus much harder to implement page segmentation which is so important for larger ranges in order to take advantage of maximum CPU capabilities as to memory access speeds.

  2. The SoA likely loses on CPU clock cycles by at least 20% due to slightly more complex inner loop code and higher numbers of outer loops.

  3. The SoA has slightly more operations as compared to a highly wheel factorized SoE (higher wheel factorization than as per the above discussion) for a reasonably large range of one to four billion.

  4. The SoA wins for very large ranges by having an improvement factor of log (log n) as ranges get bigger, but that factor would only allow the SoA to catch up to the highly wheel factorized SoE at a very high number range (perhaps at about 10 raised to the nineteenth power) if ever, and then only if the relative complexity stays the same, which it isn’t likely to do.

So I say again: ” Why use the Sieve of Atkin when the Sieve of Eratosthenes is better in almost every respect and especially is less complex while having less of a constant computational factor overhead? ” Based on my conclusions here that I feel the SoA is an inferior sieve for any kind of reasonable sieve range and will likely never write a page segmented version of the Sieve of Atkin but have written many of them in various computer languages for the Sieve of Eratosthenes. While the Sieve of Atkin is interesting intellectually and mathematically, it isn’t really practical and the Atkin and Bernstein comparison was flawed in overly restricting the wheel factorization of their implementation of the Sieve of Eratosthenes used in their comparison to show an advantage of about 30%, which should have actually been about 60% other than for the constant implementation overheads of the Sieve of Atkin.

    Intereting Posts