¿Por qué una conversión de ida y vuelta a través de una cadena no es segura para un doble?

185

Recientemente tuve que serializar un doble en texto y luego recuperarlo. El valor parece no ser equivalente:

double d1 = 0.84551240822557006;
string s = d1.ToString("R");
double d2 = double.Parse(s);
bool s1 = d1 == d2;
// -> s1 is False

Pero según MSDN: cadenas de formato numérico estándar , se supone que la opción "R" garantiza la seguridad de ida y vuelta.

El especificador de formato de ida y vuelta ("R") se utiliza para garantizar que un valor numérico que se convierte en una cadena se analice nuevamente en el mismo valor numérico

¿Por qué pasó esto?

Philip Ding
fuente
66
Depuré mi VS y está volviendo cierto aquí
Neel
19
Lo he reproducido devolviendo falso. Pregunta muy interesante
Jon Skeet
40
.net 4.0 x86 - verdadero, .net 4.0 x64 - falso
Ulugbek Umirov
25
Felicitaciones por encontrar un error tan impresionante en .net.
Aron
14
@Casperah El viaje de ida y vuelta está específicamente destinado a evitar inconsistencias de coma flotante
Gusdor

Respuestas:

178

Encontré el error.

.NET hace lo siguiente en clr\src\vm\comnumber.cpp:

DoubleToNumber(value, DOUBLE_PRECISION, &number);

if (number.scale == (int) SCALE_NAN) {
    gc.refRetVal = gc.numfmt->sNaN;
    goto lExit;
}

if (number.scale == SCALE_INF) {
    gc.refRetVal = (number.sign? gc.numfmt->sNegativeInfinity: gc.numfmt->sPositiveInfinity);
    goto lExit;
}

NumberToDouble(&number, &dTest);

if (dTest == value) {
    gc.refRetVal = NumberToString(&number, 'G', DOUBLE_PRECISION, gc.numfmt);
    goto lExit;
}

DoubleToNumber(value, 17, &number);

DoubleToNumberes bastante simple: solo llama _ecvt, que está en el tiempo de ejecución de C:

void DoubleToNumber(double value, int precision, NUMBER* number)
{
    WRAPPER_CONTRACT
    _ASSERTE(number != NULL);

    number->precision = precision;
    if (((FPDOUBLE*)&value)->exp == 0x7FF) {
        number->scale = (((FPDOUBLE*)&value)->mantLo || ((FPDOUBLE*)&value)->mantHi) ? SCALE_NAN: SCALE_INF;
        number->sign = ((FPDOUBLE*)&value)->sign;
        number->digits[0] = 0;
    }
    else {
        char* src = _ecvt(value, precision, &number->scale, &number->sign);
        wchar* dst = number->digits;
        if (*src != '0') {
            while (*src) *dst++ = *src++;
        }
        *dst = 0;
    }
}

Resulta que _ecvtdevuelve la cadena 845512408225570.

¿Nota el cero final? ¡Resulta que hace toda la diferencia!
Cuando el cero está presente, el resultado vuelve a analizar0.84551240822557006, que es sunúmero original , por lo que se compara igual y, por lo tanto, solo se devuelven 15 dígitos.

Sin embargo, si trunco ​​la cadena en ese cero 84551240822557, entonces vuelvo 0.84551240822556994, que no es su número original, y por lo tanto devolvería 17 dígitos.

Prueba: ejecute el siguiente código de 64 bits (la mayoría de los cuales extraje de Microsoft Shared Source CLI 2.0) en su depurador y examine val final de main:

#include <stdlib.h>
#include <string.h>
#include <math.h>

#define min(a, b) (((a) < (b)) ? (a) : (b))

struct NUMBER {
    int precision;
    int scale;
    int sign;
    wchar_t digits[20 + 1];
    NUMBER() : precision(0), scale(0), sign(0) {}
};


