4

Tengo un problema con un programa que he hecho para resolver un problema matemático, que consiste en encontrar números enteros y que cumplan que el resultado x = sqrt(1141.0*y^2+1.0) sea un número entero. Si no estoy equivocado la máxima precisión se obtiene con long double como formato de variable, que es el que he usado en el programa, y, en teoría, alcanza 10^308, sin embargo como no me funcionaba correctamente he hecho una depuración del programa y me sale que sólo está manejando cifras de 18 dígitos, por lo tanto me está devolviendo como posibles soluciones valores de y que realmente no son soluciones.

Por ejemplo me devuelve como posible solución el valor y=23052875, cuando para ese valor de y la x=7.7869595299999994670577156…×10^8 realmente, por lo que se ve que C está redondeando el resultado.

¿Es una limitación de C, o de la función sqrt, y, por lo tanto, no se puede resolver este problema con C? o ¿estoy utilizando mal algo y eso me está limitando la precisión?

#include <stdio.h>
#include <math.h>

int is_int(long double x);

int main() {

long double x;
long int y, MAX=5000000;
printf("\n\tSolutions: ");

for (y = 1; y < MAX; y++ ) {
    x = sqrt(1141.0*y*y+1.0);
    if (is_int(x)) {
        printf("\t%ld", y);
    }
}
printf("\n\n");
return 0;
}

int is_int (long double x) {

if (floor(x) != x) {
    return 0;
}
else{
        return 1;
    }
}

(aclaración sobre el programa: la variable MAX tiene un valor elegido arbitrario, pero se puede cambiar sin problemas para realizar pruebas)

2
  • Ya está insertado el código en el mensaje. Era por no hacerlo muy larga. Un saludo.
    – Vacendak
    Commented el 18 jul. 2018 a las 15:49
  • Eliminado. Algo más? Y en cuanto al argumento de sqrt , si te refieres a los decimales, como podrás observar en el código están con decimales, y si te refieres a otra cosa, si especificas un poco igual me entero de a qué te refieres. Gracias.
    – Vacendak
    Commented el 18 jul. 2018 a las 15:54

1 respuesta 1

5

Problema

El problema no es del C, sino del formato de coma flotante.

En un formato de coma flotante se usan parte de los bits para contener las "cifras significativas" de la cantidad a representar (en binario) y otra parte de los bits para contener un "exponente". Esto permite que, cuando se trabaja con números muy pequeños se pueda tener mucha precisión (pues los exponentes serían negativos), y que cuando se trabaja con números grandes se pueda tener mucho rango, al ser los exponentes positivos y grandes, pero a cambio de perder precisión.

En un double, efectivamente puedes representar número tan grandes como 10^308, pero eso no implica que dispongas de 308 cifras significativas. Cuanto más grande es el número, mayor es el error de representación, y una gran parte de las cifras del mismo se sustituyen por ceros (en binario), por lo que en realidad no estás representando el número que pensabas, sino una aproximación del mismo.

En tu caso, cuando evalúas la expresión 1141*y*y+1 para y=23052875, el resultado sería (manteniendo todos los dígitos significativos) 601053036760921876 (lo he calculado con python, que tiene precisión arbitraria con enteros). En binario ese número es 100001010111010111101000110101101011010011011111011100010100.

En cambio al representarlo como un double, todos esos dígitos no caben (son 62 bits, y un double sólo tiene 51 bits para la mantisa). Por tanto parte de ellos se eliminan (se consideran cero), y estarás representando entonces la aproximación 601053036760921856, que en binario es 100001010111010111101000110101101011010011011111011100000000, distinto del que queríamos en las últimas 5 cifras binarias.

A partir de ahí, lógicamente, todo saldrá mal.

Solucion

Creo que para este problema, que además parece matemáticamente un problema de aritmética entera, el tipo apropiado sería el entero.

Python tiene enteros de precisión arbitraria (es decir, la cantidad de bits que dedica para guardar un entero, puede crecer tanto como te permita la memoria).

Pero sin necesidad de cambiar de lenguaje, si te mantienes en C, el número anterior nos cabría dentro de un unsigned long long int, pues este tipo de datos tiene 64 bits y ya hemos visto que el número en cuestión requería 62.

