Mostrando entradas con la etiqueta C. Mostrar todas las entradas
Mostrando entradas con la etiqueta C. Mostrar todas las entradas

viernes, 1 de marzo de 2013

El factorial de las mil y una noches.

Leyendo en Internet algunas cosas relacionadas con lo que se trata en la entrada anterior llegué a una de las "delicias digitales endiabladamente difíciles" (en palabras del autor). Se trata de una cuestión incluida en el capítulo 37 del libro La maravilla de los números, de Clifford A. Pickover. Así que desempolvo el libro, que lo tengo en mi estantería.


Incluyo aquí el fragmento al que me refiero:

"Desde la edad de trece años, cada 1001 días el doctor Googol lee Las mil y una noches, lo que significa que lee esa obra una vez cada 2,74 años. Con la excepción del Corán, ninguna otra obra literaria árabe se conoce mejor ni tiene más influencia en Occidente que Las mil y una noches. Esta colección de relatos está agrupada en derredor de uno central relativo a un sultán y a sus amantes. Después de descubrir que su esposa le había sido infiel, el sultán hace la promesa de tomar una nueva novia cada día y hacerla ejecutar al día siguiente.

Cuando Sahrazad fue elegida como su nueva esposa, cada atardecer ella le contaba un cuento al sultán, pero no lo terminaba, prometiendo hacerlo a la siguiente noche si sobrevivía. Esto continuó durante mil y una noches, hasta que el sultán quedó profundamente enamorado de Sahrazad y olvidó sus crueles planes de ejecución.

Una noche, después de haber leído Las mil y una noches, el doctor Googol empezó a preguntarse acerca de un número especial llamado el factorial de las mil y una noches. Este número está definido con el valor de x tal que x! tiene 1001 dígitos. Los factoriales crecen bastante rápidamente: 5!=120, 10!=3.682.800 y 15!=1.307.674.368.000. ¿Cuál es el factorial de las mil y una noches?"


Aprovechando lo que expliqué en la entrada anterior, lo que estamos buscando es el menor número x tal que x! tenga 1001 cifras, entonces, utilizando la expresión cifras de x! = [log(x!)] + 1 obtenemos:
1001 = [log(x!)]+1 que es equivalente a [log(x!)]=1000.

Como log(x!)=log(1)+log(2)+...+log(x) lo que podemos hacer es ir sumando los logaritmos decimales de los números naturales hasta que la parte entera del resultado sea 1000 (teniendo en cuenta que podríamos pasarnos de largo). O equivalentemente, hasta sumar el logaritmo decimal del número que haga que el resultado sea (por primera vez) mayor que 1000.

Con pequeñas modificaciones al programa de la entrada anterior construimos el siguiente código:


#include < stdlib.h >
#include < iostream >
#include < math.h >
main()
{
    double r=0;
    unsigned long int i=2, numero=1001;
    while(r < numero-1)
    {
        r+=log10(i);
        i++;
    }
    i--;  // quitamos el último incremento
    std::cout < < "El primer factorial cuyo resultado tiene por lo menos " < < numero < < " cifras es " < < i < < "!." < < std::endl;
    return 0;
}

La ejecución del programa anterior devuelve:
El primer factorial cuyo resultado tiene por lo menos 1001 cifras es 450!.

Podría ser que 450! tuviera más de 1001 cifras, dado que la condición que imponemos es que por lo menos tenga esa cantidad de cifras. Si ejecutamos el programa de la entrada anterior obtenemos:
El factorial de 450 tiene 1001 cifras.
O si alguien no se fía de mis dotes de programador, también puede preguntárselo a Wolfram Alpha :-)

Problema resuelto. Bueno, este en particular y el general de hallar el menor número x cuyo factorial tenga por lo menos n cifras (siempre que el tamaño de n no desborde la variable).


Curiosidad: Al igual que en la obra El prodigio de los números, Pickover dedica el libro al cuadrado mágico apocalíptico.

Clifford A. Pickover

¿Cuántas cifras tiene un número? Por ejemplo, el factorial de 1.000.000

El sistema de numeración que utilizamos es posicional y decimal, es decir que utilizando 10 dígitos (0, 1, 2, 3, 4, 5, 6, 7, 8 y 9) podemos escribir cualquier número entero.

¿Cuántas cifras tiene el número 92.368.556? Una pregunta bastante sencilla. Contamos y listo. 8 cifras.

¿Cuántas cifras tiene el factorial de 1.001? Pues ahora la cosa no es tan sencilla. Calcular el factorial de 1.001 y contar sus cifras es una posibilidad, pero en vista de que el número es bastante elevado vamos a intentar encontrar una alternativa que nos evite calcular el factorial.

