elhacker.net cabecera Bienvenido(a), Visitante. Por favor Ingresar o Registrarse
¿Perdiste tu email de activación?.

 

 


Tema destacado: Rompecabezas de Bitcoin, Medio millón USD en premios


+  Foro de elhacker.net
|-+  Programación
| |-+  Programación C/C++ (Moderadores: Eternal Idol, Littlehorse, K-YreX)
| | |-+  Calculo de pi en alta precisión (aporte)
0 Usuarios y 1 Visitante están viendo este tema.
Páginas: [1] 2 3 Ir Abajo Respuesta Imprimir
Autor Tema: Calculo de pi en alta precisión (aporte)  (Leído 11,133 veces)
engel lex
Moderador Global
***
Desconectado Desconectado

Mensajes: 15.514



Ver Perfil
Calculo de pi en alta precisión (aporte)
« en: 9 Abril 2014, 08:16 am »

Buenas, en estos días estaba por aquí un alguien preguntando sobre un programa básico con pi, por cosas de la vida, decidí ahondar en el tema... calcular pi con alta precisión...

deja claro que el siguiente tema no lo hago en pro de discutir la naturaleza de pi, ni nada al respecto, es solo en pro del calculo de alta precisión ;)

aquí básicamente solo hay 2 temas importantes que abordar

-¿qué formula usar?
-¿como calculo números en alta precisión?

con lo primero... la formula, decidí usar la de euler

escogí esta por su precisión y su relativamente fácil aplicación

para aplicarla la tenemos que separar en 2 partes por su sumatoria y su factorial...


Código
  1. sumatoria = factorial = 1;
  2. for(i=1;i<ciclos;i++){
  3.    factorial *= i/(i*2+1);
  4.    sumatoria += factorial;
  5. }
  6. sumatoria *=2;
  7. pi = sumatoria
  8.  

ahí resolvemos esa formula básicamente... pero nos damos cuenta de un problema... los numeros double están limitados a 64 bits... unos 20 decimales para nuestro caso...

en este caso viene la aritmetica de alta precisión, GMP... aquí en el foro buscaba como usarlo... al final conseguí y publiqué el como...

este usa unidades especificas y operaciones especificas para alta precisión...

doy una ligera explicación de lo usado... no voy a caer en mucho detalle de las funciones...
mpf_t es una variable de tipo "float" propia
mpf_set_default_prec da el valor en bytes que la variable usará
mpf_init_set_str inicializa las variables desde un string... por que un string? porque quise -.-... se hay que inicializar las variables obligatoriamente

mpf_set_ui da un valor desde un entero sin signo a un float
div= division, mul = multiplicaccion, add= adición

mpf_get_str convierte de float de GMP a char* (string)

como no habia una multiplicacion / asignacion directa, me tocó usar una variable de intercambio


ahora la parte interesante... el codigo  ;-)

hice un pequeño/bonito formato de impresión de 10 en 10 dígitos... en  lineas de 50 y bloques de 500 con el numero de dígito como indicador... probé solo hasta 40.000 dígitos, comprobando en internet precisos... pero en 5min 9 seg! :P es mi único intento en este campo... no esperen romper el record mundial de decimales ;P por cierto... de ahí en adelante, los dígitos son cuesta arriba... así que cuidado con el tiempo... pueden usar el programa pasando el numero de dígitos como argumento... los dígitos los toma en cuanta de 100 en 100

pd: la precisión la hago como  "precision = digitos * 100 * 3 * 1.12" porque descubrí que se requieren 3 ciclos por dígito correcto sobre los 100 con un 12% de error (calculado a ojo) :P


---------------------------------------------------------------codigo final con formulas a partir de aqui---------------------------------------------------------------

librerias
Código
  1. #include <iostream>
  2. #include <stdlib.h>
  3. #include <string>
  4. #include <gmp.h>
  5. #include <gmpxx.h>
  6. #include <time.h>
  7. #include <sstream>
  8. using namespace std; // :P
  9.  

