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

En primer lugar, esta no es una pregunta de punto flotante para novatos. 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 terminantes no se pueden representar exactamente como números binarios de coma flotante.

Dicho esto, cada posible valor de coma flotante corresponde exactamente a un racional diádico (un número racional p/q donde q es una potencia de 2), que a su vez tiene un decimal exacto representatividad.

Mi pregunta es: ¿Cómo encuentra 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 manera única el valor de coma 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. He aquí un esquema:

  1. Convierte la mantisa a un entero decimal. Puedes puede hacer esto separando los bits para leer la mantisa directamente, o puede escribir un bucle de coma flotante desordenado que primero multiplique el valor por una potencia de dos para ponerlo en el rango 1
  2. Aplique el exponente multiplicando o dividiendo repetidamente por 2. Esta es una operación en la cadena de dígitos decimales que generó. Cada ~ 3 multiplicaciones añadirá un dígito adicional a la izquierda. Cada dividión agregará un dígito adicional a la derecha.

¿Es esto realmente lo mejor posible? Lo dudo, pero no soy un experto en punto flotante y no puedo encontrar una manera de hacer los cálculos de base-10 en la representación de punto flotante del número sin correr en una posibilidad de 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).

Author: R.., 2010-07-09
Source

10 answers

Esta pregunta tiene una parte burocrática y una parte algorítmica. Un número de coma flotante se almacena internamente como (2e × 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 más interesado en la parte algorítmica de la pregunta, a saber, convertir (2e × 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 conocía antes de hoy).

Es cierto que, a primera vista, toma O(e2) de su trabajo para escribir 2e en decimal, y más tiempo aún para la mantisa. Pero, gracias a la magia del algoritmo de multiplicación rápida Schonhage-Strassen , puede hacerlo en Õ(e) tiempo, donde la tilde significa "hasta log factores". Si ves Schonhage-Strassen como magia, entonces no es tan difícil pensar en qué hacer. Si e es par, puede recursivamente calcule 2e/2, y luego cuadrado usando multiplicación rápida. Por otro lado, si e es impar, puede calcular recursivamente 2e-1 y luego duplicarlo. Hay que tener cuidado para comprobar que hay una versión de Schonhage-Strassen en base 10. Aunque no está ampliamente documentado, 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 2k + b. Luego convierta recursivamente a y b a base 10, convierta 2^k a base 10, y haga otra multiplicación rápida para calcular m en base 10.

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

Si la pregunta es sobre números de coma flotante estándar de 64 bits, entonces es demasiado pequeña para el sofisticado algoritmo Schonhage-Strassen. En este rango puede guardar el trabajo con varios trucos. Un enfoque es almacene todos los 2048 valores de 2e en una tabla de búsqueda, y luego trabaje en la mantisa con multiplicación asimétrica (entre multiplicación larga y multiplicación corta). Otro truco es trabajar en base 10000 (o una potencia superior de 10, dependiendo de la arquitectura) en lugar de base 10. Pero, como R. señala en los comentarios, los números de coma flotante de 128 bits ya permiten exponentes lo suficientemente grandes como para cuestionar tanto las tablas de búsqueda como 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 después de eso una variación de Schonhage-Strassen es mejor no solo en teoría sino también en la práctica.

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

 34
Author: Greg Kuperberg,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2017-03-16 06:54:29

Veo que ya ha aceptado una respuesta, pero aquí hay un par de implementaciones de código abierto de esta conversión que puede que desee ver:

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

  2. La función ___printf_fp() en file /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 pida en un %f tipo printf (como he escrito acerca de 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/).

 16
Author: Rick Regan,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2013-08-23 12:00:48

Ha habido mucho trabajo en la impresión de números de coma flotante. El estándar de oro es imprimir un equivalente decimal de longitud mínima tal que cuando el equivalente decimal se lee de nuevo, se obtiene el mismo número de coma flotante con el que comenzó, no importa cuál sea el modo de redondeo durante la lectura. Puedes leer sobre el algoritmo en el excelente documento de Burger y Dybvig.

 5
Author: Norman Ramsey,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2010-07-09 19:52:30

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

De un vistazo rápido, no parece ser demasiado difícil portar a C, e incluso más fácil escribir en C++.

Tras una reflexión posterior, 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 coma flotante), y no estoy seguro de que pueda convertir de base 2 a base 10 en un tiempo mejor que log-cuadrado.

 3
