Pruebas condicionales en primalidad por división de prueba

1 minuto de lectura

Pruebas condicionales en primalidad por division de prueba
bosón z

Mi pregunta es sobre la prueba condicional en la división de juicio. Parece haber cierto debate sobre qué prueba condicional emplear. Veamos el código para esto de Código Rosetta.

int is_prime(unsigned int n)
{
    unsigned int p;
    if (!(n & 1) || n < 2 ) return n == 2;

    /* comparing p*p <= n can overflow */
    for (p = 3; p <= n/p; p += 2)
        if (!(n % p)) return 0;
    return 1;
}

La factorización de la rueda o el uso de una lista predeterminada de números primos no cambiará la esencia de mi pregunta.

Hay tres casos en los que puedo pensar para hacer la prueba condicional:

  1. p<=n/p
  2. p*p<=n
  3. int corte = sqrt(n); para (p = 3; p <= corte; p += 2)

Caso 1: funciona para todos los n pero tiene que hacer una división extra cada iteración (editar: en realidad no necesita una división adicional, pero aún es más lento. No estoy seguro de por qué. Vea el resultado del ensamblaje a continuación). Descubrí que es el doble de lento que el caso 2 para valores grandes de n que son primos (en mi sistema Sandy Bridge).

Caso 2: es significativamente más rápido que el caso 1 pero tiene el problema de que se desborda para n grande y entra en un bucle infinitivo. El valor máximo que puede manejar es

(sqrt(n) + c)^2 = INT_MAX  //solve
n = INT_MAX -2*c*sqrt(INT_MAX) + c^2
//INT_MAX = 2^32 -> n = 2^32 - c*s^17 + c^2; in our case c = 2

Por ejemplo, para uint64_t, el caso 2 entrará en un bucle infinito para x =-1L-58 (x^64-59), que es un número primo.

Caso 3: no es necesario realizar divisiones ni multiplicaciones en cada iteración y no se desborda como en el caso 2. También es un poco más rápido que el caso 2. La única pregunta es si sqrt(n) es lo suficientemente preciso.

¿Puede alguien explicarme por qué el caso 2 es mucho más rápido que el caso 1? El caso 1 NO usa una división adicional como pensé, pero a pesar de eso, sigue siendo mucho más lento.

Estos son los tiempos para el primo 2^56-5;

case 1 9.0s
case 2 4.6s
case 3 4.5s

Aquí está el código que usé para probar esto http://coliru.stacked-crooked.com/a/69497863a97d8953. También agregué las funciones al final de esta pregunta.

Aquí está el formulario de salida de ensamblaje GCC 4.8 con -O3 para el caso 1 y el caso 2. Ambos solo tienen una división. El caso 2 también tiene una multiplicación, por lo que mi primera suposición es que el caso 2 sería más lento, pero es aproximadamente el doble de rápido tanto en GCC como en MSVC. no sé por qué

Caso 1:

.L5:
  testl %edx, %edx
  je  .L8
.L4:
  addl  $2, %ecx
  xorl  %edx, %edx
  movl  %edi, %eax
  divl  %ecx
  cmpl  %ecx, %eax
  jae .L5

Caso 2:

.L20:
  xorl  %edx, %edx
  movl  %edi, %eax
  divl  %ecx
  testl %edx, %edx
  je  .L23
.L19:
  addl  $2, %ecx
  movl  %ecx, %eax
  imull %ecx, %eax
  cmpl  %eax, %edi
  jae .L20

Aquí están las funciones que estoy usando para probar el tiempo:

int is_prime(uint64_t n)
{
    uint64_t p;
    if (!(n & 1) || n < 2 ) return n == 2;

    /* comparing p*p <= n can overflow */
    for (p = 3; p <= n/p; p += 2)
        if (!(n % p)) return 0;
    return 1;
}

int is_prime2(uint64_t n)
{
    uint64_t p;
    if (!(n & 1) || n < 2 ) return n == 2;

    /* comparing p*p <= n can overflow */
    for (p = 3; p*p <= n; p += 2)
        if (!(n % p)) return 0;
    return 1;
}

