Calculando números primos

De vez em quando alguns algoritmos ficam pipocando na minha cabeça. Alguns úteis, outros inúteis. E como estava a muito tempo sem programar, resolvi fazer um programinha em C para calcular números primos.

Abaixo vai o programa, e algumas das observações sobre o algoritmo empregado.

Atenção, este artigo tem um nível alto de "nerdice", portanto leia por sua conta e risco.

Existem várias observações, e otimizações, sobre como se descobrir números primos. A primeira é que o único número par primo é o 2. Todos os outros números pares não são primos, pois são divisíveis por 2. Outra observação é que todos os números com o último dígito 5 não são primos, exceto o próprio 5. Se pudermos eliminar estes números dos testes ganhamos agilidade. Em cada 10 números, passamos a testar somente 4, reduzindo a quantidade de testes em 60%. Mas tem um número primo entre 2 e 5, que é o 3. Se pudermos retirar os divisíveis por 3, em cada 30 números teremos 8 testes de primos, 73.333% de testes a menos.

Como o ciclo de divisíveis por 2, 3 e 5 se repete de 30 em 30 (Mínimo Múltiplo Comum, MMC = 2*3*5), então basta fazer um loop de 30 em 30 e testar os números dentro de cada grupo de 30 que não são divisíveis por estes e números primos. Poderia se acrescentar o 7? Sim, mas o grupo passaria a ser de 210 (MMC = 2*3*5*7) ao invés de só 30, complicando o código para um ganho que seria pequeno.

Como se testa se um número é primo? Basta dividir pelos primos anteriores e ver se tem resto zero em alguma(s) das divisões. Se tiver resto zero em qualquer divisão por um primo, o número é divisível por aquele primo, portanto não é primo.

Outro truque está em como parar o teste. O normal é testar com todos os números primos, mas isto testará demais. Não adianta testar dividindo por um primo maior que a metade dele, não? Eu digo o mesmo quanto ao terço e ao quinto. Vejamos a seguinte seguinte formula: a=b*c. A variável a é igual ao produto de b por c. Se a é constante, então se b aumentar c tem que diminuir, e vice-versa. Existe uma situação na qual b e c são iguais, e esta será quando b e c forem a raiz quadrada de a. Isto faz pensar: "Vale a pena testar os primos conhecidos maiores que a raiz quadrada do número que se testa?". Não, pois os outros componentes, no exemplo a outra variável do produto, será menor.

Um exemplo é 28, que é divisível por 2 e 7 (Por 4 também, mas este não é primo.). A raiz quadrada de 28 é um pouco maior que 5 (5.2915026221292 aproximadamente), e os primos que seriam testados pelo critério acima são 2 3 e 5, mas 2 já mostraria que este número não é primo. Com 14, a raiz quadrada é aproximadamente 3.741657386774, e os primos testados seriam 2 e 3, e o 2 já mostraria que o número não é primo, não precisando testar se é divisível por 7. Se um número tem um primo maior que a raiz quadrada dele, o outro fator do produto será um primo menor que a raiz quadrada dele.

Re-explicando. Um número não primo pode ser decomposto como um produto de números primos (28=2*2*7, 15=3*5, 33=3*11, 27=3*3*3 etc). Se tem um primo que o divide maior que a raiz quadrada dele, tem um outro primo, que o divide, menor que a raiz quadrada dele.

Tem que se testar inclusive a raiz quadrada, como é o caso do 49, que é 7 vezes 7.

Na realidade, a maioria dos números não primos deve ser detectado nos primeiros primos testados, bem antes de chegar na raiz quadrada dele. Então a condição de parada é mais importante para números primos, para não fazer testes desnecessários.

O que é mais rápido, calcular o quadrado ou a raiz quadrada? A raiz quadrada de um número qualquer, mesmo que seja inteiro, é um numero real. Isto implica em possíveis arredondamentos e comparações mais complexas. Enquanto o quadrado de um número inteiro é um número inteiro, e pode ser obtido somente com a multiplicação dele por ele mesmo. Desde o tempo do 486, se não antes, a multiplicação de dois números inteiros é muito rápida.

Já que o loop é de 30 em 30, pode-se colocar os primeiros primos "hardcoded" e podemos começar a calcular mesmo à partir do 30. Poderia-se começar do 10, mas resolvi fazer ficar bonitinho. E ainda, como o loop elimina os múltiplos de 2, 3 e 5, então começo o teste de primos pulando estes valores. De que adianta testar o que não vai acontecer mesmo?

Como já temos a explicação básica do algoritmo, eis o programa:

#include        <stdio.h>

#define MAX_PRIMOS      (128000000)

