atofeisel.go

     1  // Copyright 2020 The Go Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package strconv
     6  
     7  // This file implements the Eisel-Lemire ParseFloat algorithm, published in
     8  // 2020 and discussed extensively at
     9  // https://nigeltao.github.io/blog/2020/eisel-lemire.html
    10  //
    11  // The original C++ implementation is at
    12  // https://github.com/lemire/fast_double_parser/blob/644bef4306059d3be01a04e77d3cc84b379c596f/include/fast_double_parser.h#L840
    13  //
    14  // This Go re-implementation closely follows the C re-implementation at
    15  // https://github.com/google/wuffs/blob/ba3818cb6b473a2ed0b38ecfc07dbbd3a97e8ae7/internal/cgen/base/floatconv-submodule-code.c#L990
    16  //
    17  // Additional testing (on over several million test strings) is done by
    18  // https://github.com/nigeltao/parse-number-fxx-test-data/blob/5280dcfccf6d0b02a65ae282dad0b6d9de50e039/script/test-go-strconv.go
    19  
    20  import (
    21  	"math/bits"
    22  )
    23  
    24  func eiselLemire64(man uint64, exp10 int, neg bool) (f float64, ok bool) {
    25  	// The terse comments in this function body refer to sections of the
    26  	// https://nigeltao.github.io/blog/2020/eisel-lemire.html blog post.
    27  
    28  	// Exp10 Range.
    29  	if man == 0 {
    30  		if neg {
    31  			f = float64frombits(0x8000000000000000) // Negative zero.
    32  		}
    33  		return f, true
    34  	}
    35  	pow, exp2, ok := pow10(exp10)
    36  	if !ok {
    37  		return 0, false
    38  	}
    39  
    40  	// Normalization.
    41  	clz := bits.LeadingZeros64(man)
    42  	man <<= uint(clz)
    43  	retExp2 := uint64(exp2+64-float64Bias) - uint64(clz)
    44  
    45  	// Multiplication.
    46  	xHi, xLo := bits.Mul64(man, pow.Hi)
    47  
    48  	// Wider Approximation.
    49  	if xHi&0x1FF == 0x1FF && xLo+man < man {
    50  		yHi, yLo := bits.Mul64(man, pow.Lo)
    51  		mergedHi, mergedLo := xHi, xLo+yHi
    52  		if mergedLo < xLo {
    53  			mergedHi++
    54  		}
    55  		if mergedHi&0x1FF == 0x1FF && mergedLo+1 == 0 && yLo+man < man {
    56  			return 0, false
    57  		}
    58  		xHi, xLo = mergedHi, mergedLo
    59  	}
    60  
    61  	// Shifting to 54 Bits.
    62  	msb := xHi >> 63
    63  	retMantissa := xHi >> (msb + 9)
    64  	retExp2 -= 1 ^ msb
    65  
    66  	// Half-way Ambiguity.
    67  	if xLo == 0 && xHi&0x1FF == 0 && retMantissa&3 == 1 {
    68  		return 0, false
    69  	}
    70  
    71  	// From 54 to 53 Bits.
    72  	retMantissa += retMantissa & 1
    73  	retMantissa >>= 1
    74  	if retMantissa>>53 > 0 {
    75  		retMantissa >>= 1
    76  		retExp2 += 1
    77  	}
    78  	// retExp2 is a uint64. Zero or underflow means that we're in subnormal
    79  	// float64 space. 0x7FF or above means that we're in Inf/NaN float64 space.
    80  	//
    81  	// The if block is equivalent to (but has fewer branches than):
    82  	//   if retExp2 <= 0 || retExp2 >= 0x7FF { etc }
    83  	if retExp2-1 >= 0x7FF-1 {
    84  		return 0, false
    85  	}
    86  	retBits := retExp2<<float64MantBits | retMantissa&(1<<float64MantBits-1)
    87  	if neg {
    88  		retBits |= 0x8000000000000000
    89  	}
    90  	return float64frombits(retBits), true
    91  }
    92  
    93  func eiselLemire32(man uint64, exp10 int, neg bool) (f float32, ok bool) {
    94  	// The terse comments in this function body refer to sections of the
    95  	// https://nigeltao.github.io/blog/2020/eisel-lemire.html blog post.
    96  	//
    97  	// That blog post discusses the float64 flavor (11 exponent bits with a
    98  	// -1023 bias, 52 mantissa bits) of the algorithm, but the same approach
    99  	// applies to the float32 flavor (8 exponent bits with a -127 bias, 23
   100  	// mantissa bits). The computation here happens with 64-bit values (e.g.
   101  	// man, xHi, retMantissa) before finally converting to a 32-bit float.
   102  
   103  	// Exp10 Range.
   104  	if man == 0 {
   105  		if neg {
   106  			f = float32frombits(0x80000000) // Negative zero.
   107  		}
   108  		return f, true
   109  	}
   110  	pow, exp2, ok := pow10(exp10)
   111  	if !ok {
   112  		return 0, false
   113  	}
   114  
   115  	// Normalization.
   116  	clz := bits.LeadingZeros64(man)
   117  	man <<= uint(clz)
   118  	retExp2 := uint64(exp2+64-float32Bias) - uint64(clz)
   119  
   120  	// Multiplication.
   121  	xHi, xLo := bits.Mul64(man, pow.Hi)
   122  
   123  	// Wider Approximation.
   124  	if xHi&0x3FFFFFFFFF == 0x3FFFFFFFFF && xLo+man < man {
   125  		yHi, yLo := bits.Mul64(man, pow.Lo)
   126  		mergedHi, mergedLo := xHi, xLo+yHi
   127  		if mergedLo < xLo {
   128  			mergedHi++
   129  		}
   130  		if mergedHi&0x3FFFFFFFFF == 0x3FFFFFFFFF && mergedLo+1 == 0 && yLo+man < man {
   131  			return 0, false
   132  		}
   133  		xHi, xLo = mergedHi, mergedLo
   134  	}
   135  
   136  	// Shifting to 54 Bits (and for float32, it's shifting to 25 bits).
   137  	msb := xHi >> 63
   138  	retMantissa := xHi >> (msb + 38)
   139  	retExp2 -= 1 ^ msb
   140  
   141  	// Half-way Ambiguity.
   142  	if xLo == 0 && xHi&0x3FFFFFFFFF == 0 && retMantissa&3 == 1 {
   143  		return 0, false
   144  	}
   145  
   146  	// From 54 to 53 Bits (and for float32, it's from 25 to 24 bits).
   147  	retMantissa += retMantissa & 1
   148  	retMantissa >>= 1
   149  	if retMantissa>>24 > 0 {
   150  		retMantissa >>= 1
   151  		retExp2 += 1
   152  	}
   153  	// retExp2 is a uint64. Zero or underflow means that we're in subnormal
   154  	// float32 space. 0xFF or above means that we're in Inf/NaN float32 space.
   155  	//
   156  	// The if block is equivalent to (but has fewer branches than):
   157  	//   if retExp2 <= 0 || retExp2 >= 0xFF { etc }
   158  	if retExp2-1 >= 0xFF-1 {
   159  		return 0, false
   160  	}
   161  	retBits := retExp2<<float32MantBits | retMantissa&(1<<float32MantBits-1)
   162  	if neg {
   163  		retBits |= 0x80000000
   164  	}
   165  	return float32frombits(uint32(retBits)), true
   166  }
   167
View as plain text