int is_prime3(uint64_t n)
{
    uint64_t p;
    if (!(n & 1) || n < 2 ) return n == 2;

    /* comparing p*p <= n can overflow */
    uint32_t cut = sqrt(n);
    for (p = 3; p <= cut; p += 2)
        if (!(n % p)) return 0;
    return 1;
}

Contenido agregado después de la recompensa.

Aean descubrió que en el caso 1, guardar el cociente y el resto es tan rápido (o un poco más rápido) que el caso 2. Llamemos a este caso 4. El siguiente código es el doble de rápido que el caso 1.

int is_prime4(uint64_t n)
{
    uint64_t p, q, r;
    if (!(n & 1) || n < 2 ) return n == 2;

    for (p = 3, q=n/p, r=n%p; p <= q; p += 2, q = n/p, r=n%p)
        if (!r) return 0;
    return 1;
}

No estoy seguro de por qué esto ayuda. En cualquier caso, ya no es necesario utilizar el caso 2. Para el caso 3, la mayoría de las versiones del sqrt La función en hardware o software obtiene los cuadrados perfectos para que sea seguro de usar en general. El caso 3 es el único caso que funcionará con OpenMP.

  • El caso 1 en realidad no debería implicar una división adicional: su procesador probablemente sea capaz de calcular tanto el cociente como el resto en una sola instrucción, y eso está expuesto a C (al menos en C99) en forma de div función. ¿Has probado a usar eso?

    –Mark Dickinson

    21 de marzo de 2014 a las 10:58

  • Referencia para div y ldiv: ver apartado 7.20.6.2 de open-std.org/jtc1/sc22/wg14/www/docs/n1256.pdf

    –Mark Dickinson

    21 de marzo de 2014 a las 10:59

  • @MarkDickinson, no, no lo he probado. Supuse que el compilador haría esta optimización por mí (MSVC2012 y GCC 4.8+)…

    – bosón Z

    21 de marzo de 2014 a las 11:49

  • Mmm; sí, supongo que es posible, dependiendo de los niveles de optimización, etc. Podría echar un vistazo a la salida del ensamblado (compilar con -S bandera) para comprobar.

    –Mark Dickinson

    21 de marzo de 2014 a las 11:51

  • @MarkDickinson, investigué esas funciones. Llaman a funciones externas y no alinean una sola instrucción. Sin embargo, observé el resultado del ensamblaje de los casos 1 y 2 y ambos solo dividen una vez. No sé por qué el caso 1 es mucho más lento. Actualicé mi respuesta con la salida del ensamblado.

    – bosón Z

    21 de marzo de 2014 a las 12:50

Pruebas condicionales en primalidad por division de prueba
Aean

UPD: Este es un problema de optimización del compilador, obviamente. Mientras que MinGW usó solo una div instrucción en el cuerpo del bucle, tanto GCC en Linux como MSVC no pudieron reutilizar el cociente de la iteración anterior.

Creo que lo mejor que podemos hacer es definir explícitamente quo y rem y calcularlos en el mismo bloque de instrucciones básico, para mostrarle al compilador que queremos tanto el cociente como el resto.

int is_prime(uint64_t n)
{
    uint64_t p = 3, quo, rem;
    if (!(n & 1) || n < 2) return n == 2;

    quo = n / p;
    for (; p <= quo; p += 2){
        quo = n / p; rem = n % p;
        if (!(rem)) return 0;
    }
    return 1;
}

Probé tu código de http://coliru.stacked-crooked.com/a/69497863a97d8953 en un compilador MinGW-w64, case 1 es más rápido que case 2.

ingrese la descripción de la imagen aquí

Asique adivinar está compilando para una arquitectura de 32 bits y utiliza uint64_t escribe. Su ensamblaje muestra que no usa ningún registro de 64 bits.

Si lo hice bien, ahí está la razón.

En la arquitectura de 32 bits, los números de 64 bits se representan en dos registros de 32 bits, su compilador realizará todos los trabajos de concatenación. Es simple hacer sumas, restas y multiplicaciones de 64 bits. Pero el módulo y la división se realizan mediante una pequeña llamada de función que se nombra como ___umoddi3 y ___udivdi3 en CCG, aullrem y aulldiv en MSVC.