#define I64(x) x##LL
static const unsigned long long rgval64Power10[] = {
    // powers of 10
    /*1*/ I64(0xa000000000000000),
    /*2*/ I64(0xc800000000000000),
    /*3*/ I64(0xfa00000000000000),
    /*4*/ I64(0x9c40000000000000),
    /*5*/ I64(0xc350000000000000),
    /*6*/ I64(0xf424000000000000),
    /*7*/ I64(0x9896800000000000),
    /*8*/ I64(0xbebc200000000000),
    /*9*/ I64(0xee6b280000000000),
    /*10*/ I64(0x9502f90000000000),
    /*11*/ I64(0xba43b74000000000),
    /*12*/ I64(0xe8d4a51000000000),
    /*13*/ I64(0x9184e72a00000000),
    /*14*/ I64(0xb5e620f480000000),
    /*15*/ I64(0xe35fa931a0000000),

    // powers of 0.1
    /*1*/ I64(0xcccccccccccccccd),
    /*2*/ I64(0xa3d70a3d70a3d70b),
    /*3*/ I64(0x83126e978d4fdf3c),
    /*4*/ I64(0xd1b71758e219652e),
    /*5*/ I64(0xa7c5ac471b478425),
    /*6*/ I64(0x8637bd05af6c69b7),
    /*7*/ I64(0xd6bf94d5e57a42be),
    /*8*/ I64(0xabcc77118461ceff),
    /*9*/ I64(0x89705f4136b4a599),
    /*10*/ I64(0xdbe6fecebdedd5c2),
    /*11*/ I64(0xafebff0bcb24ab02),
    /*12*/ I64(0x8cbccc096f5088cf),
    /*13*/ I64(0xe12e13424bb40e18),
    /*14*/ I64(0xb424dc35095cd813),
    /*15*/ I64(0x901d7cf73ab0acdc),
};

static const signed char rgexp64Power10[] = {
    // exponents for both powers of 10 and 0.1
    /*1*/ 4,
    /*2*/ 7,
    /*3*/ 10,
    /*4*/ 14,
    /*5*/ 17,
    /*6*/ 20,
    /*7*/ 24,
    /*8*/ 27,
    /*9*/ 30,
    /*10*/ 34,
    /*11*/ 37,
    /*12*/ 40,
    /*13*/ 44,
    /*14*/ 47,
    /*15*/ 50,
};

static const unsigned long long rgval64Power10By16[] = {
    // powers of 10^16
    /*1*/ I64(0x8e1bc9bf04000000),
    /*2*/ I64(0x9dc5ada82b70b59e),
    /*3*/ I64(0xaf298d050e4395d6),
    /*4*/ I64(0xc2781f49ffcfa6d4),
    /*5*/ I64(0xd7e77a8f87daf7fa),
    /*6*/ I64(0xefb3ab16c59b14a0),
    /*7*/ I64(0x850fadc09923329c),
    /*8*/ I64(0x93ba47c980e98cde),
    /*9*/ I64(0xa402b9c5a8d3a6e6),
    /*10*/ I64(0xb616a12b7fe617a8),
    /*11*/ I64(0xca28a291859bbf90),
    /*12*/ I64(0xe070f78d39275566),
    /*13*/ I64(0xf92e0c3537826140),
    /*14*/ I64(0x8a5296ffe33cc92c),
    /*15*/ I64(0x9991a6f3d6bf1762),
    /*16*/ I64(0xaa7eebfb9df9de8a),
    /*17*/ I64(0xbd49d14aa79dbc7e),
    /*18*/ I64(0xd226fc195c6a2f88),
    /*19*/ I64(0xe950df20247c83f8),
    /*20*/ I64(0x81842f29f2cce373),
    /*21*/ I64(0x8fcac257558ee4e2),

    // powers of 0.1^16
    /*1*/ I64(0xe69594bec44de160),
    /*2*/ I64(0xcfb11ead453994c3),
    /*3*/ I64(0xbb127c53b17ec165),
    /*4*/ I64(0xa87fea27a539e9b3),
    /*5*/ I64(0x97c560ba6b0919b5),
    /*6*/ I64(0x88b402f7fd7553ab),
    /*7*/ I64(0xf64335bcf065d3a0),
    /*8*/ I64(0xddd0467c64bce4c4),
    /*9*/ I64(0xc7caba6e7c5382ed),
    /*10*/ I64(0xb3f4e093db73a0b7),
    /*11*/ I64(0xa21727db38cb0053),
    /*12*/ I64(0x91ff83775423cc29),
    /*13*/ I64(0x8380dea93da4bc82),
    /*14*/ I64(0xece53cec4a314f00),
    /*15*/ I64(0xd5605fcdcf32e217),
    /*16*/ I64(0xc0314325637a1978),
    /*17*/ I64(0xad1c8eab5ee43ba2),
    /*18*/ I64(0x9becce62836ac5b0),
    /*19*/ I64(0x8c71dcd9ba0b495c),
    /*20*/ I64(0xfd00b89747823938),
    /*21*/ I64(0xe3e27a444d8d991a),
};

