Dlaczego konwersja w obie strony za pomocą łańcucha nie jest bezpieczna dla podwójnego?

185

Ostatnio musiałem przekształcić tekst podwójny w szereg, a potem go odzyskać. Wartość wydaje się nie być równoważna:

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

Ale zgodnie z MSDN: Standardowe ciągi formatu liczbowego opcja „R” ma gwarantować bezpieczeństwo w obie strony.

Specyfikator formatu w obie strony („R”) służy do zapewnienia, że ​​wartość liczbowa przekonwertowana na ciąg zostanie przeanalizowana z powrotem do tej samej wartości liczbowej

Dlaczego się to stało?

Philip Ding
źródło
6
Debugowałem mój VS i jego powrót tutaj
Neel,
19
Powtórzyłem to, zwracając wartość false. Bardzo interesujące pytanie.
Jon Skeet,
40
.net 4.0 x86 - true, .net 4.0 x64 - false
Ulugbek Umirov
25
Gratulujemy znalezienia tak imponującego błędu w .net.
Aron,
14
@Casperah Round trip ma na celu uniknięcie niespójności zmiennoprzecinkowych
Gusdor

Odpowiedzi:

178

Znalazłem błąd.

.NET wykonuje następujące czynności 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);

DoubleToNumberjest dość proste - po prostu wywołuje _ecvt, co jest w środowisku uruchomieniowym 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;
    }
}

Okazuje się, że _ecvtzwraca ciąg 845512408225570.

Zwróć uwagę na końcowe zero? Okazuje się, że robi różnicę!
Kiedy zero jest obecne, wynik faktycznie analizuje z powrotem0.84551240822557006, który jest twoim pierwotnym numerem - więc porównuje równy, a zatem zwracanych jest tylko 15 cyfr.

Jeśli jednak obetnę ciąg o tym zeru do 84551240822557, to wrócę 0.84551240822556994, co nie jest twoim pierwotnym numerem, a zatem zwróci 17 cyfr.

Dowód: uruchom następujący kod 64-bitowy (większość z nich wyodrębniłem z interfejsu Microsoft Shared Source CLI 2.0) w swoim debuggerze i sprawdź vna końcu 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;
}
użytkownik541686
źródło
4
Dobre wyjaśnienie +1. Ten kod pochodzi z shared-source-cli-2.0, prawda? To jedyna myśl, jaką znalazłem.
Soner Gönül,
10
Muszę powiedzieć, że to raczej żałosne. Ciągi, które są matematycznie równe (jak jeden z końcowym zerem lub, powiedzmy, 2.1e-1 vs. 0.21), zawsze powinny dawać identyczne wyniki, a ciągi, które są matematycznie uporządkowane, powinny dawać wyniki zgodne z kolejnością.
gnasher729
4
@MrLister: Dlaczego „2.1E-1 nie powinien być taki sam jak 0.21 tak po prostu”?
user541686,
9
@ gnasher729: W pewnym sensie zgadzam się na „2.1e-1” i „0.21” ... ale ciąg z zerowym końcem nie jest dokładnie równy zeru bez - w tym pierwszym zero jest cyfrą znaczącą i dodaje precyzja.
cHao
4
@ cHao: Eee ... dodaje precyzji, ale wpływa to tylko na to, jak zdecydujesz się zaokrąglić ostateczną odpowiedź, jeśli sigfigs są dla ciebie ważne, a nie to, w jaki sposób komputer powinien obliczyć ostateczną odpowiedź. Zadaniem komputera jest obliczenie wszystkiego z najwyższą precyzją, niezależnie od faktycznych dokładności pomiaru liczb; to problem programisty, jeśli chce zaokrąglić wynik końcowy.
user541686,
107

Wydaje mi się, że to po prostu błąd. Twoje oczekiwania są całkowicie uzasadnione. Odtworzyłem go przy użyciu .NET 4.5.1 (x64), uruchamiając następującą aplikację konsolową, która korzysta z mojej DoubleConverterklasy. DoubleConverter.ToExactStringpokazuje dokładną wartość reprezentowaną przez 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);
    }
}

Wyniki w .NET:

0.84551240822557
0.845512408225570055719799711368978023529052734375
0.84551240822556994469749724885332398116588592529296875
False

Wyniki w Mono 3.3.0:

0.84551240822557006
0.845512408225570055719799711368978023529052734375
0.845512408225570055719799711368978023529052734375
True

Jeśli ręcznie określisz ciąg z Mono (który zawiera „006” na końcu), .NET przeanalizuje go z powrotem do oryginalnej wartości. Wygląda na to, że problemem jest ToString("R")obsługa, a nie parsowanie.

Jak zauważono w innych komentarzach, wygląda na to, że jest to specyficzne dla działania w CLR x64. Jeśli skompilujesz i uruchomisz powyższy kod kierujący na x86, wszystko będzie w porządku:

csc /platform:x86 Test.cs DoubleConverter.cs

... masz takie same wyniki jak w przypadku Mono. Byłoby ciekawie wiedzieć, czy błąd pojawia się w RyuJIT - sam tego nie mam. W szczególności mogę sobie wyobrazić, że jest to błąd JIT, lub całkiem możliwe, że istnieją zupełnie różne implementacje elementów wewnętrznych double.ToStringopartych na architekturze.

Sugeruję zgłoszenie błędu na stronie http://connect.microsoft.com

Jon Skeet
źródło
1
Więc Jon? Aby potwierdzić, czy jest to błąd w JITerze, który zawiera ToString()? Próbowałem zastąpić wartość zakodowaną na stałe rand.NextDouble()i nie było problemu.
Aron,
1
Tak, zdecydowanie jest w trakcie ToString("R")konwersji. Spróbuj ToString("G32")zauważyć, że drukuje prawidłową wartość.
user541686,
1
@Aron: Nie wiem, czy to błąd w JITter, czy w specyficznej dla x64 implementacji BCL. Wątpię jednak, czy jest to tak proste, jak wprowadzenie. Testowanie z losowymi wartościami tak naprawdę niewiele pomaga, IMO ... Nie jestem pewien, czego się spodziewasz.
Jon Skeet,
2
Myślę, że to, co się dzieje, polega na tym, że format „podróży w obie strony” generuje wartość, która jest o 0,489ulp większa niż powinna być, a logika analizowania czasami błędnie zaokrągla ją w górę do tej ostatniej niewielkiej części wrzodu. Nie jestem pewien, który kod winię bardziej, ponieważ uważam, że format „w obie strony” powinien generować wartość liczbową, która mieści się w zakresie ULP wynoszącym jedną czwartą ULP; parsowanie logiki, która daje wartość w zakresie 0,75ulp od podanego, jest znacznie łatwiejsze niż logika, która musi dawać wynik w granicach 0,502ulp od podanej.
supercat
1
Strona Jona Skeeta nie działa? Uważam, że tak mało prawdopodobne , że ... tracę tutaj całą wiarę.
Patrick M,
2

Ostatnio próbuję rozwiązać ten problem . Jak wskazano w kodzie , double.ToString („R”) ma następującą logikę:

  1. Spróbuj przekonwertować podwójny na ciąg z dokładnością do 15.
  2. Przekształć ciąg z powrotem na podwójny i porównaj z oryginalnym podwójnym. Jeśli są takie same, zwracamy przekonwertowany ciąg, którego dokładność wynosi 15.
  3. W przeciwnym razie przekształć podwójny w ciąg z dokładnością do 17.

W takim przypadku double.ToString („R”) nieprawidłowo wybrał wynik z dokładnością do 15, więc błąd się zdarza. Dokument MSDN zawiera oficjalne obejście:

W niektórych przypadkach podwójne wartości sformatowane za pomocą standardowego ciągu liczbowego formatu „R” nie działają poprawnie w obie strony, jeśli są kompilowane przy użyciu przełączników / platform: x64 lub / platform: anycpu i działają w systemach 64-bitowych. Aby obejść ten problem, możesz sformatować wartości podwójne przy użyciu standardowego ciągu formatu liczbowego „G17”. W poniższym przykładzie użyto ciągu formatu „R” z wartością Double, która nie powraca w obie strony, a także łańcucha „G17” w celu pomyślnego przejścia w obie strony oryginalnej wartości.

Tak więc, chyba że problem zostanie rozwiązany, należy użyć funkcji double.ToString („G17”) w celu przejścia w obie strony.

Aktualizacja : teraz istnieje konkretny problem do śledzenia tego błędu.

Jim Ma
źródło