Problema nº2 de Janeiro de 2005: Goldbach's Conjecture (II)

solução

O segundo problema para esta época em que estamos cheios de tempo para isto (coff.. exames.. coff) é para praticar os números primos. Fácil de perceber, mas não tão fácil de resolver :) Pensem sempre em inputs no limite e repetidos, o programa tem de produzir as respostas rápido!

Goldbach's Conjecture: For any even number n greater than or equal to 4, there exists at least one pair of prime numbers p1 and p2 such that n = p1 + p2.


This conjecture has not been proved nor refused yet. No one is sure whether this conjecture actually holds. However, one can find such a pair of prime numbers, if any, for a given even number. The problem here is to write a program that reports the number of all the pairs of prime numbers satisfying the condition in the conjecture for a given even number.


A sequence of even numbers is given as input. Corresponding to each number, the program should output the number of pairs mentioned above. Notice that we are interested in the number of essentially different pairs and therefore you should not count (p1, p2) and (p2, p1) separately as two different pairs.

Input 

An integer is given in each input line. You may assume that each integer is even, and is greater than or equal to 4 and less than 215. The end of the input is indicated by a number 0.

Output 

Each output line should contain an integer number. No other characters should appear in the output.

Sample Input 

6
10
12
0

Sample Output 

1
2
1

Solução

Vamos pre-calcular os valores para todos os inputs possíveis, são só mais ou menos 2^15 = 32768. Para calcular os pares, temos de gerar todos os primos até 32768. A maneira mais simples é, para cada número, tentar dividir com todos os anteriores; se não der com nenhum, é primo:
for (i = 2; i <= MAX; i++) {
	for (j = 2; j < i; j++) {
		if (i % j == 0)
			break;
	}

	if (j == i) {
		// é primo...
	}
}
Mas podemos fazer melhor que isto; o módulo é uma operação muito pesada e isto é muito lento. Podemos reduzir o número de testes para cada número testando só até sqrt(i). Porque é que isso funciona? Seja r = sqrt(i). Então, r*r = i. Um produto a*b (a,b != r) para dar i, tem a < b e b > r (ou ao contrário). Por isso, em vez de testar a divisibilide com o a e o b, basta tentar com um deles. O maior valor que temos de testar é então r = sqrt(i), e não vale a pena ir além dele porque já testámos com o outro valor mais pequeno.

Uma alteração que já acelera mais alguma coisa é testar só com os primos já encontrados. Por exemplo, testar a divisibilidade por 2 e depois por 4 é redundante, porque se fosse divisível por 4, já teria sido também por 2. Por isso, cada vez que encontramos um primo, adicionamo-lo a uma lista; para testar se um número é primo, testamos a divisibilidade só com os números que estão na lista:

num_primes = 0;
for (i = 2; i < MAX; i++) {
	for (j = 0; j < num_primes; j++) {
		if (i % primes[j] == 0)
			break;
	}

	if (j == num_primes) {
		// é primo
		primes[num_primes] = i;
		num_primes++;
	}
}
Podemos ainda usar a paragem ao atingir um valor superior a sqrt(i), e podemos também começar logo com o 2 na lista, num_primes = 1 e começar de i=3 até MAX com incrementos de 2 (metade dos ciclos). Mas há uma maneira ainda melhor!

