float error arithmetic c floating-point floating-point-precision pow

error - pow() parece estar fuera por uno aquí



floating point representation (4)

Que está pasando aqui:

#include <stdio.h> #include <math.h> int main(void) { printf("17^12 = %lf/n", pow(17, 12)); printf("17^13 = %lf/n", pow(17, 13)); printf("17^14 = %lf/n", pow(17, 14)); }

Obtengo esta salida:

17^12 = 582622237229761.000000 17^13 = 9904578032905936.000000 17^14 = 168377826559400928.000000

13 y 14 no coinciden con wolfram alpa cf:

12: 582622237229761.000000 582622237229761 13: 9904578032905936.000000 9904578032905937 14: 168377826559400928.000000 168377826559400929

Además, no está mal por alguna fracción extraña, ¡está mal por exactamente una!

Si me toca llegar a los límites de lo que pow() puede hacer por mí, ¿hay alguna alternativa que pueda calcular esto? Necesito una función que pueda calcular x^y , donde x^y siempre es menor que ULLONG_MAX.


Los números que obtienes son demasiado grandes para ser representados con un double precisión. Un número de punto flotante de doble precisión tiene esencialmente 53 dígitos binarios significativos y puede representar todos los enteros hasta 2^53 o 9,007,199,254,740,992.

Para números más altos, los últimos dígitos se truncan y el resultado de su cálculo se redondea al siguiente número que se puede representar como un double . Para 17^13 , que está solo ligeramente por encima del límite, este es el número par más cercano. Para números mayores que 2^54 este es el número más cercano que es divisible por cuatro, y así sucesivamente.


Si sus argumentos de entrada son enteros no negativos, entonces puede implementar su propio poder.

Recursivamente:

unsigned long long pow(unsigned long long x,unsigned int y) { if (y == 0) return 1; if (y == 1) return x; return pow(x,y/2)*pow(x,y-y/2); }

Iterativamente

unsigned long long pow(unsigned long long x,unsigned int y) { unsigned long long res = 1; while (y--) res *= x; return res; }

Eficientemente:

unsigned long long pow(unsigned long long x,unsigned int y) { unsigned long long res = 1; while (y > 0) { if (y & 1) res *= x; y >>= 1; x *= x; } return res; }


Una pequeña adición a otras buenas respuestas: bajo la arquitectura x86, generalmente está disponible el formato extendido x87 de 80 bits , que es compatible con la mayoría de los compiladores de C a través del tipo long double . Este formato permite operar con números enteros hasta 2^64 sin espacios.

Hay un análogo de pow() en <math.h> que está diseñado para operar con números long double - powl() . También debe tenerse en cuenta que el especificador de formato para los valores long double no es para los double - %Lf . Así que el programa correcto que usa el tipo long double ve así:

#include <stdio.h> #include <math.h> int main(void) { printf("17^12 = %Lf/n", powl(17, 12)); printf("17^13 = %Lf/n", powl(17, 13)); printf("17^14 = %Lf/n", powl(17, 14)); }

Como señaló Stephen Canon en los comentarios, no hay garantía de que este programa dé resultados exactos.


pow funciona con números double . Estos representan números de la forma s * 2 ^ e donde s es un entero de 53 bits. Por lo tanto, double puede almacenar todos los enteros por debajo de 2 ^ 53, pero solo algunos enteros por encima de 2 ^ 53. En particular, solo puede representar números pares> 2 ^ 53, ya que para e> 0 el valor es siempre un múltiplo de 2.

17 ^ 13 necesita 54 bits para representar exactamente, por lo que e se establece en 1 y, por tanto, el valor calculado se convierte en un número par. El valor correcto es impar, por lo que no es sorprendente que esté desactivado en uno. Asimismo, 17 ^ 14 toma 58 bits para representar. El hecho de que también esté desactivado por uno es una coincidencia afortunada (siempre que no aplique demasiada teoría de números), resulta que está fuera de un múltiplo de 32 , que es la granularidad con la que double números de esa magnitud son redondeados

Para la exponenciación exacta de enteros, debe usar enteros hasta el final. Escribe tu propia rutina de double exponencia. Use la exponenciación cuadrada si y puede ser grande, pero supongo que siempre es menor que 64, lo que hace que este problema sea discutible.