np.around uses an inexact algorithm, not mentioned in docstring #14391

Activity

1. ahaldane commented on Aug 29, 2019

   

   ahaldane

   on Aug 29, 2019

   MemberAuthor

   More actions

   Some brainstorming on if/how np.around might actually be fixed, rather than just changing the docstring. The goal is to get np.around to return the nearest floating-point value to the correctly rounded float.

   Instead of falling back to the (slow) string-formatting code like CPython, maybe there is a simpler computation to do: As noted in the link, the basic computation we need is done in decimal string-printing code (eg our dragon4 code) which uses BigInt manipulation to do decimal arithmetic of the binary float mantissa/exponent. That involves a large number of BigInt multiplies and divides.

   It feels like we can greatly simplify that if all we want is decimal rounding.

   Instead of so many mul/div in dragon4, it seems like we should be able to get away with one of each, as the basic computation we want to do is rint(x * (10**d))/10**d which has 1 mul and 1 div. Lets unpack a float64 number into integer exponent e and mantissa m, so x = m/2^(52-e) (lets ignore sign for now). Then x*10**d in the formula becomes (m*10**d)/2^(52-e). Then the rint involves flooring the numerator to the nearest multiple of 2^(52-e), which means just masking the lower bits of the numerator since the denom is a power of 2. Then we divide by 10**d, which I think we can do on the numerator only since any lower bits would get lost by renormalization anyway. So the computation is something like m' = maskbits(m * 10**d, 52-e)/10**d using integer arithmetic only, and we leave e alone.

   Here's a toy implementation using Python's builtin BigInts. It also has some extra steps to get rounding right.

   def round(x, d):
       # use Python's BigInt implementation to round mantissa
   
       s, m, e = unpack_sme(x)
   
       m *= 10**d
       m += 2**(52-e-1) # add 0.5 for correct rounding (can optimize: lower bits unchanged)
       m &= ~(2**(52-e)-1) # rint (round to integer, ie multiple of 2**(52-e))
       m += (10**d)//2  # add another 0.5 prec to get right // rounding (?)
       m //= 10**d
   
       # Uses: one MulPow10, one DivPow10, one shift/add op, one mask op, 
       # and one AddPow10/2.
   
       # may need some extra steps here to normalize float, especially for
       # negative d (this function designed for positive d so far)
   
       return pack_sme(s, m, e)
   
   def unpack_sme(x):
       "return sign, mantissa, exponent"
       x = np.float64(x).view(np.int64)
       m = (x&0xfffffffffffff)|(1<<52)
       e = ((x>>52)&0x7ff) - 1023
       s = x>>63
       return int(s), int(m), int(e)
   
   def pack_sme(s, m, e):
       x = np.int64(0)
       x |= s<<63
       x |= (e+1023)<<52
       x |= m&~(1<<52)
       return x.view(float64)
   
   round(12.34567, 2)  # 12.35
   round(56294995342131.5, 3) # 56294995342131.5

   I haven't worked through it rigorously to make sure it is mathematically correct, this is just to sketch the idea. We might be able to re-use the BigInt_MultiplyPow10 and the like from dragon4, or at least all the PowerOf10 tables of values.

   Or, if we limit the allowed values of decimals, we might get away with a very simple BigInt implementation.
2. mattip

   closed this as completedin #14392on Sep 1, 2019
3. maxwell-aladago

   added a commit that references this issue on Sep 6, 2019

   DOC: update np.around docstring with note about floating-point error

   ...

   This issue will close once commit a7a98fb is merged into the 'master' branch.

   a7a98fb