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;
    }


}
Advertisements

Mi tesis

Decidí escribir sobre mi tesis.

Primero hablaré el porque tome este camino y no pasantía. Lo primero, trabajo para mi mismo y mi compañero, es decir, yo dispongo de mi tiempo. Segundo, trabajar usualemente es algo muy mecánico y en la cual uno no toca temas interesantes  de la teoría de la computación, la cual me parece full interesante. Con éstas razones, sumándole de tener un buen compañero y una profesora dispuesta a ser tutora, se dió el proyecto.

Mi tesis trata sobre clustering de datos, el cual hoy en día es usado en diversas aplicaciones, como lo son Data Mining, preprocesamiento de imágenes, etc. Es una búsqueda no supervisada (no se tiene información previa y se clasifica desde cero). Es un problema NP-Hard.

La manera para atacarlo fue mediante metaheurísticas. Éstas son algoritmos, génericos, por ello la palabra meta, que usan suposiciones/intuiciones que no siempre dan resultados exactos, de ahí la palabra heurística. Hay montones: genético, búsqueda Tabu, algoritmo hormiga, PSO, simulated anieling, etc. Mi compañero y yo implementamos 5 de ellas, basadas todas en población e inteligencia colectiva.

La implementación fue en C++. Para este tipo de problemas se quiere un lenguaje que compile a código de máquina. Se busca eficiencia. La diferencia entre C y C++ es mínima en tiempo, siendo este último preferible por ser orientado a objetos.

De todo lo que se hizo, lo que mejor quedó fue el algoritmo genético. Dió los mejores resultados, es muy flexible, eficiente, etc. Mi compañero y yo hemos estados haciendo pruebas, cambiando y mejorando el código para ver si con la tutora, que tenemos abandonada, logramos escribir una publicación. El código ha sido reescrito y mejorado. Ahora es super entendible y a su vez bastante efieciente. Ya en uno de mis post toque el tema de mis intentos y logros por hacer el código mas efieciente.

Es bastante divertido saber la cantidad de lenguajes que usamos. Además de C++, se usan Python, Haskell y Latex, cada uno para cosas diferentes. Para el manejo de versiones se usa GIT. Ya veremos que sucede, posteo unas imágenes para demostrar lo que es capaz de hacer:

Antes:

Peppers

Después:

Peppers, Resultado

Antes:

Lenna

Después:

Lenna, Resultado

Antes:

Fotografía

Después:

Fotografía, Resultado

Finalmente una imagen de la cual se pueden conseguir patrones escondidos:

John Lennon

Eso es todo, hablamos !

Tratando de optimizar el código

Ya casi no tengo mucho tiempo para escribir cosas interesantes :S!

En los siguientes post me gustaria hablar de javascript y Google Maps.

Mi tesis(https://github.com/fedep3/Metaheuriscticas-Data-Clustering) está hecha en C++ y es sobre data clustering  de datos numéricos usando metaheurísticas. Esta enfocada en las imagenes principalmente (las imágenes no son mas que conjuntos de vectores de uno o mas colores).

Actualmente ando investigando y mejorando el código para hacer una poblicacion sobre el algoritmo genético que cree yo con mi companero.

En el código he hecho dos cambios:

  1. Cambiar el generador de numeros aleatorios.
  2. Intentar optimizar la función de la distancia euclideana.

El motivo de l primer cambio era que el generador de numeros aleatorios de C++, es un http://en.wikipedia.org/wiki/Linear_congruential_generator , los cuales son muy eficiente, pero no son tan buenos. Habiendo usado python y ya de hace un tiempo que habia leído del tema decidí colocar el http://en.wikipedia.org/wiki/Mersenne_twister que es de los mejores, por no decir el mejor actualmente. Ahora el reto era buscar que implementecion usar o si hacerla uno mismo. Buscando y probando consegui uno creado por un invetigador de Pixar que hace uso de las instrucciones SSE2 y que es increiblemente eficiente(http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/cpp_sse2.zip).

Al principio no me compilo asi que busqué lo que estaba pasando. Después de como 1 hora de leer y probar di con la solucion. El código hace uso de unos defines o funciones que por lo visto no estan incluidos en el emmintrin.h. Basto con cambiar el código y listo compilo.  Lo arregle para usarlo y la diferencia de tiempo/resultados con lo que tenia antes ni se notaba.

Averiguando todo este, me surgió un interés por esas instrucciones SSE2. Quería ver si podía programar en ellas la distancia euclideana. Encontre una especia de API de ellas (http://software.intel.com/sites/products/documentation/studio/composer/en-us/2011/compiler_c/index.htm#intref_cls/common/intref_sse2_int_load.htm) y leyendo otras cosas(http://fastcpp.blogspot.comhttps://developer.apple.com/hardwaredrivers/ve/sse.html) entendí su funcionamiento. Las intrucciones son capaces(segun leí y entendí) de procesar 4 flotantes a la vez o 2 doubles.

Mi tesis usa floats, asi que programé lo siguiente:

float Metaheuristic::euclideanDistance(float* v1, float* v2)
{
    __m128 sum = _mm_setzero_ps(), square;
    for(register int k = 0; k < M; k+=4){
        square = _mm_sub_ps( _mm_load_ps(v1 + k), _mm_load_ps(v2 +k));
        sum = _mm_add_ps(_mm_mul_ps(square,square), sum);
    }
    sum = _mm_sqrt_ps(sum);
    _mm_store_ps(eDistance, sum);
    return eDistance[0];
}

Surgieron un montón de problemas para lograr que éesto funcionara. El principal fue que los arreglos que usan estas instrucciones deben tener tamaño divisible entre 4, ya que se cargan los 4 flotantes a los registros. Cambie gran parte del codigo para que a la hora de crear los objetos se cumpliera esta condición y probé de no tener leaks o errors en la memoria usando valgrind. Cuando ya no tuve problemas y los resultados de las corridas tenían coherencia, observé los tiempo. Una imagen que me corria antes en 1 – 2 segundos pasó a 10 segundos !

Mi intento por optimizar el codigo fue fallido, pero aprendi algo. Para mi la razón de que no funcionara muy bien es que los arreglo en realidad eran arreglos de tamano 1 o 3 transformados a 4. Más tiempo se perdía en la carga y cálculo que el que ganaba de tener la capacidad de procesar 4 cosas a la vez. Eso que hasta intente tener las cosas aleanadas a 16 bytes y nada ! Por supuesto que es posible que halla hecho algo mal, pero ahora no tengo mucho tiempo para seguir jugando con el código, ando full de otras cosas.

Espero que ésto sea de ayuda