Así que en realidad necesitas uno ___umoddi3 y uno ___udivdi3 para cada iteración en case 1una ___udivdi3 y una concatenación de multiplicación de 64 bits en case 2. Es por eso case 1 parece el doble de lento que case 2 en tu prueba.

Lo que realmente obtienes case 1:

L5:
    addl    $2, %esi
    adcl    $0, %edi
    movl    %esi, 8(%esp)
    movl    %edi, 12(%esp)
    movl    %ebx, (%esp)
    movl    %ebp, 4(%esp)
    call    ___udivdi3         // A call for div
    cmpl    %edi, %edx
    ja  L6
    jae L21
L6:
    movl    %esi, 8(%esp)
    movl    %edi, 12(%esp)
    movl    %ebx, (%esp)
    movl    %ebp, 4(%esp)
    call    ___umoddi3        // A call for modulo.
    orl %eax, %edx
    jne L5

Lo que realmente obtienes case 2:

L26:
    addl    $2, %esi
    adcl    $0, %edi
    movl    %esi, %eax
    movl    %edi, %ecx
    imull   %esi, %ecx
    mull    %esi
    addl    %ecx, %ecx
    addl    %ecx, %edx
    cmpl    %edx, %ebx
    ja  L27
    jae L41
L27:
    movl    %esi, 8(%esp)
    movl    %edi, 12(%esp)
    movl    %ebp, (%esp)
    movl    %ebx, 4(%esp)
    call    ___umoddi3         // Just one call for modulo
    orl %eax, %edx
    jne L26

MSVC no pudo reutilizar el resultado de div. La optimización se rompe por return. Pruebe estos códigos:

__declspec(noinline) int is_prime_A(unsigned int n)
{
    unsigned int p;
    int ret = -1;
    if (!(n & 1) || n < 2) return n == 2;

    /* comparing p*p <= n can overflow */
    p = 1;
    do {
        p += 2;
        if (p >= n / p) ret = 1; /* Let's return latter outside the loop. */
        if (!(n % p)) ret = 0;
    } while (ret < 0);
    return ret;
}

__declspec(noinline) int is_prime_B(unsigned int n)
{
    unsigned int p;
    if (!(n & 1) || n < 2) return n == 2;

    /* comparing p*p <= n can overflow */
    p = 1;
    do {
        p += 2;
        if (p > n / p) return 1; /* The common routine. */
        if (!(n % p)) return 0;
    } while (1);
}