prototipos
Código
  1. string metodo_euler(unsigned long int digitos, mp_exp_t &exponente);
  2. string metodo_ramanujan(unsigned long int digitos, mp_exp_t &exponente);
  3. void imprimir_pi_aux(string pi, mp_exp_t exponente, int digitos);

main
Código
  1. int main(int argc, char **argv) {
  2.    unsigned long int digitos = strtol(argv[1], NULL, 10); //argv a int
  3.    mpf_set_default_prec(32 * 10 * (digitos / 100) * 1.11); //precision + 10%
  4.    long exponente;
  5.    clock_t start = clock();
  6.    string pi;
  7.    int metodo = 0;
  8.    int i = 1;
  9.    for (i = 0; i < argc; i++) {
  10.        if (argv[i][0] == '-') {
  11.            switch (argv[i][1]) {
  12.                case 'e':
  13.                    metodo = 0;
  14.                    break;
  15.                case 'r':
  16.                    metodo = 1;
  17.                    break;
  18.            }
  19.        }
  20.    }
  21.    switch (metodo) {
  22.        case 0:
  23.            cout << "imprimiendo por metodo Euler, " << digitos * 3 * 1.12 << " ciclos"<<endl;
  24.            pi = metodo_euler(digitos, exponente);
  25.            break;
  26.        case 1:
  27.            cout << "imprimiendo por metodo Ramanujan, " << (int) (digitos / 8 * 1.01) << " ciclos"<<endl;
  28.            pi = metodo_ramanujan(digitos, exponente);
  29.            break;
  30.    }
  31.  
  32.    clock_t end = clock();
  33.  
  34.    imprimir_pi_aux(pi, exponente, digitos);
  35.    cout << "\ntiempo total de ejecucion: " << (float) (end - start) / CLOCKS_PER_SEC << "\n";
  36.    cout << "tiempo total de impresion: " << (float) (clock() - end) / CLOCKS_PER_SEC << "\n\n";
  37.    return 0;
  38. }

tiene 2 argumentos... primero los digitos (obligatorio)... luego -r para ramanujan o -e para euler (opcional, euler predeterminado)

imprimir_pi_aux

para tener un formato único de impresión para todo

Código
  1. void imprimir_pi_aux(string pi, long exponente, int digitos) {
  2.    cout << "Pi: ";
  3.    stringstream aux;
  4.    unsigned long int i;
  5.    if (exponente <= 0) {
  6.        cout << "0.";
  7.        while (exponente < 0) {
  8.            cout << 0;
  9.            exponente++;
  10.        }
  11.    } else {
  12.        for (i = 0; i < exponente; i++) {
  13.            cout << pi[i];
  14.        }
  15.        cout << ".\n\n1:\t";
  16.    }
  17.    for (i = exponente; i < digitos + 1; i++) {
  18.        aux << pi[i];
  19.        if ((i - exponente + 1) % 10 == 0) aux << ' ';
  20.        if ((i - exponente + 1) % 500 == 0) aux << endl;
  21.        if ((i - exponente + 1) % 50 == 0 && i < digitos) aux << endl << i + 1 << ":\t";
  22.    }
  23.    aux << endl;
  24.    cout << aux.str();
  25. }
aux por consejo de amchacon

metodo_euler

corto, simple, pero largo... para 10.000 dígitos 33.600 ciclos o para mi unos 11 seg

basado en la formula



Código
  1. string metodo_euler(unsigned long int digitos, long &exponente) {
  2.    digitos /= 100; //trabaja con lotes de 100
  3.    unsigned long int precision = digitos * 100 * 3 * 1.12;
  4.    unsigned long int x, i;
  5.    mpf_t factorial, sumatoria, buff_factorial; //variables GMP
  6.    mpf_init_set_ui(factorial, 1); //factorial=1
  7.    mpf_init_set_ui(sumatoria, 1); //sumatoria=1
  8.    mpf_init_set_ui(buff_factorial, 1); //buff_factorial=1
  9.  
  10.    for (x = 1; x < precision; x++) {
  11.        i = 2 * x + 1;
  12.        mpf_set_ui(buff_factorial, x); //buff_factorial=x
  13.        mpf_div_ui(buff_factorial, buff_factorial, i); //buff_factorial*=i
  14.        mpf_mul(factorial, factorial, buff_factorial); //factorial*=buff_factorial
  15.        mpf_add(sumatoria, sumatoria, factorial); //sumatoria+=factorial
  16.    }
  17.    mpf_mul_ui(sumatoria, sumatoria, 2); //sumatoria*=2
  18.    return mpf_get_str(NULL, &exponente, 10, 0, sumatoria);
  19. }

