Foro de elhacker.net

Programación => Programación C/C++ => Mensaje iniciado por: engel lex en 9 Abril 2014, 08:16 am



Título: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex 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
(http://upload.wikimedia.org/math/6/3/1/63133aa55473a0e625a4684655831463.png)
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í (http://foro.elhacker.net/programacion_cc/duda_matematica_de_precision_arbitraria-t412080.0.html;msg1933292#msg1933292) 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

(http://upload.wikimedia.org/math/6/3/1/63133aa55473a0e625a4684655831463.png)

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

(http://upload.wikimedia.org/math/3/c/c/3cca05380343cb6f78d776d5e40b3a9d.png)

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


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: 1mpuls0 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:


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: ivancea96 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


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: dato000 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


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: amchacon 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.  


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex 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


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: amchacon 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


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: MCKSys Argentina 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/ (http://www.mscs.dal.ca/~selinger/md5collision/)

Saludos!


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: z3nth10n 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.


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex 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?


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: do-while en 10 Abril 2014, 00:06 am
¡Buenas!

Si te interesa el tema de pi, hay por ahí una formula de Ramanujan (http://es.wikipedia.org/wiki/Srinivasa_Aiyangar_Ramanujan) que calcula 8 decimales por iteración (era un fiera).

¡Saludos!


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex en 10 Abril 2014, 01:25 am
gracias do-while! :P vamos con esa! voy a extender el programa... ya llegué a mi casa, tengo 1,5l de te negro listos y 250g de maní japonés! ese código sale en un rato! :P


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: 1mpuls0 en 10 Abril 2014, 01:39 am
próximo a probar:salida concatenada a archivo, tiempo de ejecución, alguna sugerencia?

No sé mucho de c++, apenas estoy en pañales pero tal vez esto te pueda servir.
De hecho la idea viene de un programa que hice en otro lenguaje pero tratando de implementarlo en c++, aunque tal vez necesitas algo más "exacto" (microsegundos, aunque creo que los ordenadores no trabajan a ese nivel)

Además parece haber un fallo en mi método xD

Código
  1. #include <iostream>
  2. #include <ctime>
  3.  
  4. using namespace std;
  5.  
  6. int main() {
  7.  
  8.    clock_t tiempo_inicial = clock();
  9.  
  10.    //Alguna función
  11.    for(int indice=0;indice<=10000;indice=indice+1) {
  12.        std::cout << indice << endl;
  13.    }
  14.    clock_t tiempo_final = clock();
  15.  
  16.    std::cout << float( tiempo_final - tiempo_inicial ) /  CLOCKS_PER_SEC << " segundos";
  17.  
  18.    return 0;
  19. }
  20.  

Lo que me causa intriga es lo siguiente:  :¬¬
Citar
....
10000
1.283 segundos
Process returned 0 (0x0)   execution time : 1.323 s


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex en 10 Abril 2014, 04:51 am
la  ecuación de Ramanujan es una locura @.@ es un poco dificil de aplicar... sin embargo dudo de su eficiencia al momento de procesar, ya que a diferencia de la de euler es muy compleja, tiene 2 factoriales diferentes (y el de arriba es grande, con k = 3 el valor es de 479.001.600) la de abajo tiene 2 potencias y la segunda es de 4*k, haciéndola sin funciones matemáticas precompiladas sería

Código
  1. for(i=1;i<precision;i++){
  2. for(j=1;j<=4*i;j++){
  3.  factorial1*=j;
  4. }
  5. for(j=1;j<=i;j++){
  6.  factorial2*=j;
  7. }
  8. for(j=1;j<=4;j++){
  9.  factorial2*=factorial2;
  10. }
  11. for(j=1;j<=4*i;j++){
  12.  divisor*=396;
  13. }
  14. sumatoria+=(factorial1*(1103+26390*i))/(factorial2*divisor)
  15.  
  16. }
  17. sumatoria=1/sumatoria
  18.  

el código sobre eso se puede optimizar pero igual son muchos ciclos dentro de cada ciclo (cuando precision = 1.000; tendrá que hacer por lo menos 4.000.000 de ciclos secundarios los cuales multiplicarán valores muy grandes, esto sin contar que sería solo para obtener 8.000 dígitos teóricos...)


la de euler es simple y hasta bonita XD solo requiere en ciclo principal y unas pocas operaciones

aqui la formula de Ramanujan
(http://upload.wikimedia.org/math/3/c/c/3cca05380343cb6f78d776d5e40b3a9d.png)

cambié un poco el programa para hacerlo más modular...

ya se me hizo tarde... así que ni argumentos (más que ramanujan), ni corrida silenciosa, ni salida a archivo! :P sorry por ahora... estoy cansado  :(

agregado por ahora, ramanujan (no da pi... si alguien me ayuda con esto, sobre donde me equivoco...), tiempo de calculo e impresion


por cierto amchacon sobre el tiempo de impresión incluso a 40.000 decimales la diferencia es instrumental, no llega al medio segundo de diferencia :P

aquí el código... se está tornando largo D:
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.  
  9. using namespace std;
  10. string metodo_euler(unsigned long int digitos, mp_exp_t &exponente);
  11. string metodo_ramanujan(unsigned long int digitos, mp_exp_t &exponente);
  12. void imprimir_pi_aux(string pi, mp_exp_t exponente, int digitos);
  13.  
  14. int main(int argc, char **argv) {
  15.  
  16.    unsigned long int digitos = strtol(argv[1], NULL, 10); //argv a int
  17.    mpf_set_default_prec(32 * 10 * (digitos / 100) * 1.11); //precision + 10%
  18.    long exponente;
  19.    clock_t start = clock();
  20.    string pi;
  21.    if(argc>2 && argv[2]=="ramanujan"){
  22.        pi = metodo_ramanujan(digitos, exponente);
  23.    }else{
  24.        pi = metodo_euler(digitos, exponente);
  25.    }
  26.    clock_t end = clock();
  27.    imprimir_pi_aux(pi, exponente, digitos);
  28.    cout << "\ntiempo total de ejecucion: " << (float) (end - start) / CLOCKS_PER_SEC << "\n";
  29.    cout << "tiempo total de impresion: " << (float) (clock() - end) / CLOCKS_PER_SEC << "\n\n";
  30.    return 0;
  31. }
  32.  
  33. void imprimir_pi_aux(string pi, long exponente, int digitos) {
  34.    cout << "PI: ";
  35.    stringstream aux;
  36.    unsigned long int i;
  37.    if (exponente < 0) {
  38.        cout << "0.";
  39.        while (exponente < 0) {
  40.            cout << 0;
  41.            exponente++;
  42.        }
  43.    } else {
  44.        for (i = 0; i < exponente; i++) {
  45.            cout << pi[i];
  46.        }
  47.        cout << ".\n\n1:\t";
  48.    }
  49.    for (i = exponente; i < digitos + 1; i++) {
  50.        aux << pi[i];
  51.        if (i % 10 == 0) aux << ' ';
  52.        if (i % 500 == 0) aux << endl;
  53.        if (i % 50 == 0 && i < digitos) aux << endl << i + 1 << ":\t";
  54.    }
  55.    aux << endl;
  56.    cout << aux.str();
  57. }
  58. string metodo_euler(unsigned long int digitos, long &exponente) {
  59.    digitos /= 100; //trabaja con lotes de 100
  60.    unsigned long int precision = digitos * 100 * 3 * 1.12;
  61.    unsigned long int x, i;
  62.    mpf_t factorial, sumatoria, buff_factorial; //variables GMP
  63.    mpf_init_set_ui(factorial, 1); //factorial=1
  64.    mpf_init_set_ui(sumatoria, 1); //sumatoria=1
  65.    mpf_init_set_ui(buff_factorial, 1); //buff_factorial=1
  66.  
  67.    for (x = 1; x < precision; x++) {
  68.        i = 2 * x + 1;
  69.        mpf_set_ui(buff_factorial, x); //buff_factorial=x
  70.        mpf_div_ui(buff_factorial, buff_factorial, i); //buff_factorial*=i
  71.        mpf_mul(factorial, factorial, buff_factorial); //factorial*=buff_factorial
  72.        mpf_add(sumatoria, sumatoria, factorial); //sumatoria+=factorial
  73.    }
  74.    mpf_mul_ui(sumatoria, sumatoria, 2); //sumatoria*=2
  75.    return mpf_get_str(NULL, &exponente, 10, 0, sumatoria);
  76. }
  77.  
  78. string metodo_ramanujan(unsigned long int digitos, long &exponente) {
  79.    unsigned long int i, j;
  80.    mpf_t sumatoria, buff_sumatoria, buff_ciclo_superior, buff_ciclo_inferior, primera_parte,
  81.            buff_ciclo_inferior2, factorial_superior, factorial_inferior;
  82.    mpf_init_set_ui(sumatoria, 1103); //sumatoria=0
  83.    mpf_init_set_ui(buff_sumatoria, 0); //buff_sumatoria=0
  84.    mpf_init_set_ui(buff_ciclo_superior, 0); //buff_ciclo_superior=0
  85.    mpf_init_set_ui(buff_ciclo_inferior, 0); //buff_ciclo_inferior=0
  86.    mpf_init_set_ui(primera_parte, 0); //primera_parte=0
  87.    mpf_init_set_ui(buff_ciclo_inferior2, 0); //buff_ciclo_inferior2=0
  88.    mpf_init_set_ui(factorial_superior, 1); //factorial_superior=1
  89.    mpf_init_set_ui(factorial_inferior, 1); //factorial_inferior=1
  90.    mpf_sqrt_ui(primera_parte, 2); //primera_parte=sqrt(2)
  91.    mpf_mul_ui(primera_parte, primera_parte, 2); //primera_parte*=2
  92.    mpf_div_ui(primera_parte, primera_parte, 9801); //primera_parte/=9801
  93.  
  94.    for (i = 1; i <= digitos; i++) {
  95.        mpf_set_ui(factorial_superior, 1);
  96.        for (j = 2; j <= i * 4; j++) {
  97.            mpf_mul_ui(factorial_superior, factorial_superior, j);
  98.        }
  99.        mpf_set_ui(buff_ciclo_superior, 26390);
  100.        mpf_mul_ui(buff_ciclo_superior, buff_ciclo_superior, i);
  101.        mpf_add_ui(buff_ciclo_superior, buff_ciclo_superior, 1103);
  102.        mpf_mul(buff_ciclo_superior, buff_ciclo_superior, factorial_superior);
  103.  
  104.        mpf_mul_ui(factorial_inferior, factorial_inferior, i);
  105.        mpf_pow_ui(buff_ciclo_inferior, factorial_inferior, 4);
  106.        mpf_set_ui(buff_ciclo_inferior2, 396);
  107.        mpf_pow_ui(buff_ciclo_inferior2, buff_ciclo_inferior2, 4 * i);
  108.        mpf_mul(buff_ciclo_inferior, buff_ciclo_inferior, buff_ciclo_inferior2);
  109.  
  110.        mpf_div(buff_sumatoria, buff_ciclo_superior, buff_ciclo_inferior);
  111.        mpf_add(sumatoria, sumatoria, buff_sumatoria);  
  112.    }
  113.    mpf_mul(sumatoria, sumatoria, primera_parte);
  114.    mpf_ui_div(sumatoria, 1, sumatoria);
  115.    return mpf_get_str(NULL, &exponente, 10, 0, sumatoria);
  116. }


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: do-while en 10 Abril 2014, 11:22 am
¡Buenas!

Te ha faltado calcular el factorial inferior:
Tu código:
Código
  1.        mpf_mul_ui(factorial_inferior, factorial_inferior, i); //linea 104
  2.       mpf_pow_ui(buff_ciclo_inferior, factorial_inferior, 4);
  3.  

Corregido:
Código
  1.        mpf_set_ui(factorial_inferior, 1);
  2.       for (j = 2; j <= i ; j++) {
  3.           mpf_mul_ui(factorial_inferior, factorial_inferior, j);
  4.       }
  5.       mpf_pow_ui(buff_ciclo_inferior, factorial_inferior, 4);
  6.  

Creo que eso era todo.

¡Saludos!


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex en 10 Abril 2014, 14:17 pm
no, el factorial inferior si se calcula, lo hago igual que el de euler, es decir como es un factorial que crece de 1 en 1 uso el ciclo principal para evitar hacer más ciclos

fíjate que el factorial superior inicia reasignándose y el inferior no... pero me diste una idea para hacer algo similar con el superior


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: do-while en 10 Abril 2014, 14:44 pm
no, el factorial inferior si se calcula, lo hago igual que el de euler, es decir como es un factorial que crece de 1 en 1 uso el ciclo principal para evitar hacer más ciclos


El factorial inferior lo tienes que recalcular en cada iteracion. No es el factorial de n, siendo n el numero de iteraciones. Para cada iteración k es el factorial de k. Si dices que no te da el resultado, en algún punto tiene que haber un error. ¿No?

Además, de la forma en la que lo haces no estas calculando el factorial de nada. Estas multiplicando un acumulador por el numero de ciclo y luego lo elevas a 4, lo multiplicas por el siguiente contador de ciclo y lo vuelves a elevar a 4. Eso no es lo que indica la fórmula...


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex en 10 Abril 2014, 14:59 pm
mira como lo hice en euler :P

si tengo un error, pero hice ya el debug porque justamente se que esa es la principal posible primera falla :P


factorial_inferior = 1 // 0!
ciclo 1:
factorial_inferior *=1 //1*1 = 1
ciclo 2:
factorial_inferior *=2 //2*1*1 = 2
ciclo 3:
factorial_inferior *=3 //3*2*1*1 = 6
ciclo 4:
factorial_inferior *=4 //4*3*2*1*1 = 24
ciclo 5:
factorial_inferior *=5 //5*4*3*2*1*1 = 120

el arrastra el valor del factorial y lo multiplica cada vez, por eso no le reasigno un valor entero :P (así evito hacer un ciclo extra con una variable de 1kb o más) el superior lo voy a hacer similar corriéndolo solo de i*4 y x< (i+1)*4 así solo hace 4 ciclos por ciclo :P pero esa parte si está funcionando... comprobado (hace 15 segundos XD)

sospecho que tendré que hacer debug con buena calma ver los valores con lápiz en mano... porque en donde sospecho que puede salir el error es en uno de esos números grandes


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: Gh057 en 11 Abril 2014, 06:37 am
buenas! entré un segundo justo antes de ir a dormir, ando complicado con los tiempos jejej y recién veo el tema! pero no quería dejar el sitio sin antes felicitarte engel lx por el post, es muy interesante!

tan solo a titulo informativo, el excelente algoritmo de Ramanujan ha sido modificado ligeramente obteniéndose el algoritmo de Chudnovsky, que permite algo así como 15 decimales de exactitud creo, luego de cada término XD una locura!

(http://s.wordpress.com/latex.php?latex=%20\cfrac{1}{\pi}%20%3D%2012%20\%3B%20\displaystyle{\sum^\infty_{k%3D0}%20\cfrac{%28-1%29^k%20%286k%29!%20%2813591409%20%2B%20545140134k%29}{%283k%29!%28k!%29^3%20640320^{3k%20%2B%203%2F2}}}&bg=ffffff&fg=000000&s=0)

pd: acabo de buscarlo, son 14 decimales; de 9 que se obtenía con el anterior. tengo entendido que los últimos records que indicaban más arriba fueron con dicha implementación. saludos


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex en 11 Abril 2014, 07:26 am
do-while: lo que es la arrogancia! XD al final si era el factorial el error! XD (eso creo)
XD pero ahorita está más feo aún pero funciona! y es increiblemente eficiente!

de ahora en adelante publicaré el codigo en el post inicial...  no lo publicaré como un gran bloque sino, cada sección... olvidense de imprimir en archivo y de la modularidad! XD eso aún lo postergo

Gh057: lo había visto! XD pero me tomará otro dia adaptarlo, aunque como tiene similaridad con Ramanujan no se me debría ser tan difícil...


alguna vez vi algo de comentario de código cómicos... este merece uno que vi

Código
  1. //cuando empeze a hacer esto solo dios y yo sabiamos lo que hacia, en este momento solo dios lo sabe

aqui la correccion de Ramanujan

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

si no fuera tan flojo (y estuviera claro en el tema), intentara hacer operator overload para usar los operadores normales! D:

a partir de ahora modifico el post principal


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: scott_ en 14 Agosto 2014, 02:33 am
do-while: lo que es la arrogancia! XD al final si era el factorial el error! XD (eso creo)
XD pero ahorita está más feo aún pero funciona! y es increiblemente eficiente!

de ahora en adelante publicaré el codigo en el post inicial...  no lo publicaré como un gran bloque sino, cada sección... olvidense de imprimir en archivo y de la modularidad! XD eso aún lo postergo

Gh057: lo había visto! XD pero me tomará otro dia adaptarlo, aunque como tiene similaridad con Ramanujan no se me debría ser tan difícil...


alguna vez vi algo de comentario de código cómicos... este merece uno que vi

Código
  1. //cuando empeze a hacer esto solo dios y yo sabiamos lo que hacia, en este momento solo dios lo sabe

aqui la correccion de Ramanujan

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

si no fuera tan flojo (y estuviera claro en el tema), intentara hacer operator overload para usar los operadores normales! D:

a partir de ahora modifico el post principal

Disculpa por revivir el post. Muy buen post, en definitiva Pi tiene muchos decimales, que en calculos escolares me basta con 3.1416 xD, Pero jamás entendí la ecuacuón de Ramanujan.
Me volví loco y opté por enender otras, ¿me la puedes explicar con peras y manzanas?  :P

Un Cordial Saludo  ;D


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex en 14 Agosto 2014, 02:48 am
Disculpa por revivir el post. Muy buen post, en definitiva Pi tiene muchos decimales, que en calculos escolares me basta con 3.1416 xD, Pero jamás entendí la ecuacuón de Ramanujan.
Me volví loco y opté por enender otras, ¿me la puedes explicar con peras y manzanas?  :P

Un Cordial Saludo  ;D


que cosa no entendiste? las operaciones en el código?


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: scott_ en 14 Agosto 2014, 02:50 am

que cosa no entendiste? las operaciones en el código?

Si  :huh:


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex en 14 Agosto 2014, 03:13 am
lo voy a traducir a codigo en lenguaje "normal"(c++) para que sea más visible

recordaré la formula por cosas paracticas
(http://upload.wikimedia.org/math/3/c/c/3cca05380343cb6f78d776d5e40b3a9d.png)

Código
  1. unsigned long int i, j;
  2.  
  3. float sumatoria, buff_sumatoria, buff_ciclo_superior, buff_ciclo_inferior, primera_parte, buff_ciclo_inferior2, factorial_superior, factorial_inferior;
  4.  
  5. sumatoria=1103;
  6. buff_sumatoria=0;
  7. buff_ciclo_superior=0;
  8. buff_ciclo_inferior=0;
  9. primera_parte=0;
  10. buff_ciclo_inferior2=0;
  11. factorial_superior=1;
  12. factorial_inferior=1;
  13.  
  14. //sqrt(2)
  15. primera_parte=sqrt(2);
  16.  
  17. //sqrt(2) * 2
  18. primera_parte*=2;
  19.  
  20. //sqrt(2)*2 / 9801
  21. primera_parte/=9801
  22.  
  23. for (i = 1; i <= digitos; i++) {//sabiendo que ramanujan produce 8 digitos con un margen de error de 1% por ciclo, según mis pruebas
  24.        for (j = (i - 1)*4 + 1; j <= i * 4; j++) {
  25.            factorial_superior*=*j;//ciclo acumulador del factorial (4k)!
  26.        }
  27.        //29690 * k
  28.        buff_ciclo_superior= 26390;
  29.        buff_ciclo_superior *= i;
  30.  
  31.        //(1103 + 29690*k)
  32.        buff_ciclo_superior+=1103;
  33.  
  34.        //(4k)! * (1103+29690*k)
  35.        buff_ciclo_superior*=factorial_superior;
  36.  
  37.        //(k!)^4
  38.        factorial_inferior*=i;
  39.        buff_ciclo_inferior = pow(factorial_inferior, 4);
  40.  
  41.        //396^4k
  42.        buff_ciclo_inferior2 = 396;
  43.        buff_ciclo_inferior2 = pow(buff_ciclo_inferior2, 4 * i);
  44.  
  45.        //(k!)^4 * 396^4
  46.        buff_ciclo_inferior*=buff_ciclo_inferior2;
  47.  
  48.        // (4k)!*(1103+29690*k) / (k!)^4*396^4
  49.        buff_sumatoria = buff_ciclo_superior / buff_ciclo_inferior;
  50.  
  51.        sumatoria+=buff_sumatoria;//la sumatoria...
  52.  
  53.    }
  54.  
  55.    //sqrt(2)*2/9801 * SUMATORIA (4k)!*(1103+29690*k) / (k!)^4*396^4
  56.    sumatoria*=primera_parte;
  57.  
  58.    // el resultado original es 1/pi, entonces con esto hago la inversión
  59.    sumatoria = 1 / sumatoria;
  60.  
  61.    //convierte a texto para poder imprimirlo
  62.    return mpf_get_str(NULL, &exponente, 10, 0, sumatoria);
  63.  
  64.  
  65.  

espero que así sea más visible


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: scott_ en 14 Agosto 2014, 03:14 am
Mejor  ;D ;D ;D

Que lenguaje C++ si se una gran parte de él  ;-)

Gracias y un Cordial Saludo  ;D


Título: Re: Calculo de pi en alta precisión (aporte)
Publicado por: engel lex en 14 Agosto 2014, 03:20 am
Mejor  ;D ;D ;D

Que lenguaje C++ si se una gran parte de él  ;-)

Gracias y un Cordial Saludo  ;D

es que siempre fue en c++ solo que para poder manejar grandes cifras (al final manejaba 10 millones de digitos (en menos de 5 minutos creo...) descargate la librería y prueba el código e intenta a ver si puedes hacer el que no pude, el de Chudnovsky que lo intenté y se me hizo muy complicado (luego estuve muy saturado de trabajo en la oficina para pensar en esto XD)

(http://s.wordpress.com/latex.php?latex=%20\cfrac{1}{\pi}%20%3D%2012%20\%3B%20\displaystyle{\sum^\infty_{k%3D0}%20\cfrac{%28-1%29^k%20%286k%29!%20%2813591409%20%2B%20545140134k%29}{%283k%29!%28k!%29^3%20640320^{3k%20%2B%203%2F2}}}&bg=ffffff&fg=000000&s=0)