los is_prime_B será dos veces más lento que is_prime_A en MSVC/ICC para Windows.

  • Interesante. Lo ejecuté en MSVC201 versión en modo de 64 bits y obtuve: caso 1 1,8 s, caso 2 1,2 s, caso 3 1,1 segundos. Entonces, el caso dos sigue siendo aproximadamente el doble de lento con MSVC2012 en modo de 64 bits. Déjame intentarlo con GCC en Linux.

    – bosón Z

    25 de marzo de 2014 a las 8:36

  • @Z boson MSVC no pudo reutilizar el resultado de div, te explico más. Por favor, espere un momento.

    – Aean

    25 de marzo de 2014 a las 8:38


  • No, simplemente inicié Linux y compilé con `g++ -O3 -fopenmp prime.C y obtuve los mismos resultados: case 1 2.0s, case 2 1.0 s, case 3 1.1 s. El caso 1 sigue siendo el doble de lento. Tengo un E5-1620 a 3,60 GHz (núcleos Xeon Sandy Bridge).

    – bosón Z

    25 de marzo de 2014 a las 8:44

  • ¿Has comprobado la salida del ensamblado?

    – Aean

    25 de marzo de 2014 a las 8:48

  • No. Lo comprobaré más tarde. Solo para aclarar. Arranqué dualmente, así que probé MSVC y GCC en la misma máquina, solo con un sistema operativo diferente.

    – bosón Z

    25 de marzo de 2014 a las 8:57

sqrt(n) es lo suficientemente preciso siempre y cuando su sqrt es monótono creciente, obtiene cuadrados perfectos correctos, y cada unsigned int puede representarse exactamente como un double. Los tres son el caso en todas las plataformas que conozco.

Puede sortear estos problemas (si considera que son problemas) implementando una función unsigned int sqrti(unsigned int n) que devuelve el suelo de la raíz cuadrada de un unsigned int utilizando el método de Newton. (¡Este es un ejercicio interesante si nunca lo has hecho antes!)

  • ¿Tiene una fuente sobre sqrt obteniendo los cuadrados perfectos? Sé que sin iniciar sesión se puede representar exactamente como doble, pero uint64_t no se puede representar exactamente como doble, por lo tanto, no veo por qué debería garantizarse que los cuadrados perfectos sean correctos. Buen punto sobre hacer esto yo mismo con el método de Newton. He estado pensando en hacer eso.

    – bosón Z

    21 de marzo de 2014 a las 17:32

  • sqrt en la mayoría de las plataformas principales (¿todas?) está “redondeado correctamente” — en el modo de redondeo al más cercano, siempre está dentro de 1/2 ulp de estar en lo correcto, y en los modos de redondeo dirigido, todavía está dentro de 1 ulp de estar en lo cierto . Eso implica que acierta en los casos exactos (ya que hacerlo mal implicaría una diferencia de al menos 1 ulp). IEEE 754 requiere que sqrt proporcionarse y redondearse correctamente. Sin embargo, IEEE 754 no es el estándar C.

    – tmyklebu

    21 de marzo de 2014 a las 22:38


Una respuesta a solo una pequeña parte de esta publicación.

Corrección del caso 2 para lidiar con el desbordamiento.

#include <limits.h>

int is_prime(unsigned n) {
  unsigned p;
  if (!(n & 1) || n < 2)
    return n == 2;

  #define UINT_MAX_SQRT (UINT_MAX >> (sizeof(unsigned)*CHAR_BIT/2))
  unsigned limit = n;
  if (n >= UINT_MAX_SQRT * UINT_MAX_SQRT)
    limit = UINT_MAX_SQRT * UINT_MAX_SQRT - 1;

  for (p = 3; p * p < limit; p += 2)
    if (!(n % p))
      return 0;

  if (n != limit)
    if (!(n % p))
      return 0;
  return 1;
}

El cálculo del límite falla si ambos sizeof(unsigned) y CHAR_BIT son extraños – una situación rara.

Acerca de su primera pregunta: ¿por qué (2) es más rápido que (1)?
Bueno, esto depende del compilador, tal vez.
Sin embargo, en general se podría esperar que una división sea una operación más costosa que una multiplicación.

Acerca de su segunda pregunta: ¿es sqrt () una función precisa?

En general, es exacto.
El único caso que podría darte problemas es el que sqrt(n) es un número entero.
Por ejemplo, si n == 9 y sqrt(n) == 2.9999999999999 en su sistema, entonces tiene problemas porque la parte entera es 2, pero el valor exacto es 3.
Sin embargo, estos casos raros pueden manejarse fácilmente agregando una constante doble no tan pequeña como 0.1, digamos.
Así, puedes escribir:

  double stop = sqrt(n) + 0.1;  
  for (unsigned int d = 2; d <= stop; d += 2)
       if (n % d == 0)
           break;  /*  not prime!! */

El término agregado 0.1 podría agregar una iteración a su algoritmo, lo cual no es un gran problema en absoluto.

Finalmente, la opción obvia para su algoritmo es (3), es decir, el sqrt() enfoque, porque no hay ningún cálculo (multiplicaciones o divisiones), y el valor stop se calcula una sola vez.

Otra mejora que puedes tener es la siguiente:

  • Tenga en cuenta que todo número primo p >= 5 tiene la forma 6n - 1 o bien 6n + 1.

Entonces, puedes alternar los incrementos de la variable d siendo 2, 4, 2, 4, y así sucesivamente.

¿Ha sido útil esta solución?

Esta web utiliza cookies propias y de terceros para su correcto funcionamiento y para fines analíticos y para mostrarte publicidad relacionada con sus preferencias en base a un perfil elaborado a partir de tus hábitos de navegación. Al hacer clic en el botón Aceptar, acepta el uso de estas tecnologías y el procesamiento de tus datos para estos propósitos. Configurar y más información
Privacidad