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