From 03509643cb1a9d207fb790e30006a1d8d03d7efb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9my=20Oudompheng?= Date: Fri, 18 Jan 2019 00:02:26 +0100 Subject: [PATCH] =?UTF-8?q?Extend=20Ry=C5=AB=20algorithm=20to=20be=20usabl?= =?UTF-8?q?e=20in=20ParseFloat.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The new function reuses the same elements to compute a floating point representation of an input decimal number. The allowed decimal mantissas are those which fit in 55 bits (16-17 decimal digits). Change-Id: I1c607f93069971b9efb7b1154d754aacd8a7b2eb --- src/strconv/atof.go | 29 +++++-- src/strconv/export_test.go | 47 +++++++++++ src/strconv/extfloat2.go | 146 ++++++++++++++++++++++++++++++++++ src/strconv/extfloat2_test.go | 127 ++++++++++++++++++++++++++++- src/strconv/ftoahard_test.go | 4 +- 5 files changed, 341 insertions(+), 12 deletions(-) diff --git a/src/strconv/atof.go b/src/strconv/atof.go index ada85e9fed6926..ba979e25d7dea6 100644 --- a/src/strconv/atof.go +++ b/src/strconv/atof.go @@ -490,14 +490,29 @@ func atof64(s string) (f float64, err error) { } } // Try another fast path. - ext := new(extFloat) - if ok := ext.AssignDecimal(mantissa, exp, neg, trunc, &float64info); ok { - b, ovf := ext.floatBits(&float64info) - f = math.Float64frombits(b) - if ovf { - err = rangeError(fnParseFloat, s) + if RyuEnabled && !trunc { + b, ovf, ok := RyuFromDecimal(mantissa, exp, &float64info) + if ok { + f = math.Float64frombits(b) + if neg { + f = -f + } + if ovf { + err = rangeError(fnParseFloat, s) + } + return f, err + } + } + if !RyuEnabled { + ext := new(extFloat) + if ok := ext.AssignDecimal(mantissa, exp, neg, trunc, &float64info); ok { + b, ovf := ext.floatBits(&float64info) + f = math.Float64frombits(b) + if ovf { + err = rangeError(fnParseFloat, s) + } + return f, err } - return f, err } } } diff --git a/src/strconv/export_test.go b/src/strconv/export_test.go index bc41225c051955..489f9422abf100 100644 --- a/src/strconv/export_test.go +++ b/src/strconv/export_test.go @@ -4,6 +4,10 @@ package strconv +import ( + "math" +) + var ( BitSizeError = bitSizeError BaseError = baseError @@ -52,3 +56,46 @@ func ShowDecimal(d *decimalSlice) string { exp := d.dp - 1 return string(d.d[0]) + "." + string(d.d[1:d.nd]) + "e" + Itoa(exp) } + +func OldAtof(mant uint64, exp int) float64 { + if f, ok := atof64exact(mant, exp, false); ok { + return f + } + + // Try another fast path. + ext := new(extFloat) + if ok := ext.AssignDecimal(mant, exp, false, false, &float64info); ok { + b, _ := ext.floatBits(&float64info) + return math.Float64frombits(b) + } + + var d decimal + d.Assign(mant) + d.dp += exp + b, _ := d.floatBits(&float64info) + return math.Float64frombits(b) +} + +// FastAtof optimistically performs Atof, and only tries the float64 and Grisu fast paths. +func FastAtof(mant uint64, exp int) (float64, bool) { + if f, ok := atof64exact(mant, exp, false); ok { + return f, true + } + + // Try another fast path. + ext := new(extFloat) + if ok := ext.AssignDecimal(mant, exp, false, false, &float64info); ok { + b, _ := ext.floatBits(&float64info) + return math.Float64frombits(b), true + } + + return 0, false +} + +func ReadFloat(s string) (mant uint64, exp int) { + mantissa, exp, _, trunc, ok := readFloat(s) + if !ok || trunc { + panic("readFloat failure") + } + return mantissa, exp +} diff --git a/src/strconv/extfloat2.go b/src/strconv/extfloat2.go index 962ddeb678a692..aa6b35994fb352 100644 --- a/src/strconv/extfloat2.go +++ b/src/strconv/extfloat2.go @@ -308,6 +308,128 @@ func RyuFixed(d *decimalSlice, mant uint64, exp int, prec int, flt *floatInfo) { return } +func RyuFromDecimal(mant uint64, exp int, flt *floatInfo) (fbits uint64, ovf, ok bool) { + // Conversion from decimal to binary floating-point + // can be achieved by reusing the same building blocks + // as the Ryū algorithm. + // + // Given a decimal mantissa, we can multiply by the requested + // power of ten using the same routines. The 64 bit result + // is guaranteed to be correctly truncated (floored), when + // the decimal mantissa fits in 55 bits. + // + // This covers 16-digit mantissas, and a few 17-digits values. + + bitlen := bits.Len64(mant) + if bitlen > 55 { + return 0, false, false // cannot handle values that large. + } + // Shift mantissa to be exactly 55 bits. + mant <<= uint(55 - bitlen) + e2 := bitlen - 55 + + // multiply by a power of 10. It is required to know + // whether the computation is exact. + var pow *extfloat128 // a representation of 10^q + switch { + case exp > 309: + return 0x7ff << 52, true, true + case exp < -342: + return 0, false, true + case exp > 0: + pow = &RyuPowersOfTen[exp] + case exp == 0: + // no multiply + case exp < 0: + pow = &RyuInvPowersOfTen[-exp] + } + // Is it an exact computation? + exact := false + switch { + case exp > 55: + // large positive powers of ten are not exact + case 54 >= exp && exp >= 0: + exact = true + case 0 > exp && exp >= -25: + // division by a power of ten might be exact + // if mantissas are multiples of 5 + if divisibleByPower5(mant, -exp) { + exact = true + } + default: + // division by 10^25 cannot be exact + // as 5^25 has 59 bits. + } + + // Compute Floor(x*10^q) + var di uint64 + var d0 bool + if exp == 0 { + di, d0 = mant, true + exact = true + } else { + di, d0 = ryuMultiply(mant, pow.Hi, pow.Lo) + e2 += pow.Exp + 64 + 55 + } + // If computation was an exact division, lower bits must be ignored. + if exp < 0 && exact { + d0 = true + } + // Is exponent too low? Shrink mantissa for denormals. + blen := bits.Len64(di) + e2 += blen - 1 + extra := uint(blen - 53) // number of lower bits to remove + if e2 < flt.bias+1 { + extra += uint(flt.bias + 1 - e2) + e2 = flt.bias + 1 + } + if extra > uint(blen) { + return 0.0, false, true + } + // Compute correct rounding. + extramask := uint64(1<>extra, di&extramask + roundUp := false + if exact { + // If we computed an exact product, d + 1/2 + // should round to d+1 if 'd' is odd. + roundUp = dfrac > 1<<(extra-1) || + (dfrac == 1<<(extra-1) && !d0) || + (dfrac == 1<<(extra-1) && d0 && di&1 == 1) + } else { + // otherwise, d+1/2 always rounds up because + // we truncated below. + roundUp = dfrac>>(extra-1) == 1 + } + if dfrac != 0 { + d0 = false + } + if roundUp { + di++ + } + + // Rounding might have added a bit; shift down. + if di == 2<>= 1 + e2++ + } + + // Infinities. + if e2-flt.bias >= 1<> 1 // clear sign bit + if bits>>52 == 2047 { + // only finite numbers + bits ^= (1 << 60) + } + x := math.Float64frombits(bits) + s := FormatFloat(x, 'g', -1, 64) + + // Compute decimal mantissa, exponent. + mant, exp := ReadFloat(s) + f1, ok1 := FastAtof(mant, exp) + b2, _, ok2 := RyuFromDecimal(mant, exp, &Float64info) + f2 := math.Float64frombits(b2) + if ok1 { + oldOk++ + } + if ok2 { + ryuOk++ + } + if ok1 && ok2 && f1 != f2 { + t.Fatalf("inconsistent results %s => %v %v", s, f1, f2) + } + } + t.Logf("%d successes with old fast paths", oldOk) + t.Logf("%d successes with Ryū", ryuOk) +} + func mantExp(x float64) (mant uint64, exp int) { const ( mantbits = 52 @@ -318,7 +386,7 @@ func TestRyuPowersOfTen(t *testing.T) { } } - for q := 0; q < 320; q++ { + for q := 0; q < 344; q++ { // negative exponents // Let's compute 2^128 * 8^q / 5^q, which is never an integer pow := big.NewInt(5) @@ -342,8 +410,8 @@ func TestRyuPowersOfTen(t *testing.T) { exp := sz - 256 - 4*q //t.Logf(`{Hi: 0x%016x, Lo: 0x%016x, Exp: %d},`, hi, lo, exp) expect := ExtFloat128{Hi: hi, Lo: lo, Exp: exp} - if RyuInvPowersOfTen[q] != expect { - t.Errorf("wrong entry") + if int(q) >= len(RyuInvPowersOfTen) || RyuInvPowersOfTen[q] != expect { + t.Errorf("wrong entry, wants %#v", expect) } } } @@ -544,3 +612,56 @@ func BenchmarkOldFixed(b *testing.B) { }) } } + +var ryuAtofBenches = []struct { + name string + mant uint64 + exp int +}{ + {"64Decimal", 33909, 0}, + {"64Float", 3397784, -4}, + {"64FloatExp", 509, 73}, + {"64Denormal", 6226662346353213, -324}, + // Almost halfway, with less than 1e-16 ulp difference + // with only 16 decimal digits. + {"64HalfwayHard1", 6808957268280643, 116}, // from ftoahard + {"64HalfwayHard2", 4334126125515466, -225}, // from ftoahard + // Only 3e-13*ulp larger than halfway between denormals, + {"64HalfwayDenormal", 168514038588815, -323}, + // Few digits, but 9.11691642378e-312 = 0x1ada385d67b.7fffffff5d9...p-1074 + // so naive, rounded 64-bit arithmetic is not enough to round it correctly. + {"64HalfwayDenormalShort", 911691642378, -323}, + // 1.62420278e-315 = 0x1398359e.7fffe022p-1074, + // should parsable using 64-bit arithmetic. + {"64HalfwayDenormalVeryShort", 162420278, -323}, + // https://www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/ + {"64Denormal", 22250738585072011, -324}, +} + +func BenchmarkRyuAtof(b *testing.B) { + for _, c := range ryuAtofBenches { + b.Run(fmt.Sprintf("%de%d", c.mant, c.exp), func(b *testing.B) { + var v uint64 + for i := 0; i < b.N; i++ { + f, ovf, ok := RyuFromDecimal(c.mant, c.exp, &Float64info) + if ovf || !ok { + b.Fatal("could not parse") + } + v = f + } + b.Logf("parsed to %v", math.Float64frombits(v)) + }) + } +} + +func BenchmarkOldAtof(b *testing.B) { + for _, c := range ryuAtofBenches { + b.Run(fmt.Sprintf("%de%d", c.mant, c.exp), func(b *testing.B) { + var f float64 + for i := 0; i < b.N; i++ { + f = OldAtof(c.mant, c.exp) + } + b.Logf("parsed to %v", f) + }) + } +} diff --git a/src/strconv/ftoahard_test.go b/src/strconv/ftoahard_test.go index 21dab35442d732..0f4b9a165dfdc3 100644 --- a/src/strconv/ftoahard_test.go +++ b/src/strconv/ftoahard_test.go @@ -62,7 +62,7 @@ func GenerateHardFloat64s() []float64 { if bits.Len64(y) == 54 { f := math.Ldexp(float64(y>>1), e) - fmt.Printf("f=ldexp(%d,%d)=%v, f+=(%d+%.3e)e%d\n", + _ = fmt.Sprintf("f=ldexp(%d,%d)=%v, f+=(%d+%.3e)e%d\n", y>>1, e, f, x, prec, -q) hards = append(hards, f) } @@ -77,7 +77,7 @@ func GenerateHardFloat64s() []float64 { x, y, prec = findFrac(a, b, bitlen) f := math.Ldexp(float64(y>>1), e) - fmt.Printf("f=ldexp(%d,%d)=%v, f+=(%d+%.3e)e%d\n", + _ = fmt.Sprintf("f=ldexp(%d,%d)=%v, f+=(%d+%.3e)e%d\n", y>>1, e, f, x, prec, -q) hards = append(hards, f) }