¿Cómo se imprime el valor EXACTO de un número de coma flotante?

En primer lugar, esta no es una pregunta para principiantes de coma flotante. Sé que los resultados de la aritmética de coma flotante (sin mencionar las funciones trascendentales) generalmente no se pueden representar exactamente, y que la mayoría de los decimales de terminación no se pueden representar exactamente como números binarios de coma flotante.

Dicho esto, cada posible valor de coma flotante corresponde exactamente a una razón diadica (un número racional p/q donde q es una potencia de 2), que a su vez tiene una representación decimal exacta.

Mi pregunta es: ¿cómo se puede encontrar esta representación decimal exacta de manera eficiente? sprintf y funciones similares generalmente solo se especifican hasta un número de dígitos significativos para determinar de forma exclusiva el valor del punto flotante original; no necesariamente imprimen la representación decimal exacta. Conozco un algoritmo que he usado, pero es muy lento, O(e^2) donde e es el exponente. Aquí hay un resumen:

  1. Convierta la mantisa en un entero decimal. Puede hacer esto separando los bits para leer la mantisa directamente, o puede escribir un bucle de punto flotante desordenado que primero multiplica el valor por una potencia de dos para ponerlo en el rango 1 <= x <10, luego tira de un dígito a la vez al convertir a int, restar y multiplicar por 10.
  2. Aplica el exponente multiplicando o dividiendo repetidamente por 2. Esta es una operación en la cadena de dígitos decimales que has generado. Cada ~ 3 multiplicaciones agregará un dígito extra a la izquierda. Cada división individual agregará un dígito extra a la derecha.

¿Es esto realmente lo mejor posible? Lo dudo, pero no soy un experto en coma flotante y no puedo encontrar una forma de hacer los cálculos de base 10 en la representación en coma flotante del número sin correr el riesgo de obtener resultados inexactos (multiplicar o dividir por cualquier cosa menos una potencia de 2 es una operación con pérdida en números de punto flotante a menos que sepa que tiene bits libres para trabajar).

Esta pregunta tiene una parte burocrática y una parte algorítmica. Un número de punto flotante se almacena internamente como (2 e × m), donde e es un exponente (en sí mismo en binario) y m es una mantisa. La parte burocrática de la pregunta es cómo acceder a estos datos, pero R. parece estar más interesado en la parte algorítmica de la pregunta, concretamente la conversión (2 e × m) a una fracción (a / b) en forma decimal. La respuesta a la pregunta burocrática en varios idiomas es “frexp” (que es un detalle interesante que no sabía antes de hoy).

Es cierto que a primera vista, se necesita O (e 2 ) para escribir 2 e en decimal y más tiempo para la mantisa. Pero, gracias a la magia del algoritmo de multiplicación rápida Schonhage-Strassen , puedes hacerlo en Õ (e) tiempo, donde la tilde significa “hasta factores de registro”. Si ves Schonhage-Strassen como magia, entonces no es tan difícil pensar en qué hacer. Si e es par, puede calcular recursivamente 2 e / 2 , y luego cuadrarlo usando una multiplicación rápida. Por otro lado, si e es impar, puede calcular recursivamente 2 e-1 y luego duplicarlo. Debes tener cuidado para verificar que haya una versión de Schonhage-Strassen en la base 10. Aunque no está ampliamente documentada, se puede hacer en cualquier base.

Convertir una mantisa muy larga de binaria a base 10 no es exactamente la misma pregunta, pero tiene una respuesta similar. Puedes dividir la mantisa en dos mitades, m = a 2 k + b. Luego convierta recursivamente ayb en la base 10, convierta 2 ^ k en base 10 y realice otra multiplicación rápida para calcular m en la base 10.

El resultado abstracto detrás de todo esto es que puedes convertir números enteros de una base a otra en Õ (N) tiempo.