static const signed short rgexp64Power10By16[] = {
    // exponents for both powers of 10^16 and 0.1^16
    /*1*/ 54,
    /*2*/ 107,
    /*3*/ 160,
    /*4*/ 213,
    /*5*/ 266,
    /*6*/ 319,
    /*7*/ 373,
    /*8*/ 426,
    /*9*/ 479,
    /*10*/ 532,
    /*11*/ 585,
    /*12*/ 638,
    /*13*/ 691,
    /*14*/ 745,
    /*15*/ 798,
    /*16*/ 851,
    /*17*/ 904,
    /*18*/ 957,
    /*19*/ 1010,
    /*20*/ 1064,
    /*21*/ 1117,
};

static unsigned DigitsToInt(wchar_t* p, int count)
{
    wchar_t* end = p + count;
    unsigned res = *p - '0';
    for ( p = p + 1; p < end; p++) {
        res = 10 * res + *p - '0';
    }
    return res;
}
#define Mul32x32To64(a, b) ((unsigned long long)((unsigned long)(a)) * (unsigned long long)((unsigned long)(b)))

static unsigned long long Mul64Lossy(unsigned long long a, unsigned long long b, int* pexp)
{
    // it's ok to losse some precision here - Mul64 will be called
    // at most twice during the conversion, so the error won't propagate
    // to any of the 53 significant bits of the result
    unsigned long long val = Mul32x32To64(a >> 32, b >> 32) +
        (Mul32x32To64(a >> 32, b) >> 32) +
        (Mul32x32To64(a, b >> 32) >> 32);

    // normalize
    if ((val & I64(0x8000000000000000)) == 0) { val <<= 1; *pexp -= 1; }

    return val;
}

