Números primos

El fin estaba algo aburrido y venía de leer el siguiente artículo:

http://tech-queries.blogspot.com/2009/04/amazon-interview.html

Lo conseguí ya que me dio curiosidad como es el Fibonacci con eficiencia (n*log(n)) (http://tech-queries.blogspot.com/2010/09/nth-fibbonacci-number-in-ologn.html)

En eso me dio por hacer un problema en http://www.spoj.pl. Entre en mi cuenta, en la cual tengo como 7 problema hechos y vi cual hacer. Elegí uno de generar números primos( http://www.spoj.pl/problems/PRIME1/). Se da una serie de intervalos y hay que imprimir los números primos que hay en ellos.

Empece a programar en Python y yo ya conocía de la criba de Eratóstenes (http://es.wikipedia.org/wiki/Criba_de_Erat%C3%B3stenes) para generar números primos. Fui por una solución simple, buscaba el máximo número en el cual un intervalo terminaba y calculaba los primos hasta ahí y luego imprimía los que había en cada uno de los intervalos.

#!/usr/bin/python
# -*- coding: utf-8 -*-

from sys import stdin

def sieve(end):
    primes = range(2, end+1)
    for i in primes:
        j = 2
        while i * j <= primes[-1]:
            if i * j in primes:
	            primes.remove(i*j)
            j += 1
    return primes


if __name__ == "__main__":
    end = 0
    ncases = int(stdin.readline())
    cases = []

    for ncase in xrange(ncases):
        case = stdin.readline()
        case_start, case_end = map(int, case.split(' '))
        cases.append({'start' : case_start, 'end' : case_end})
        end = max(end, case_end)

    primes = sieve(end)

    count = 1

    for case in cases:
        for prime in primes:
            if prime < case['start']:
                continue
            elif prime <= case['end']:
                print prime
            else:
                break

        if count != ncases:
            print ''

        count += 1

No cumplió por problemas de memoria, ya que el arreglo de primos era muy grande, así que decido hacer lo mismo, pero en C++. Me encargue de hacerlo super eficiente en memoria, usando 1 bit por cada primo, en vez de 1 byte. Con mi computadora me tomó unos minutos sacar todos los primos hasta 1000000000. Me generó un archivo de 500 megas.

#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>
#include <stdlib.h>
using namespace std;

const char masks [] = {-2, -3, -5, -9, -17, -33, -65, 127};

char *sieve(int max)
{
    int max_index = max/8 + 1;
    char *primes = new char[max_index];

    for(int i = 0; i < max_index; i++)
        primes [i] = 255;

    for(int i = 2; i <= max; i++)
    {
        int j = 2;
        int mult =  0;
        while((mult =  i * j) <= max)
        {
            int index = (mult - 1)/8;
            primes[index] &= masks[(mult - 1) % 8];
            j++;
        }

    }

    return primes;

}

int main()
{

    int ntest;
    int case_start, case_end;
    int max_end = 0, test = 1;
    vector<pair<int, int> > cases;

    cin >> ntest;

    for(; test <= ntest; test++)
    {
        cin >> case_start;
        cin >> case_end;

        pair<int, int> test_case(case_start, case_end);

        cases.push_back(test_case);

        max_end = (case_end > max_end ? case_end : max_end);

    }

    char *primes = sieve(max_end);

    vector<pair<int, int> >::iterator it;
    test = 1;

    for ( it=cases.begin() ; it != cases.end(); it++, test++ )
    {
        pair<int, int> test_case = *it;
        case_start = test_case.first;
        case_end = test_case.second;

        if(!(case_start & 1))
            case_start++;

        for(case_start; case_start <= case_end; case_start += 2)
        {
            if(case_start == 1) continue;
            int index = (case_start -1)/8;
            if(primes[index] & (1 << ((case_start -1) % 8)))
                cout << case_start << endl;
        }

        if(test != ntest)
            cout << endl;
    }

}

Esta vez se cayó por tiempo. Era claro que no use el modo correcto para calcular las cosas. Pensé un rato y leí un poco por Internet. Resulta que existe una versión segmentada de la criba de Eratóstenes.  Se calculan los primos hasta la raíz cuadrada del número (en este caso 32000) y luego con esos primos, se revisa en el intervalo si algo de esos es múltiplo de él, sino es un primo. Lo implemente y por tratar de programar de modo legible, entendible y bonito, se me cayo por tiempo. A la final elimine funciones e hice casi todo en la función principal y lo logre. Perdí un buen tiempo, pero fue interesante.

#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>
#include <stdlib.h>
#include <string.h>
using namespace std;

#define MAX_FIRST_SIEVE_INDEX 32000
#define MAX_FIRST_SIEVE_PRIME 32000
#define MAX_PRIMES 4000

bool sieve[MAX_FIRST_SIEVE_INDEX];
int primes[MAX_PRIMES];

void gen_first()
{
    int i = 0, k = 0;

    memset(sieve, true, MAX_FIRST_SIEVE_INDEX);

    for(i = 2; i <= MAX_FIRST_SIEVE_PRIME; i++)
    {
        int j = 2;
        int mult =  0;
        if(sieve[i-1]) 
            primes[k++] = i;

        while((mult =  i * j) <= MAX_FIRST_SIEVE_PRIME)
        {
            sieve[mult-1] = false;
            j++;
        }

    }

}

int main()
{
    int ntest;

    cin >> ntest;
    gen_first();

    for(int test = 1; test <= ntest; test++)
    {
        int case_start, case_end;

        cin >> case_start;
        cin >> case_end;

        for(int i = (case_start > 2 ? case_start : 2); i <= case_end; i++)
        {
            bool print = true;

            for(int j = 0; primes[j]*primes[j] <= i; j++)
            {
                if(i % primes[j] == 0 )
                {
                    print = false;
                    break;
                }
            }
            if(print) cout << i << endl;

        }

        if(test != ntest)
            cout << endl;
    }


}