Si volvemos al hecho de que usamos un sistema decimal, podemos concluir fácilmente que se añade una cifra más cada vez que damos un salto a la siguiente potencia de 10.
En 10¹ = 10 comienzan los números de 2 cifras.
En 10² = 100 comienzan los números de 3 cifras.
En 10³ = 1.000 comienzan los números de 4 cifras
... si generalizamos esta propiedad, podemos formularla como ...
En 10^n comienzan los números de n+1 cifras.

Si el cambio de cifras se produce en las potencias de 10, para conocer cuántas cifras tiene un número podemos utilizar la operación inversa a la exponenciación de base 10: los logaritmos decimales.


Por ejemplo, si tenemos en cuenta que log(10)=1 y log(100)=2, al calcular el logaritmo decimal de cualquier número entre 10 y 100 obtenemos, dado que el logaritmo es una función estrictamente creciente, un resultado decimal entre 1 y 2.

Entonces, la cantidad de cifras de un número cualquiera se obtiene sumando 1 a la parte entera de su logaritmo decimal.

Ejemplo: ¿Cuántas cifras tiene el número 92.368.556?
log(92368556) = 7,96552415434
La parte entera, a partir de ahora denotada por [], es [7,96552415434]=7.
Por tanto, el número 92.368.556 tiene 8 (7+1) cifras.

Si aplicamos esto mismo a la pregunta de cuántas cifras tiene 1001!, tendríamos:
cantidad de cifras de 1001! = [log(1001!)]+1

1001! es un número elevado para las calculadoras, pero wxMaxima puede calcular la expresión [log(1001!)]+1 directamente con el siguiente comando, sin que se produzca ningún desbordamiento:
entier(log(1001!)/log(10))+1;
Téngase en cuenta que en wxMaxima la función log corresponde al logaritmo neperiano y no al decimal, por eso aparece en la expresión la división entre log(10). El resultado que devuelve wxMaxima es:
2571
Bien, ya sabemos cuántas cifras tiene el factorial de 1001, pero además de que hemos hecho la gran trampa de que wxMaxima sí ha calculado el factorial de 1001, ¿qué pasa si el número es mayor? Pongamos por ejemplo que queremos saber cuántas cifras tiene 1.000.000!
Aquí wxMaxima ya no arroja resultados. Podríamos recurrir a la programación para definir una función que calcule el factorial multiplicando tipos enteros, pero la mayoría de lenguajes de programación tiene un rango limitado para sus tipos de variables. Por ejemplo, un tipo unsigned long long en el lenguaje de programación C tiene un rango de 0 a 18.446.744.073.709.551.615. Cualquier operación cuyo resultado sea mayor provoca un desbordamiento. Otros lenguajes de programación, como python, ya llevan implementadas en sus librerías funciones que efectúan las operaciones de este tipo sin provocar desbordamiento.
 
En cualquier caso, lo que estamos pensando es si puedo evitar el cálculo del factorial porque lo único que me interesa es saber cuántas cifras tiene. Así que volvemos a utilizar la expresión general:
cifras de n! = [log(n!)]+1

Nuestro problema era que no queremos/podemos calcular el factorial y si observamos la propiedad anterior resulta que ahora no sólo tengo que calcular un factorial sino que además tengo que calcular un logaritmo. Menudo invento. Bueno, tranquilidad. El factorial son productos y sabemos que los logaritmos tienen la propiedad de "transformar las multiplicaciones en sumas": log(a·b)=log(a)+log(b).
Así pues,
log(1000000!) = log(1000000·999999·...·2·1)=log(1000000)+log(999999)+...+log(2)+log(1).

¿Qué ganamos con esto? Pues reducir considerablemente el orden de los números con los que se trabaja. Los números decimales que aparecen en los sumandos de log(1000000)+...+log(2)+log(1) están comprendidos entre 0 y 6 (incluidos). La precisión de los tipos decimales de cualquier lenguaje de programación es más que suficiente para el cálculo anterior, especialmente si tenemos en cuenta que lo que nos interesa ahora es la parte entera del resultado del sumatorio.


Por ejemplo, en C++ podemos escribir el siguiente programa:



#include < stdlib.h >
#include < iostream >
#include < math.h >
main()
{
    double r=0;
    unsigned long int numero=1000000;
    for(unsigned long int n=1; n<=numero; n++)
    {
        r+=log10(n);
    } //for n
    unsigned long int cifras_factorial= (unsigned long int) r;
    cifras_factorial++;
    std::cout < < "El factorial de " < < numero < < " tiene " < < cifras_factorial < < " cifras." < < std::endl;
    return 0;
}


La ejecución del programa anterior devuelve la siguiente sentencia:
El factorial de 1000000 tiene 5565709 cifras.

Como curiosidad final, hice un programa en python que calcula todas las cifras del factorial de 1000000.


Entrada relacionada:
- El factorial de las mil y una noches.

sábado, 17 de noviembre de 2012

La historieta de Martin y Gala: lo esperado no es tan probable.