void NumberToDouble(NUMBER* number, double* value)
{
    unsigned long long val;
    int exp;
    wchar_t* src = number->digits;
    int remaining;
    int total;
    int count;
    int scale;
    int absscale;
    int index;

    total = (int)wcslen(src);
    remaining = total;

    // skip the leading zeros
    while (*src == '0') {
        remaining--;
        src++;
    }

    if (remaining == 0) {
        *value = 0;
        goto done;
    }

    count = min(remaining, 9);
    remaining -= count;
    val = DigitsToInt(src, count);

    if (remaining > 0) {
        count = min(remaining, 9);
        remaining -= count;

        // get the denormalized power of 10
        unsigned long mult = (unsigned long)(rgval64Power10[count-1] >> (64 - rgexp64Power10[count-1]));
        val = Mul32x32To64(val, mult) + DigitsToInt(src+9, count);
    }

    scale = number->scale - (total - remaining);
    absscale = abs(scale);
    if (absscale >= 22 * 16) {
        // overflow / underflow
        *(unsigned long long*)value = (scale > 0) ? I64(0x7FF0000000000000) : 0;
        goto done;
    }

    exp = 64;

    // normalize the mantisa
    if ((val & I64(0xFFFFFFFF00000000)) == 0) { val <<= 32; exp -= 32; }
    if ((val & I64(0xFFFF000000000000)) == 0) { val <<= 16; exp -= 16; }
    if ((val & I64(0xFF00000000000000)) == 0) { val <<= 8; exp -= 8; }
    if ((val & I64(0xF000000000000000)) == 0) { val <<= 4; exp -= 4; }
    if ((val & I64(0xC000000000000000)) == 0) { val <<= 2; exp -= 2; }
    if ((val & I64(0x8000000000000000)) == 0) { val <<= 1; exp -= 1; }

    index = absscale & 15;
    if (index) {
        int multexp = rgexp64Power10[index-1];
        // the exponents are shared between the inverted and regular table
        exp += (scale < 0) ? (-multexp + 1) : multexp;

        unsigned long long multval = rgval64Power10[index + ((scale < 0) ? 15 : 0) - 1];
        val = Mul64Lossy(val, multval, &exp);
    }

    index = absscale >> 4;
    if (index) {
        int multexp = rgexp64Power10By16[index-1];
        // the exponents are shared between the inverted and regular table
        exp += (scale < 0) ? (-multexp + 1) : multexp;

        unsigned long long multval = rgval64Power10By16[index + ((scale < 0) ? 21 : 0) - 1];
        val = Mul64Lossy(val, multval, &exp);
    }

    // round & scale down
    if ((unsigned long)val & (1 << 10))
    {
        // IEEE round to even
        unsigned long long tmp = val + ((1 << 10) - 1) + (((unsigned long)val >> 11) & 1);
        if (tmp < val) {
            // overflow
            tmp = (tmp >> 1) | I64(0x8000000000000000);
            exp += 1;
        }
        val = tmp;
    }
    val >>= 11;

    exp += 0x3FE;

    if (exp <= 0) {
        if (exp <= -52) {
            // underflow
            val = 0;
        }
        else {
            // denormalized
            val >>= (-exp+1);
        }
    }
    else
        if (exp >= 0x7FF) {
            // overflow
            val = I64(0x7FF0000000000000);
        }
        else {
            val = ((unsigned long long)exp << 52) + (val & I64(0x000FFFFFFFFFFFFF));
        }

        *(unsigned long long*)value = val;

done:
        if (number->sign) *(unsigned long long*)value |= I64(0x8000000000000000);
}

int main()
{
    NUMBER number;
    number.precision = 15;
    double v = 0.84551240822557006;
    char *src = _ecvt(v, number.precision, &number.scale, &number.sign);
    int truncate = 0;  // change to 1 if you want to truncate
    if (truncate)
    {
        while (*src && src[strlen(src) - 1] == '0')
        {
            src[strlen(src) - 1] = 0;
        }
    }
    wchar_t* dst = number.digits;
    if (*src != '0') {
        while (*src) *dst++ = *src++;
    }
    *dst++ = 0;
    NumberToDouble(&number, &v);
    return 0;
}
usuario541686
fuente
44
Buena explicación +1. Este código es de shared-source-cli-2.0 ¿verdad? Este es el único pensamiento que encontré.
Soner Gönül
10
Debo decir que es bastante patético. Las cadenas que son matemáticamente iguales (como una con un cero al final, o digamos 2.1e-1 vs. 0.21) siempre deberían dar resultados idénticos, y las cadenas que están matemáticamente ordenadas deberían dar resultados consistentes con el orden.
gnasher729
44
@ MrLister: ¿Por qué no debería "2.1E-1 ser igual a 0.21 así"?
user541686
9
@ gnasher729: Estoy un poco de acuerdo con "2.1e-1" y "0.21" ... pero una cadena con un cero al final no es exactamente igual a uno sin: en el primero, el cero es un dígito significativo y agrega precisión.
cHao
44
@ cHao: Er ... agrega precisión, pero eso solo afecta la forma en que decides redondear la respuesta final si las señales te importan, no cómo la computadora debe calcular la respuesta final en primer lugar. El trabajo de la computadora es calcular todo con la máxima precisión, independientemente de las precisiones de medición reales de los números; Es el problema del programador si quiere redondear el resultado final.
user541686
107

Me parece que esto es simplemente un error. Sus expectativas son completamente razonables. Lo reproduje usando .NET 4.5.1 (x64), ejecutando la siguiente aplicación de consola que usa mi DoubleConverterclase. DoubleConverter.ToExactStringmuestra el valor exacto representado por a double:

using System;