metodo_ramanujan

por consejo de do-while, basado en la formula

eficiente, pero dificil de aplicar (si no, miren el codigo de abajo)
para 10.000 digitos 1.262 ciclos en 1.1 segundos



Código
  1. string metodo_ramanujan(unsigned long int digitos, long &exponente) {
  2.    digitos /= 8;
  3.    digitos *= 1.01;
  4.    unsigned long int i, j;
  5.    mpf_t sumatoria, buff_sumatoria, buff_ciclo_superior, buff_ciclo_inferior, primera_parte,
  6.            buff_ciclo_inferior2, factorial_superior, factorial_inferior;
  7.    mpf_init_set_ui(sumatoria, 1103); //sumatoria=0
  8.    mpf_init_set_ui(buff_sumatoria, 0); //buff_sumatoria=0
  9.    mpf_init_set_ui(buff_ciclo_superior, 0); //buff_ciclo_superior=0
  10.    mpf_init_set_ui(buff_ciclo_inferior, 0); //buff_ciclo_inferior=0
  11.    mpf_init_set_ui(primera_parte, 0); //primera_parte=0
  12.    mpf_init_set_ui(buff_ciclo_inferior2, 0); //buff_ciclo_inferior2=0
  13.    mpf_init_set_ui(factorial_superior, 1); //factorial_superior=1
  14.    mpf_init_set_ui(factorial_inferior, 1); //factorial_inferior=1
  15.    mpf_sqrt_ui(primera_parte, 2); //primera_parte=sqrt(2)
  16.    mpf_mul_ui(primera_parte, primera_parte, 2); //primera_parte*=2
  17.    mpf_div_ui(primera_parte, primera_parte, 9801); //primera_parte/=9801
  18.  
  19.    for (i = 1; i <= digitos; i++) {
  20.        for (j = (i - 1)*4 + 1; j <= i * 4; j++) {
  21.            mpf_mul_ui(factorial_superior, factorial_superior, j);
  22.        }
  23.        mpf_set_ui(buff_ciclo_superior, 26390);
  24.        mpf_mul_ui(buff_ciclo_superior, buff_ciclo_superior, i);
  25.        mpf_add_ui(buff_ciclo_superior, buff_ciclo_superior, 1103);
  26.        mpf_mul(buff_ciclo_superior, buff_ciclo_superior, factorial_superior);
  27.        mpf_mul_ui(factorial_inferior, factorial_inferior, i);
  28.        mpf_pow_ui(buff_ciclo_inferior, factorial_inferior, 4);
  29.        mpf_set_ui(buff_ciclo_inferior2, 396);
  30.        mpf_pow_ui(buff_ciclo_inferior2, buff_ciclo_inferior2, 4 * i);
  31.        mpf_mul(buff_ciclo_inferior, buff_ciclo_inferior, buff_ciclo_inferior2);
  32.        mpf_div(buff_sumatoria, buff_ciclo_superior, buff_ciclo_inferior);
  33.        mpf_add(sumatoria, sumatoria, buff_sumatoria);
  34.  
  35.    }
  36.    mpf_mul(sumatoria, sumatoria, primera_parte);
  37.    mpf_ui_div(sumatoria, 1, sumatoria);
  38.    return mpf_get_str(NULL, &exponente, 10, 0, sumatoria);
  39. }


« Última modificación: 11 Abril 2014, 08:02 am por engel lex » En línea

El problema con la sociedad actualmente radica en que todos creen que tienen el derecho de tener una opinión, y que esa opinión sea validada por todos, cuando lo correcto es que todos tengan derecho a una opinión, siempre y cuando esa opinión pueda ser ignorada, cuestionada, e incluso ser sujeta a burla, particularmente cuando no tiene sentido alguno.
1mpuls0