Si la pregunta es acerca de los números de punto flotante estándar de 64 bits, entonces es demasiado pequeño para el sofisticado algoritmo Schonhage-Strassen. En este rango, en su lugar puede guardar el trabajo con varios trucos. Un enfoque es almacenar todos los 2048 valores de 2 e en una tabla de búsqueda, y luego trabajar en la mantisa con multiplicación asimétrica (entre multiplicación larga y multiplicación corta). Otro truco es trabajar en la base 10000 (o una potencia mayor de 10, dependiendo de la architecture) en lugar de la base 10. Pero, como señala R. en los comentarios, los números de punto flotante de 128 bits ya permiten exponentes suficientemente grandes para llamar cuestionar ambas tablas de búsqueda y la multiplicación larga estándar. Como cuestión práctica, la multiplicación larga es la más rápida hasta un puñado de dígitos, luego en un rango medio significativo se puede usar la multiplicación de Karatsuba o la multiplicación de Toom-Cook , y luego una variación de Schonhage-Strassen es mejor no solo en teoría pero también en la práctica.

En realidad, el gran paquete entero GMP ya tiene una conversión de Õ (N) radix de tiempo, así como una buena heurística para la elección del algoritmo de multiplicación. La única diferencia entre su solución y la mía es que en lugar de hacer una gran aritmética en la base 10, computan grandes potencias de 10 en la base 2. En esta solución, también necesitan división rápida, pero eso se puede obtener de la multiplicación rápida en cualquier de varias maneras.