Tienes el problema de calcular la raíz cuadrada, ya que la función sqrt() opera con double. Y naturalmente tarde o temprano volverás a quedarte sin bits a medida que y crezca.

Podrías indagar en alguna de las muchas bibliotecas para precisión arbitraria para C, o cambiar de lenguaje. En python queda bastante sencillo (aunque no será tan rápido como una versión C).

Actualización

Pongo seguidamente una implementación en C que usa tipo de dato unsigned long long y por tanto puede llegar más lejos que la versión con double. En concreto, podría llegar hasta más de MAX=90000000, pues para este valor de y el resultado de 1141*y*y+1 aún cabe en 64 bits.

#include <stdio.h>

unsigned long long int isqrt(unsigned long long int n) {
    unsigned long long int x = n;
    unsigned long long int y = (x+1)/2;
    while (y<x) {
        x = y;
        y = (x+n/x)/2;
    }
    return x;
}


void main()
{
    unsigned long long int MAX=90000000;
    unsigned long long y, x2, x;
    for (y=0; y<MAX; y++) {
        x2 = 1141*y*y+1;
        x = isqrt(x2);
        if (x*x == x2){
            printf("\t%llu\n", y);
        }
    }
}

Observa que he tenido mucho cuidado de trabajar siempre con aritmética entera (he puesto 1141 y no 1141.0 por ejemplo). La función isqrt() calcula una aproximación entera de una raiz cuadrada, por el método de Newton. La raiz calculada será exacta si es entera (por ejemplo, isqrt(25) retornaría 5), pero se queda por debajo si no (por ejemplo, isqrt(26) retornaría también 5). Por eso luego elevo al cuadrado el valor retornado y lo comparo con lo que debía ser, para decidir si encontré o no una solución.

En mi máquina tarda un minuto en examinar los 90000000 casos, y lamentablemente no encuentra ninguna solución.

UPDATE. Como comparativa, la versión python del mismo programa, tarda 26 minutos en recorrer los 90000000 casos, cuando se ejecuta con el intérprete estándar.

Ejecutándo el mismo programa python con pypy, que es otra implementación mucho más rápida especialmente para cómputo numérico, el tiempo de ejecución desciende a 1m16'!! 😱

Este tiempo es prácticamente equivalente a C, de modo que podemos dejar que python siga buscando por encima de 90000000, pues en cuanto el dato no le quepa en 64 bits, automáticamente usará más, sin límite...

Esta es la versión python, que podrás comprobar es muy similar al código C anterior, pero sin declaraciones de tipos y sin llaves :-) (el operador // es el de división entera, ya que / hace división flotante, en python3)

def isqrt(n):
    x = n
    y = (x + 1) // 2
    while y < x:
        x = y
        y = (x + n // x) // 2
    return x


MAX=90000000
for y in range(MAX):
    x2 = 1141*y*y+1
    x = isqrt(x2)
    if x*x == x2:
        print("   %d" % y)
8
  • ok. Gracias por tu explicación. Veo que no se puede hacer en C entonces. Pensaba que era más potente, porque es el que estamos estudiando en físicas. Un saludo.
    – Vacendak
    Commented el 18 jul. 2018 a las 16:09
  • 1
    @Vacendak No he dicho yo que no se pueda en C. Sólo que no se puede con el tipo de datos double, y usar otro tipo mediante bibliotecas de terceros ya se va poniendo más complejo.
    – abulafia
    Commented el 18 jul. 2018 a las 16:11
  • Ya, pero eso queda fuera de mi nivel. Acabo de hacer otro programa que hace lo mismo sin usar sqrt , y usando el formato unsigned long long , como me propusiste en tu respuesta, pero es muy lento, porque tiene que realizar un montón de bucles...
    – Vacendak
    Commented el 18 jul. 2018 a las 16:32
  • 1
    Puedes utilizar el método de Newton para aproximar la raiz entera. De todas formas te anticipo. El problema no tiene solución entera. Lo he probado con python :-) Si quieres pongo el código.
    – abulafia
    Commented el 18 jul. 2018 a las 16:35
  • 1
    @Vacendak Mira mi actualización para la versión Python y el tiempo que requiere su ejecución.
    – abulafia
    Commented el 18 jul. 2018 a las 18:01

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.