Como disse, a divisão (ou o módulo) é uma operação muito lenta. Um algoritmo que faça a mesma coisa sem divisões é muito mais rápido. O que usei foi o Sieve; a ideia é muito simples. Temos um array de variáveis booleanas que dizem se um número é primo ou não, e ao inicío todos são primos (excepto o 0 e o 1). Em cada ciclo, verificamos se o número em que vamos (i) é primo; se não for, avançamos para o próximo i. Se for primo, então adicionamo-lo à lista (se quisermos ter a lista, pode não ser necessário) e eliminamos todos os múltiplos desse primo. Por exemplo, o primeiro valor de i é 2. Então, adicionamos o 2 à lista e corremos outro ciclo com j, de 2+2 até ao limite e de 2 em 2, e para cada número destes marcamos no array como não primo. No caso do 2, seriam o 4, 6, 8, 10, 12, (...). De seguida vem o 3. Como não foi marcado, é primo, e por isso marcamos os múltiplos como não primos: 6 (outra vez, não faz mal), 9, 12, 15, 18, (...). A seguir é o 4, mas já está marcado e passamos para o 5, que não está e é primo; marcamos os múltiplos 10, 15, 20, 25, etc. Simples e rápido (só adições). A desvantagem do Sieve em relação ao outro método é que o Sieve tem de parar de eliminar números num certo sítio, um limite; para este problema é 32768 que cabe bem num array. Noutros problemas o limite é demasiado grande para usar o Sieve. Outro problema é que depois de feitos os cálculos, não podemos aumentar o intervalo (teríamos de repetir as eliminações para cada primo...), enquanto que no outro é só voltar a testar com os primos anteriores.

Bem, cada problema tem a sua solução e esta serve para este. Depois de gerar os primos, percorremos a lista com uma variável e com outra, fazemos a soma, e o contador do resultado é aumentado em 1, e problema resolvido :)

Código em C:

#include <stdio.h>

#define MAX	32769

int ans[MAX];
int isprime[MAX];
int primes[MAX];
int num_primes;

int main(int argc, char *argv[])
{
	int i, j;
	int tmp, n;
	
	num_primes = 0;
	for (i = 2; i < MAX; i++) {
		if (isprime[i] == 0) {
			for (j = i+i; j < MAX; j+=i) {
				isprime[j] = 1;
			}

			primes[num_primes] = i;
			num_primes++;
		}
	}
	
	for (i = 0; i < num_primes; i++) {
		for (j = i; j < num_primes; j++) {
			tmp = primes[i] + primes[j];
			if (tmp < MAX) {
				ans[tmp]++;
			} else {
				break;
			}
		}
	}

	for (;;) {
		scanf("%d\n", &n);
		if (n == 0)
			break;
		printf("%d\n", ans[n]);
	}

	return 0;
}

O valor de MAX é 32769 porque o máximo é 32768, e como os valores dos arrays vão de 0 ao tamanho-1, uso um tamanho um bocado maior. Como os arrays são globais (variáveis estáticas), são inicializados a 0, o que é muito útil; por isso, vamos fazer a batota de usar o 0 para querer dizer "é primo" e o 1 para dizer "não é primo". Isto é o contrário do que costuma ser usado.

Portanto, no primeiro ciclo testamos se o número ainda não foi marcado (se é 0, ou seja, é primo). Então corremos outro ciclo para marcar com 1 (não primo) os múltiplos, e adicionamos este número à lista.

Agora para resolver o problema, para cada primo, executamos outro ciclo que vai desde esse primo até ao fim, soma-os, e aumenta o contador do resultado. Esse contador é iniciado a 0 automaticamente porque, mais uma vez, o array usado é estático e por isso, iniciado a 0. Repare-se que o 2º ciclo (do j) começa desde i e não do início; isto porque 3+5 e 5+3 é a mesma coisa e não deve ser contado duas vezes. Começa desde i porque 5+5 é válido. Assim que se excede o limite, podemos acabar o ciclo do j.

Há outras maneiras de acelerar o Sieve, como fazer fora do ciclo principal a marcação dos múltiplos de 2 e assim o ciclo pode começar em 3 e andar de 2 em 2, passa para metade dos ciclos no total. Também podemos "sievar" só até sqrt(MAX), isto garante que todos os número não primos até MAX já foram marcados. Por exemplo, se um problema pede para testar se um número qualquer até 1.000.000 é primo, basta sievar até 1000. Podemos também sacar a lista dos primos se percorrermos o isprime[].