Veo que ya ha aceptado una respuesta, pero aquí hay un par de implementaciones de código abierto de esta conversión que tal vez quiera consultar:

  1. La dtoa() de dtoa() David Gay en dtoa.c ( http://www.netlib.org/fp/dtoa.c ).

  2. La función ___printf_fp() en el archivo /stdio-common/printf_fp.c en glibc ( http://ftp.gnu.org/gnu/glibc/glibc-2.11.2.tar.gz , por ejemplo).

Ambos imprimirán tantos dígitos como pidas en un printf tipo %f (como he escrito aquí: http://www.exploringbinary.com/print-precision-of-dyadic-fractions-varies-by-language / y http://www.exploringbinary.com/print-precision-of-floating-point-integers-varies-too/ ).

Ha habido mucho trabajo para imprimir números de coma flotante. El estándar de oro consiste en imprimir un equivalente decimal de longitud mínima, de modo que cuando se vuelva a leer el equivalente decimal, obtenga el mismo número de punto flotante con el que comenzó, independientemente del modo de redondeo durante la lectura. Puede leer sobre el algoritmo en el excelente artículo de Burger y Dybvig .

Aunque es C # y su pregunta está etiquetada con C, Jon Skeet tiene un código para convertir un double en su representación exacta como una cadena: http://www.yoda.arachsys.com/csharp/DoubleConverter.cs

Desde un vistazo rápido, no parece ser demasiado difícil de portar a C, y aún más fácil de escribir en C ++.

Tras una reflexión adicional, parece que el algoritmo de Jon también es O (e ^ 2), ya que también gira sobre el exponente. Sin embargo, eso significa que el algoritmo es O (log (n) ^ 2) (donde n es el número de punto flotante), y no estoy seguro de que pueda convertir de la base 2 a la base 10 en un tiempo mejor que el log-cuadrado.

Al no ser experto en coma flotante, preferiría utilizar una biblioteca de código abierto bien probada.

El GNU MPFR es bueno.

La biblioteca MPFR es una biblioteca C para cálculos de punto flotante de precisión múltiple con redondeo correcto. El objective principal de MPFR es proporcionar una biblioteca para el cómputo de punto flotante de precisión múltiple que sea eficiente y tenga una semántica bien definida.

sprintf y funciones similares generalmente solo se especifican hasta un número de dígitos significativos para determinar de forma exclusiva el valor del punto flotante original; no necesariamente imprimen la representación decimal exacta.

Puede solicitar dígitos más significativos que los predeterminados:

 printf("%.100g\n", 0.1); 

imprime 0.1000000000000000055511151231257827021181583404541015625 .

Si desea resultados más exactos, ¿por qué no utilizar matemáticas de punto fijo en su lugar? Las conversiones son rápidas. Se conoce un error y se puede solucionar. No es una respuesta exacta a su pregunta, sino una idea diferente para usted.

Por la parte superior de mi cabeza, ¿por qué no dividir el exponente en una sum de exponentes binarios primero, y luego todas tus operaciones no tienen pérdida?

Es decir

 10^2 = 2^6 + 2^5 + 2^2 

Entonces sum:

 mantissa<<6 + mantissa<<5 + mantissa<<2 

Estoy pensando que dividirlo sería en O (n) en el número de dígitos, el desplazamiento es O (1) y la sum es O (n) dígitos ...

Tendría que tener una clase entera lo suficientemente grande como para almacenar los resultados, por supuesto ...

Avíseme - Tengo curiosidad sobre esto, realmente me hizo pensar. 🙂

Hay 3 formas

  1. imprimir números en el bin o en el hex

    Esta es la forma más precisa. Prefiero el hex porque es más como la base 10 para leer / sentir como por ejemplo F.8h = 15.5 no hay pérdida de precisión aquí.

  2. impresión en dec pero solo los dígitos relevantes

    Con esto me refiero solo a los dígitos que pueden tener 1 en su número representado lo más cerca posible.

    número de dígitos enteros son fáciles y precisos (sin pérdida de precisión):

     // n10 - base 10 integer digits // n2 - base 2 integer digits n10=log10(2^n2) n10=log2(2^n2)/log2(10) n10=n2/log2(10) n10=ceil(n2*0.30102999566398119521373889472449) // if fist digit is 0 and n10 > 1 then n10-- 

    num de dígitos fraccionarios son más complicados y empíricamente encontré esto:

     // n10 - base 10 fract. digits // n2 - base 2 fract. digits >= 8 n10=0; if (n02==8) n10=1; else if (n02==9) n10=2; else if (n02> 9) { n10=((n02-9)%10); if (n10>=6) n10=2; else if (n10>=1) n10=1; n10+=2+(((n02-9)/10)*3); } 

    si n02 <-> n10 una tabla de dependencia n02 <-> n10 , verá que la constante 0.30102999566398119521373889472449 todavía está presente, pero al inicio de 8 bits porque menos no puede representar 0.1 con una precisión satisfactoria ( 0.85 - 1.15 ). debido a los exponentes negativos de la base 2 la dependencia no es lineal, sino patrones. Este código funciona para el conteo de bits pequeños ( <=52 ) pero en conteos de bits más grandes puede haber un error porque el patrón utilizado no se ajusta exactamente a log10(2) o 1/log2(10) .

    para cuentas de bits más grandes uso esto:

     n10=7.810+(9.6366363636363636363636*((n02>>5)-1.0)); 

    pero esa fórmula está alineada 32bit !!! y también un mayor error de anuncios de conteo de bits

    PS más análisis de la representación binaria de números decádicos

     0.1 0.01 0.001 0.0001 ... 

    puede revelar la repetición exacta del patrón que llevaría a la cantidad exacta de dígitos relevantes para cualquier cantidad de bits.

    para mayor claridad:

     8 bin digits -> 1 dec digits 9 bin digits -> 2 dec digits 10 bin digits -> 3 dec digits 11 bin digits -> 3 dec digits 12 bin digits -> 3 dec digits 13 bin digits -> 3 dec digits 14 bin digits -> 3 dec digits 15 bin digits -> 4 dec digits 16 bin digits -> 4 dec digits 17 bin digits -> 4 dec digits 18 bin digits -> 4 dec digits 19 bin digits -> 5 dec digits 20 bin digits -> 6 dec digits 21 bin digits -> 6 dec digits 22 bin digits -> 6 dec digits 23 bin digits -> 6 dec digits 24 bin digits -> 6 dec digits 25 bin digits -> 7 dec digits 26 bin digits -> 7 dec digits 27 bin digits -> 7 dec digits 28 bin digits -> 7 dec digits 29 bin digits -> 8 dec digits 30 bin digits -> 9 dec digits 31 bin digits -> 9 dec digits 32 bin digits -> 9 dec digits 33 bin digits -> 9 dec digits 34 bin digits -> 9 dec digits 35 bin digits -> 10 dec digits 36 bin digits -> 10 dec digits 37 bin digits -> 10 dec digits 38 bin digits -> 10 dec digits 39 bin digits -> 11 dec digits 40 bin digits -> 12 dec digits 41 bin digits -> 12 dec digits 42 bin digits -> 12 dec digits 43 bin digits -> 12 dec digits 44 bin digits -> 12 dec digits 45 bin digits -> 13 dec digits 46 bin digits -> 13 dec digits 47 bin digits -> 13 dec digits 48 bin digits -> 13 dec digits 49 bin digits -> 14 dec digits 50 bin digits -> 15 dec digits 51 bin digits -> 15 dec digits 52 bin digits -> 15 dec digits 53 bin digits -> 15 dec digits 54 bin digits -> 15 dec digits 55 bin digits -> 16 dec digits 56 bin digits -> 16 dec digits 57 bin digits -> 16 dec digits 58 bin digits -> 16 dec digits 59 bin digits -> 17 dec digits 60 bin digits -> 18 dec digits 61 bin digits -> 18 dec digits 62 bin digits -> 18 dec digits 63 bin digits -> 18 dec digits 64 bin digits -> 18 dec digits 

    ¡Y por fin no te olvides de redondear los dígitos cortados! Eso significa que si el dígito después del último dígito relevante es >=5 en decimales que el último dígito relevante debe ser +1 ... y si ya es 9 , debe ir al dígito anterior y así sucesivamente ...

  3. imprimir el valor exacto

    Para imprimir el valor exacto del número binario fraccionario simplemente imprima n dígitos fraccionarios donde n es el número de bits fraccionarios porque el valor representado es la sum de estos valores por lo que el número de decimales fraccionarios no puede ser mayor que el num de dígitos fraccionarios de LSB . La información anterior (viñeta n.º 2 ) es relevante para almacenar números decimales en float (o para imprimir solo decimales relevantes).

    poderes negativos de dos valores exactos ...

     2^- 1 = 0.5 2^- 2 = 0.25 2^- 3 = 0.125 2^- 4 = 0.0625 2^- 5 = 0.03125 2^- 6 = 0.015625 2^- 7 = 0.0078125 2^- 8 = 0.00390625 2^- 9 = 0.001953125 2^-10 = 0.0009765625 2^-11 = 0.00048828125 2^-12 = 0.000244140625 2^-13 = 0.0001220703125 2^-14 = 0.00006103515625 2^-15 = 0.000030517578125 2^-16 = 0.0000152587890625 2^-17 = 0.00000762939453125 2^-18 = 0.000003814697265625 2^-19 = 0.0000019073486328125 2^-20 = 0.00000095367431640625 

    ahora potencias negativas de 10 impresas por estilo de valor exacto para doubles 64 bits:

     10^+ -1 = 0.1000000000000000055511151231257827021181583404541015625 = 0.0001100110011001100110011001100110011001100110011001101b 10^+ -2 = 0.01000000000000000020816681711721685132943093776702880859375 = 0.00000010100011110101110000101000111101011100001010001111011b 10^+ -3 = 0.001000000000000000020816681711721685132943093776702880859375 = 0.000000000100000110001001001101110100101111000110101001111111b 10^+ -4 = 0.000100000000000000004792173602385929598312941379845142364501953125 = 0.000000000000011010001101101110001011101011000111000100001100101101b 10^+ -5 = 0.000010000000000000000818030539140313095458623138256371021270751953125 = 0.000000000000000010100111110001011010110001000111000110110100011110001b 10^+ -6 = 0.000000999999999999999954748111825886258685613938723690807819366455078125 = 0.000000000000000000010000110001101111011110100000101101011110110110001101b 10^+ -7 = 0.0000000999999999999999954748111825886258685613938723690807819366455078125 = 0.0000000000000000000000011010110101111111001010011010101111001010111101001b 10^+ -8 = 0.000000010000000000000000209225608301284726753266340892878361046314239501953125 = 0.000000000000000000000000001010101111001100011101110001000110000100011000011101b 10^+ -9 = 0.0000000010000000000000000622815914577798564188970686927859787829220294952392578125 = 0.0000000000000000000000000000010001001011100000101111101000001001101101011010010101b 10^+-10 = 0.00000000010000000000000000364321973154977415791655470655996396089904010295867919921875 = 0.00000000000000000000000000000000011011011111001101111111011001110101111011110110111011b 10^+-11 = 0.00000000000999999999999999939496969281939810930172340963650867706746794283390045166015625 = 0.00000000000000000000000000000000000010101111111010111111111100001011110010110010010010101b 10^+-12 = 0.00000000000099999999999999997988664762925561536725284350612952266601496376097202301025390625 = 0.00000000000000000000000000000000000000010001100101111001100110000001001011011110101000010001b 10^+-13 = 0.00000000000010000000000000000303737455634003709136034716842278413651001756079494953155517578125 = 0.00000000000000000000000000000000000000000001110000100101110000100110100001001001011101101000001b 10^+-14 = 0.000000000000009999999999999999988193093545598986971343290729163921781719182035885751247406005859375 = 0.000000000000000000000000000000000000000000000010110100001001001101110000110101000010010101110011011b 10^+-15 = 0.00000000000000100000000000000007770539987666107923830718560119501514549256171449087560176849365234375 = 0.00000000000000000000000000000000000000000000000001001000000011101011111001111011100111010101100001011b 10^+-16 = 0.00000000000000009999999999999999790977867240346035618411149408467364363417573258630000054836273193359375 = 0.00000000000000000000000000000000000000000000000000000111001101001010110010100101111101100010001001101111b 10^+-17 = 0.0000000000000000100000000000000007154242405462192450852805618492324772617063644020163337700068950653076171875 = 0.0000000000000000000000000000000000000000000000000000000010111000011101111010101000110010001101101010010010111b 10^+-18 = 0.00000000000000000100000000000000007154242405462192450852805618492324772617063644020163337700068950653076171875 = 0.00000000000000000000000000000000000000000000000000000000000100100111001001011101110100011101001001000011101011b 10^+-19 = 0.000000000000000000099999999999999997524592683526013185572915905567688179926555402943222361500374972820281982421875 = 0.000000000000000000000000000000000000000000000000000000000000000111011000001111001001010011111011011011010010101011b 10^+-20 = 0.00000000000000000000999999999999999945153271454209571651729503702787392447107715776066783064379706047475337982177734375 = 0.00000000000000000000000000000000000000000000000000000000000000000010111100111001010000100001100100100100100001000100011b 

    ahora potencias negativas de 10 impresas por dígitos decimales relevantes solo estilo (estoy más acostumbrado a esto) para doubles 64 bits:

     10^+ -1 = 0.1 10^+ -2 = 0.01 10^+ -3 = 0.001 10^+ -4 = 0.0001 10^+ -5 = 0.00001 10^+ -6 = 0.000001 10^+ -7 = 0.0000001 10^+ -8 = 0.00000001 10^+ -9 = 0.000000001 10^+-10 = 0.0000000001 10^+-11 = 0.00000000001 10^+-12 = 0.000000000001 10^+-13 = 0.0000000000001 10^+-14 = 0.00000000000001 10^+-15 = 0.000000000000001 10^+-16 = 0.0000000000000001 10^+-17 = 0.00000000000000001 10^+-18 = 0.000000000000000001 10^+-19 = 0.0000000000000000001 10^+-20 = 0.00000000000000000001 

    Espero eso ayude 🙂

Tu no Lo más cerca que puedes llegar a eso es tirar los bytes.