class Test
{
    static void Main()
    {
        double d1 = 0.84551240822557006;
        string s = d1.ToString("r");
        double d2 = double.Parse(s);
        Console.WriteLine(s);
        Console.WriteLine(DoubleConverter.ToExactString(d1));
        Console.WriteLine(DoubleConverter.ToExactString(d2));
        Console.WriteLine(d1 == d2);
    }
}

Resultados en .NET:

0.84551240822557
0.845512408225570055719799711368978023529052734375
0.84551240822556994469749724885332398116588592529296875
False

Resultados en Mono 3.3.0:

0.84551240822557006
0.845512408225570055719799711368978023529052734375
0.845512408225570055719799711368978023529052734375
True

Si especifica manualmente la cadena de Mono (que contiene el "006" en el extremo), .NET lo analizará nuevamente al valor original. Parece que el problema está en el ToString("R")manejo y no en el análisis.

Como se señaló en otros comentarios, parece que esto es específico para ejecutarse bajo el CLR x64. Si compila y ejecuta el código anterior dirigido a x86, está bien:

csc /platform:x86 Test.cs DoubleConverter.cs

... obtienes los mismos resultados que con Mono. Sería interesante saber si el error aparece en RyuJIT: no tengo eso instalado en este momento. En particular, puedo imaginar que esto posiblemente sea ​​un error JIT, o es muy posible que haya implementaciones completamente diferentes de los componentes internos double.ToStringbasados ​​en la arquitectura.

Le sugiero que presente un error en http://connect.microsoft.com

Jon Skeet
fuente
1
Entonces Jon? Para confirmar, ¿es esto un error en el JITer, que incluye el ToString()? Cuando intenté reemplazar el valor codificado con rand.NextDouble()y no hubo ningún problema.
Aron
1
Sí, definitivamente está en la ToString("R")conversión. Intente ToString("G32")y observe que imprime el valor correcto.
user541686
1
@Aron: No puedo decir si es un error en el JITter o en una implementación específica de x64 del BCL. Sin embargo, dudo mucho que sea tan simple como en línea. Las pruebas con valores aleatorios realmente no ayudan mucho, en mi opinión ... No estoy seguro de lo que esperas que demuestre.
Jon Skeet
2
Creo que lo que está sucediendo es que el formato de "ida y vuelta" está produciendo un valor que es 0.498ulp más grande de lo que debería ser, y la lógica de análisis a veces lo redondea erróneamente en la última fracción minúscula de un ulp. No estoy seguro de qué código culpo más, ya que creo que un formato de "ida y vuelta" debería generar un valor numérico que esté dentro de un cuarto de ULP de ser numéricamente correcto; la lógica de análisis que produce un valor dentro de 0.75ulp de lo especificado es mucho más fácil que la lógica que debe producir un resultado dentro de 0.502ulp de lo especificado.
supercat
1
¿El sitio web de Jon Skeet está caído? Me parece muy poco probable que esté ... perdiendo toda la fe aquí.
Patrick M
2

Recientemente, estoy tratando de resolver este problema . Como se señaló a través del código , el double.ToString ("R") tiene la siguiente lógica:

  1. Intenta convertir el doble a cadena con una precisión de 15.
  2. Convierta la cadena de nuevo al doble y compárela con el doble original. Si son iguales, devolvemos la cadena convertida cuya precisión es 15.
  3. De lo contrario, convierta el doble a cadena con una precisión de 17.

En este caso, double.ToString ("R") eligió incorrectamente el resultado con una precisión de 15 para que ocurra el error. Hay una solución oficial en el documento de MSDN:

En algunos casos, los valores dobles formateados con la cadena de formato numérico estándar "R" no se redondea con éxito si se compila utilizando / platform: x64 o / platform: anycpu cambia y se ejecuta en sistemas de 64 bits. Para evitar este problema, puede formatear valores dobles utilizando la cadena de formato numérico estándar "G17". El siguiente ejemplo usa la cadena de formato "R" con un valor Doble que no se redondea con éxito, y también usa la cadena de formato "G17" para redondea con éxito el valor original.

Por lo tanto, a menos que se resuelva este problema, debe usar double.ToString ("G17") para los viajes de ida y vuelta.

Actualización : ahora hay un problema específico para rastrear este error.

Jim Ma
fuente