Desconectado Desconectado

Mensajes: 1.186


Ver Perfil
Re: Calculo de pi en alta precisión (aporte)
« Respuesta #1 en: 9 Abril 2014, 17:45 pm »

Interesante  ;-)

Hace algunos años (como por el 2009) tuve la idea realizar un programa para calcular PI con la mayor cantidad de dígitos pero al investigar, métodos, métodos, algoritmos, proyectos me topé con la sorpresa que en ese tiempo se habían calculado 2.5 trillones de dígitos decimales

Creo que el nuevo record está en 10 trillones de dígitos decimales  :xD

:http://www.numberworld.org/digits/Pi/

Pero me parece interesantes algunas deducciones que hiciste en tu implementación.

PD. Por ahora tengo la idea de hacer un colisionador md5  :silbar:


En línea

abc
ivancea96


Desconectado Desconectado

Mensajes: 3.414


ASMático


Ver Perfil WWW
Re: Calculo de pi en alta precisión (aporte)
« Respuesta #2 en: 9 Abril 2014, 18:39 pm »

Está muy bien. Una alternativa al problema de la velocidad, podría ser ir introduciendo cifra a cifra en un fichero. Pero eso ya sale del tema, ya que cambiaría el algoritmo xD
En línea

dato000


Desconectado Desconectado

Mensajes: 3.034



Ver Perfil
Re: Calculo de pi en alta precisión (aporte)
« Respuesta #3 en: 9 Abril 2014, 19:49 pm »

Creo que el nuevo record está en 10 trillones de dígitos decimales  :xD

Got to be Fucking Kidding me  :o :o :o :o :o :o :o

Excelente código
En línea


amchacon


Desconectado Desconectado

Mensajes: 1.211



Ver Perfil
Re: Calculo de pi en alta precisión (aporte)
« Respuesta #4 en: 9 Abril 2014, 20:04 pm »

Las asignaciones de este tipo de números son muy lentas. Si puedes evitarte usar esa variable auxiliar aumentarías velocidad.

Y bueno el cout es muyyyyyyyy lento, tardarías menos en imprimir haciendo así:
Código
  1. #include <sstream>
  2.  
  3. //....
  4.  
  5. string valor = mpf_get_str(NULL, &buff, 10, 0, sumatoria);
  6. stringstream aux;
  7.  
  8. for (i = 1; i < digitos * 100 + 1; i++) {
  9.       aux << valor[i];
  10.       if (i % 10 == 0) aux << ' ';
  11.       if (i % 500 == 0) aux << endl;
  12.       if (i % 50 == 0 && i < digitos * 100) aux <<endl << i + 1 << ":\t";
  13. }
  14. aux << endl;
  15.  
  16. cout<<aux.str();
  17.  
En línea

Por favor, no me manden MP con dudas. Usen el foro, gracias.

¡Visita mi programa estrella!

Rar File Missing: Esteganografía en un Rar
engel lex
Moderador Global
***
Desconectado Desconectado

Mensajes: 15.514



Ver Perfil
Re: Calculo de pi en alta precisión (aporte)
« Respuesta #5 en: 9 Abril 2014, 20:57 pm »

en realidad no me preocupé del cout porque la diferencia ue causa es baja, si calculan 10000 decimales, haces unos 33000 calculos con 3 variables de 35kb xD e imprimir 35kb no toma sino unas décimas de segundo por lento que sea :p pero lo voy a hacer con printf me dió pereza hacerlo anoche

por otro lado no se puede hacet cifra a cifra ya que el resultado es la combinación de factoriales sumatorias y potencias, el numero siempre debe debe ser reusado
En línea

El problema con la sociedad actualmente radica en que todos creen que tienen el derecho de tener una opinión, y que esa opinión sea validada por todos, cuando lo correcto es que todos tengan derecho a una opinión, siempre y cuando esa opinión pueda ser ignorada, cuestionada, e incluso ser sujeta a burla, particularmente cuando no tiene sentido alguno.
amchacon