typedef unsigned long long      t_primos        ;

static t_primos primos[ MAX_PRIMOS ]    ;


static  int     testa_primo( register t_primos primo )
{
        /* O loop foi feito para já ignorar os divisíveis por 2,3 e 5.  */
        register t_primos       *p_primos = primos + 3 ;


        for( ; primo >= (*p_primos)*(*p_primos) ; p_primos++ )
                if( ( primo % *p_primos ) == 0)
                        return 0 ;

        return 1 ;
}

main()
{
        register t_primos       primo_base      ,
                                *p_primo        ;


        /* Inicialmente alguns valores pré-definidos.   */
        primos[ 0 ] = 2 ;
        primos[ 1 ] = 3 ;
        primos[ 2 ] = 5 ;
        primos[ 3 ] = 7 ;
        primos[ 4 ] = 11 ;
        primos[ 5 ] = 13 ;
        primos[ 6 ] = 17 ;
        primos[ 7 ] = 19 ;
        primos[ 8 ] = 23 ;
        primos[ 9 ] = 29 ;
        p_primo = primos + 10 ;

        for( primo_base = 30 ; primo_base < MAX_PRIMOS*4 ; primo_base +=30 )
        {
                if( testa_primo( primo_base + 1 ) )
                {
                        *(p_primo++) = primo_base + 1 ;

                        printf( "%d\n",primo_base + 1 );
                }

                if( testa_primo( primo_base + 7 ) )
                {
                        *(p_primo++) = primo_base + 7 ;

                        printf( "%d\n",primo_base + 7 );
                }

                if( testa_primo( primo_base + 11 ) )
                {
                        *(p_primo++) = primo_base + 11 ;

                        printf( "%d\n",primo_base + 11 );
                }

                if( testa_primo( primo_base + 13 ) )
                {
                        *(p_primo++) = primo_base + 13 ;

                        printf( "%d\n",primo_base + 13 );
                }

                if( testa_primo( primo_base + 17 ) )
                {
                        *(p_primo++) = primo_base + 17 ;

                        printf( "%d\n",primo_base + 17 );
                }

                if( testa_primo( primo_base + 19 ) )
                {
                        *(p_primo++) = primo_base + 19 ;

                        printf( "%d\n",primo_base + 19 );
                }

                if( testa_primo( primo_base + 23 ) )
                {
                        *(p_primo++) = primo_base + 23 ;

                        printf( "%d\n",primo_base + 23 );
                }

                if( testa_primo( primo_base + 29 ) )
                {
                        *(p_primo++) = primo_base + 29 ;

                        printf( "%d\n",primo_base + 29 );
                }
        }
}

Talvez possa fazer alguma otimização a mais no programa, mas não me veio muitas à cabeça ainda.

Fiz algumas variações e testei os tempos, e foram equivalentes. As diferenças de tempo de execução foram irrisórias e aleatórias, sendo sempre na ordem de 3 por 10 mil, 0.03%, i.e., dentro dos erros de medida. O compilador GCC deve ter feito muita otimização do código de modo que fossem equivalentes os 3 fragmentos abaixo:

"
                if( testa_primo( primo_base + 29 ) )
                {
                        *(p_primo++) = primo_base + 29 ;

                        printf( "%d\n",primo_base + 29 );
                }
"

"
                if( testa_primo( tmp = primo_base + 29 ) )
                {
                        *(p_primo++) = tmp ;

                        printf( "%d\n",tmp );
                }
"

"
                if( testa_primo( primo_base + 29 ) )
                        printf( "%d\n",*(p_primo++) = primo_base + 29 );
"

Outro teste que fiz foi colocar a função testaprimo() como inline, ao invés  de static, para que o código seja inserido a cada chamada. A diferença foi irrisória, mas foi para melhor, porém tudo dentro da margem de erro da medida, o que não dá garantias de que realmente foi vantajoso. O fato de colocar a função em static já dizia que ela está restrita à este arquivo, e não seria chamada em mais nenhum outro lugar, além deste arquivo fonte. Isto dá uma dica ao compilador para já fazer um inline dela, ou qualquer outra otimização na chamada que ele achar possível. Só não tenho certeza se ele o fez por que o tamanho do arquivo objeto, que não variou com as 3 variações acima, cresceu com o inline, portanto deve ter gerado um código diferente.

Estou colocando este programa sob licença GLP 2. E não usem para enganar o professor em trabalho escolar de programação.

Nota ao professor: Se um aluno apresentou este programa como trabalho em matéria, pode dar ZERO. Ele merece. Se ele leu o texto se fez o programa dele, merce uma boa nota, mas se apresentou o meu programa como trabalho, dê zero a ele.