Tras contar la historieta de Martin y Gala dediqué dos entradas a mostrar dos posibles resoluciones de la misma: desde un enfoque probabilístico y utilizando una simulación por ordenador. Ambas sendas nos conducían a la solución de que se espera que el invitado tenga que contestar 255 preguntas de Martin antes de poder pasar a Gala.

Pero, ¿existen muchas probabilidades de que sean exactamente 255 preguntas las que conteste y acierte en mi previsión? Y es que, lo esperado no necesariamente coincide con lo más probable. Son conceptos que responden a preguntas diferentes.

Imaginemos un dado trucado de manera que la probabilidad de cada una de sus caras es 1/7, salvo la probabilidad de la cara con un 6 marcado que es igual a 2/7. La esperanza matemática de dicho dado es:
1·(1/7) + 2·(1/7) + 3·(1/7) + 4·(1/7) + 5·(1/7) + 6·(2/7) = 15/7 + 12/7 = 27/7 que redondeando a las centésimas queda 3'86.
Es decir, que si jugamos muchas (pero que muchas) veces, la media de las puntuaciones obtenidas en las tiradas tiende a 3,86. Pero eso no significa que el 4 sea el valor más probable (de hecho en nuestra situación imaginaria el más probable es el 6).

Volvemos a nuestra historia de Martin y Gala. ¿Cuál es la probabilidad de que el invitado tenga que contestar exactamente 255 preguntas antes de pasar a Gala? ¿Cuál es la probabilidad de que el invitado se libre de ese calvario y pueda pasar antes de contestar 255 preguntas?

La situación de la primera pregunta, la probabilidad de contestar exactamente 255 preguntas, se da cuando se pasa a Gala en el juego número 256, por lo que la probabilidad es (con p=(1/2)^8 y q=1-p):
p(contestar exactamente 255 preguntas) = q^255 · p = 0'00143984217604524.

Vemos que la probabilidad de contestar exactamente 255 preguntas es bastante pequeña, no llega ni al 15 .

Y para abordar la segunda pregunta, la probabilidad de tener que contestar menos preguntas, podemos plantearlo de la siguiente manera:
p(contestar menos de 255 preguntas) = p(contestar exactamente 0 preguntas) + p(contestar exactamente 1 pregunta) + ... + p(contestar exactamente 254 preguntas) =
= p + q · p + q^2 · p + ... + q^254 · p = p · (1+q+q^2+...+q^254) = p · (q^255 - 1) / (q - 1) = 1 - q^255 = 0'6314004029324185.


Es decir, hay más de un 63% de probabilidad de que el invitado pase a Gala antes de contestar las 255 preguntas que se espera que conteste.

¿Y la simulación por ordenador qué nos dice? Añadimos unas pocas líneas más (ver código al final de la entrada) y ejecutamos el programa:
Media de preguntas antes de pasar: 255'193
Frecuencia relativa de coincidencia con preguntas esperadas: 0'001457
Frecuencia relativa de menos preguntas de las esperadas: 0'63112

La aproximación que devuelve una simulación de 1.000.000 de repeticiones es bastante buena.


Podemos concluir que lo más probable es que el invitado tenga que contestar menos preguntas de las "esperadas".


//////////////////////////////////////////////////////////////////////////////////////////////////////
Código fuente del programa en C (en negrita lo añadido):

#include < stdlib.h >
#include < iostream >
#include < time.h >
#include < math.h >
 
main()
{
    srand(time(NULL));
    unsigned int puntos=0, aleatorio=0, puntosNecesarios=8;
    unsigned long int preguntas=0, sumaPreguntas=0, vecesSimulacion=1000000;
    unsigned long int preguntasEsperadas=pow(2,puntosNecesarios)-1, coincidePreguntasEsperadas=0, menorPreguntasEsperadas=0;
    for(unsigned int n=0; n < vecesSimulacion; n++)
    {
        puntos=0;
        preguntas=0;
        while(puntos < puntosNecesarios)
        {
            aleatorio=rand()%2;
            if(aleatorio==0)
            {
                puntos=0;
                preguntas++;
            }
            else
            {
                puntos++;
            }
        } //while
        if(preguntas==preguntasEsperadas)
        {
            coincidePreguntasEsperadas++;
        }
        if(preguntas < preguntasEsperadas)
        {
            menorPreguntasEsperadas++;
        }

        sumaPreguntas+=preguntas;
    } //for n
    std::cout < < "Media de preguntas antes de pasar: " < < (float)sumaPreguntas/(float)vecesSimulacion < < std::endl;
    std::cout < < "Frecuencia relativa de coincidencia con preguntas esperadas: " < < (float)coincidePreguntasEsperadas/(float)vecesSimulacion < < std::endl;
    std::cout < < "Frecuencia relativa de menos preguntas de las esperadas: " < < (float)menorPreguntasEsperadas/(float)vecesSimulacion < < std::endl;

    return 0;
}

//////////////////////////////////////////////////////////////////////////////////////////////////////