Author: Gabe,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2010-07-09 18:17:56

Al no ser yo mismo un experto en coma flotante, me opondría al uso de una biblioteca de código abierto bien probada.

El GNU MPFR es bueno.

La biblioteca MPFR es una biblioteca C para punto flotante de precisión múltiple cálculos con redondeo correcto. El objetivo principal de MPFR es proporcionar un biblioteca para precisión múltiple cálculo de coma flotante que es tanto eficiente y tiene un bien definido semántica.

 2
Author: Byron Whitlock,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2010-07-09 17:55:21

Sprintf y funciones similares son por lo general, solo se especifica hasta un número de dígitos significativos a determinar el punto flotante original valor; no necesariamente imprimen la representación decimal exacta.

Puedes pedir dígitos más significativos que los predeterminados:

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

Imprime 0.1000000000000000055511151231257827021181583404541015625.

 1
Author: dan04,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2010-07-13 06:34:21

Si desea resultados más exactos, ¿por qué no usar matemáticas de punto fijo en su lugar? Las conversiones son rápidas. El error es conocido y puede ser solucionado. No es una respuesta exacta a tu pregunta, sino una idea diferente para ti.

 0
Author: Michael Dorgan,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2010-07-09 17:55:25

En la parte superior de mi cabeza, ¿por qué no dividir el exponente en una suma de exponentes binarios primero, entonces todas sus operaciones son sin pérdidas.

Es decir, 

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

Luego suma: 

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

Estoy pensando que rompe sería en el O(n) en el número de dígitos, el cambio es O(1), y la suma es O(n) dígitos...

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

Hágamelo saber - tengo curiosidad sobre esto, realmente me hizo pensar. :-)

 0
Author: eruciform,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2010-07-09 18:03:12

Hay 3 maneras

  1. Imprimir números en bin o hex

    Esta es la forma más precisa. Prefiero hex porque es más como 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 en los dígitos correspondientes

    Con esto quiero decir solo dígitos que pueden tener 1 en su número representado lo más cerca posible.

    num de entero los dígitos 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 difíciles 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 haces una tabla de dependencias n02 <-> n10 entonces verás que la constante 0.30102999566398119521373889472449 todavía está presente, pero al principio desde 8 bits porque less no puede representar 0.1 con precisión satisfactoria (0.85 - 1.15). debido a los exponentes negativos de base 2, la dependencia no es lineal, sino que se patrocina. Este código funciona para small bit count (<=52) pero en cuentas de bits más grandes puede haber error porque el patrón utilizado no encaja log10(2) o 1/log2(10) exactamente.

    Para cuentas de bits más grandes utilizo esto: 

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

    ¡Pero esa fórmula está alineada a 32 bits !!! y también más grande bit count anuncios error a ella

    P. S. análisis adicional 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 conduciría al número exacto de dígitos relevantes para cualquier bit contar.

    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 se olvide de redondear los dígitos cortados !!! Esto significa que si el dígito después del último dígito pertinente es >=5 en dec, el último dígito pertinente debe ser +1... y si ya es 9 que usted debe ir al dígito anterior y así sucesivamente...

  3. Imprimir valor exacto

    Para imprimir el valor exacto de número binario fraccionario {[41] } simplemente imprima dígitos fraccionarios n donde n es el número de bits fraccionarios porque el valor representado es la suma 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. Cosas de arriba (bala #2) es relevante para almacenar números decimales en float (o imprimir solo decimales relevantes).

    Potencias negativas 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 el estilo de valor exacto para 64 bits doubles: 

    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 negativo potencias de 10 impresas por solo dígitos decimales relevantes estilo (estoy más acostumbrado a esto) para 64bit doubles: 

    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 que ayude:)

 0
Author: Spektre,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2016-03-25 08:40:06

No. Lo más cercano que puedes llegar a eso es tirar los bytes.

 -2
Author: Steven Sudit,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2010-07-09 17:54:00