Desconectado Desconectado

Mensajes: 1.211



Ver Perfil
Re: Calculo de pi en alta precisión (aporte)
« Respuesta #6 en: 9 Abril 2014, 21:08 pm »

El problema no es imprimir muchos datos, sino el acceso a pantalla. Cada vez que adcedes a pantalla el programa se "congela" hasta que se refresque.

Por eso he condensado toda tu salida en un solo string, para que haga un solo acceso a pantalla.

PD: Si crees que la diferencia de tiempo es baja, te invito a poner un simple cout en el bucle de cálculo. A ver la diferencia de tiempos ;D
En línea

Por favor, no me manden MP con dudas. Usen el foro, gracias.

¡Visita mi programa estrella!

Rar File Missing: Esteganografía en un Rar
MCKSys Argentina
Moderador Global
***
Desconectado Desconectado

Mensajes: 5.129


Diviértete crackeando, que para eso estamos!


Ver Perfil
Re: Calculo de pi en alta precisión (aporte)
« Respuesta #7 en: 9 Abril 2014, 21:09 pm »

PD. Por ahora tengo la idea de hacer un colisionador md5  :silbar:

http://www.mscs.dal.ca/~selinger/md5collision/

Saludos!
En línea

MCKSys Argentina

"Si piensas que algo está bien sólo porque todo el mundo lo cree, no estás pensando."

z3nth10n


Desconectado Desconectado

Mensajes: 1.584


"Jack of all trades, master of none." - Zenthion


Ver Perfil WWW
Re: Calculo de pi en alta precisión (aporte)
« Respuesta #8 en: 9 Abril 2014, 21:12 pm »

Interesante el tema... Ya solo falta que haya un contador en la consola diciendo cuantos decimales se han calculado y los segundos requeridos :)

Un saludo.
En línea


Interesados hablad por Discord.
engel lex
Moderador Global
***
Desconectado Desconectado

Mensajes: 15.514



Ver Perfil
Re: Calculo de pi en alta precisión (aporte)
« Respuesta #9 en: 9 Abril 2014, 21:30 pm »

ikillnukes tambien a eso voy, segundo argumento sin salida, solo para ver el tiempo del calculo, los decimales calculados correctos no es tan simple fijate que yo genero decimales de mas, para asegurar el valor correcto, si ves.el metodo de euler ves que casi siempre da decimales periódicos

próximo a probar:salida concatenada a archivo, tiempo de ejecución, alguna sugerencia?
En línea

El problema con la sociedad actualmente radica en que todos creen que tienen el derecho de tener una opinión, y que esa opinión sea validada por todos, cuando lo correcto es que todos tengan derecho a una opinión, siempre y cuando esa opinión pueda ser ignorada, cuestionada, e incluso ser sujeta a burla, particularmente cuando no tiene sentido alguno.
Páginas: [1] 2 3 Ir Arriba Respuesta Imprimir 

Ir a:  

Mensajes similares
Asunto Iniciado por Respuestas Vistas Último mensaje
Bitrate + Motion Search Precision
Multimedia
G3N3S1S 2 1,816 Último mensaje 13 Enero 2005, 23:23 pm
por fanny
pequeño aporte(proxy),pero aporte al fin.:D
Programación Visual Basic
Tengu 0 1,768 Último mensaje 22 Julio 2007, 17:33 pm
por Tengu
Obtener más precisión que long double en C?¿?
Programación C/C++
jhonsc 3 2,827 Último mensaje 4 Febrero 2012, 13:48 pm
por Xandrete
Problema con la precisión del double
Java
danielo- 4 4,375 Último mensaje 27 Marzo 2012, 21:07 pm
por danielo-
Prueban en Atapuerca una tecnología inalámbrica de alta precisión pionera en...
Noticias
wolfbcn 0 883 Último mensaje 27 Junio 2012, 21:21 pm
por wolfbcn
WAP2 - Aviso Legal - Powered by SMF 1.1.21 | SMF © 2006-2008, Simple Machines