Mengapa konversi bolak-balik melalui string tidak aman untuk double?

185

Baru-baru ini saya harus membuat cerita bersambung ganda menjadi teks, dan kemudian mendapatkannya kembali. Nilai tersebut tampaknya tidak setara:

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

Tetapi menurut MSDN: Standard Numeric Format Strings , opsi "R" seharusnya menjamin keamanan bolak -balik.

Penentu format round-trip ("R") digunakan untuk memastikan bahwa nilai numerik yang dikonversi ke string akan diuraikan kembali ke dalam nilai numerik yang sama

Kenapa ini terjadi?

c# double tostring precision Philip Ding
sumber
6
Saya debugged di VS saya dan itu kembali benar di sini
Neel
19
Saya telah mereproduksi itu kembali palsu. Pertanyaan yang sangat menarik.
Jon Skeet
40
.net 4.0 x86 - true, .net 4.0 x64 - false
Ulugbek Umirov
25
Selamat telah menemukan bug yang mengesankan di .net.
Aron
14
@Casperah Round trip secara khusus dimaksudkan untuk menghindari ketidakkonsistenan floating point
Gusdor

Jawaban:

178

Saya menemukan bug.

.NET melakukan hal berikut dalam 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);

DoubleToNumbercukup sederhana - itu hanya panggilan _ecvt, yang ada di runtime 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;
    }
}

Ternyata _ecvtmengembalikan string 845512408225570.

Perhatikan nol yang tertinggal? Ternyata itu membuat semua perbedaan!
Ketika nol hadir, hasilnya benar-benar diuraikan kembali0.84551240822557006, yang merupakannomor asli Anda- jadi itu sama, dan karenanya hanya 15 digit yang dikembalikan.

Namun, jika saya memotong string pada nol ke 84551240822557, maka saya kembali 0.84551240822556994, yang bukan nomor asli Anda, dan karenanya akan mengembalikan 17 digit.

Bukti: jalankan kode 64-bit berikut (sebagian besar yang saya ekstrak dari Microsoft Shared Source CLI 2.0) di debugger Anda dan periksa vdi akhir 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;
}
pengguna541686
sumber
4
Penjelasan yang bagus +1. Kode ini dari shared-source-cli-2.0 kan? Ini adalah satu-satunya pemikiran yang saya temukan.
Soner Gönül
10
Saya harus mengatakan itu agak menyedihkan. String yang secara matematis sama (seperti yang dengan nol trailing, atau katakanlah 2.1e-1 vs 0.21) harus selalu memberikan hasil yang identik, dan string yang secara matematis dipesan harus memberikan hasil yang konsisten dengan pemesanan.
gnasher729
4
@ Mr.Lister: Kenapa "2.1E-1 tidak sama dengan 0.21 begitu saja"?
user541686
9
@ gnasher729: Saya agak setuju pada "2.1e-1" dan "0.21" ... tetapi string dengan nol tambahan tidak persis sama dengan yang tanpa - di yang sebelumnya, nol adalah angka yang signifikan dan menambahkan presisi.
cHao
4
@ cHao: Er ... itu menambah presisi, tetapi itu hanya memengaruhi bagaimana Anda memutuskan untuk menyelesaikan jawaban akhir jika sigfigs penting bagi Anda, bukan bagaimana komputer harus menghitung jawaban akhir di tempat pertama. Tugas komputer adalah menghitung segala sesuatu dengan presisi paling tinggi terlepas dari ukuran pengukuran aktual dari angka-angka tersebut; itu masalah programmer jika dia ingin membulatkan hasil akhir.
user541686
107

Sepertinya saya ini hanyalah bug. Harapan Anda sepenuhnya masuk akal. Saya telah mereproduksinya menggunakan .NET 4.5.1 (x64), menjalankan aplikasi konsol berikut yang menggunakan DoubleConverterkelas saya . DoubleConverter.ToExactStringmenunjukkan nilai persis yang diwakili oleh 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);
    }
}

Hasil dalam .NET:

0.84551240822557
0.845512408225570055719799711368978023529052734375
0.84551240822556994469749724885332398116588592529296875
False

Hasil dalam Mono 3.3.0:

0.84551240822557006
0.845512408225570055719799711368978023529052734375
0.845512408225570055719799711368978023529052734375
True

Jika Anda secara manual menentukan string dari Mono (yang berisi "006" di akhir), .NET akan mem-parsing itu kembali ke nilai aslinya. Kelihatannya masalahnya ada diToString("R") handling bukan pada parsing.

Seperti disebutkan dalam komentar lain, sepertinya ini khusus untuk berjalan di bawah CLR x64. Jika Anda mengkompilasi dan menjalankan penargetan kode di atas x86, tidak apa-apa:

csc /platform:x86 Test.cs DoubleConverter.cs

... Anda mendapatkan hasil yang sama dengan Mono. Akan menarik untuk mengetahui apakah bug itu muncul di bawah RyuJIT - Saya sendiri belum menginstalnya. Secara khusus, saya bisa membayangkan ini mungkin menjadi bug JIT, atau sangat mungkin bahwa ada implementasi yang berbeda dari internal double.ToStringberdasarkan arsitektur.

Saya sarankan Anda mengajukan bug di http://connect.microsoft.com

Jon Skeet
sumber
1
Jadi, Jon? Untuk mengkonfirmasi, apakah ini bug di JITer, sebaris ToString()? Ketika saya mencoba mengganti nilai kode keras dengan rand.NextDouble()dan tidak ada masalah.
Aron
1
Ya, itu pasti dalam ToString("R")konversi. Coba ToString("G32")dan perhatikan itu mencetak nilai yang benar.
user541686
1
@Ron: Saya tidak tahu apakah itu bug di JITter atau implementasi BCL spesifik x64. Saya sangat ragu bahwa itu sesederhana seperti inlining. Menguji dengan nilai acak tidak banyak membantu, IMO ... Saya tidak yakin apa yang Anda harapkan untuk diperagakan.
Jon Skeet
2
Apa yang terjadi menurut saya adalah format "round trip" menghasilkan nilai yang 0,498ulp lebih besar dari yang seharusnya, dan parsing logic terkadang keliru membulatkannya ke fraksi ulp kecil terakhir. Saya tidak yakin kode mana yang lebih saya salahkan, karena saya akan berpikir format "pulang pergi" akan menghasilkan nilai numerik yang dalam seperempat ULP benar secara numerik; parsing logic yang menghasilkan nilai dalam 0,75ulp dari apa yang ditentukan jauh lebih mudah daripada logika yang harus menghasilkan hasil dalam 0,502ulp dari apa yang ditentukan.
supercat
1
Situs web Jon Skeet sedang down? Saya menemukan bahwa sangat tidak mungkin saya ... kehilangan kepercayaan di sini.
Patrick M
2

Baru-baru ini, saya mencoba menyelesaikan masalah ini . Seperti yang ditunjukkan melalui kode , double.ToString ("R") memiliki logika berikut:

  1. Cobalah untuk mengubah ganda menjadi string dalam ketepatan 15.
  2. Ubah string kembali menjadi dua kali lipat dan bandingkan dengan yang asli ganda. Jika mereka sama, kami mengembalikan string yang dikonversi dengan ketepatan 15.
  3. Jika tidak, konversikan dobel ke string dengan presisi 17.

Dalam hal ini, double.ToString ("R") salah memilih hasil dalam presisi 15 sehingga bug terjadi. Ada solusi resmi di dokumen MSDN:

Dalam beberapa kasus, nilai ganda yang diformat dengan string format angka "R" standar tidak berhasil pulang-pergi jika dikompilasi menggunakan / platform: x64 atau / platform: anycpu beralih dan berjalan pada sistem 64-bit. Untuk mengatasi masalah ini, Anda dapat memformat nilai ganda dengan menggunakan string format numerik "G17" standar. Contoh berikut menggunakan string format "R" dengan nilai dobel yang tidak pulang pergi dengan sukses, dan juga menggunakan string format "G17" untuk berhasil memotong-pergi nilai awal.

Jadi, kecuali masalah ini diselesaikan, Anda harus menggunakan double.ToString ("G17") untuk bolak-balik.

Pembaruan : Sekarang ada masalah khusus untuk melacak bug ini.

Jim Ma
sumber