python - Rounding floats so that they sum to precisely 1

Question

I have a rather gnarly bit of code that must more-or-less randomly generate a bunch of percentages, stored as decimal floats. That is, it decides that material one makes up 13.307 percent of the total, then stores that in a dict as 0.13307.

The trouble is, I can never get those numbers to add up to exactly one. I'm not entirely certain what the problem is, honestly. It might be something to do with the nature of floats.

Here's the offending code, in all its overcomplicated glory:

while not sum(atmosphere.values())>=1:
    #Choose a material randomly
    themat=random.choice(list(materials.values()))

    #If the randomly chosen material is gaseous at our predicted temperature...
    if themat.vapor < temp:
        #Choose a random percentage that it will make up of our planet's atmosphere, then put it in the atmos dict.
        atmosphere[themat]=round(random.uniform(0.001,0.5),5)

#Find out if the fractions add up to more than 1
difference=(sum(atmosphere.values())-1)
#If one does...
while difference > 0:
    #Choose a random constituent
    themat=random.choice(list(atmosphere.keys()))
    #If that constituent has a higher fraction value than the amount we'd need to reduce the total to 1...
    if atmosphere[themat]>(sum(atmosphere.values())-1):
        #Subtract that much from it.
        atmosphere[themat]-=difference
        #Then break the loop, since we're done and otherwise we'd go on removing bits of the atmosphere forever.
        break
    else:
        #Otherwise, halve its percentage and reduce difference by the amount we reduced the atmosphere 
        oldperc=atmosphere[themat]
        atmosphere[themat]=oldperc/2
        difference-=oldperc/2

#Then, finally, we correct any overcorrections the previous block made.
difference=(sum(atmosphere.values())-1)
if difference < 0:
    #Choose a random mat
    themat=random.choice(list(atmosphere.keys()))
    #Then add enough to it that the total is 1.
    atmosphere[themat]+=difference

Sorry if I've missed something obvious, or am not providing an important bit of information, but I'm tired at the moment, and I've been trying to figure this out for days now.

Answer 1

From your code it looks like you're randomly generating planet atmospheres, presumably for some kind of game or something. At any rate, the randomness of it is convincing me it doesn't need to be too accurate.

So i'd suggest you don't use floats, just use ints and go up to 100. Then you'll get your exact summing. For any maths you want to use them in just cast.

Is this not an option?

If you insist on using floats, then read on...

The problem you have using floats is as follows:

A floating point (in this case a double) is represented like this:

which corresponds to a double of value:

So,

your number is (1+M) * 2**(E) (where E = e-offset)

1+M is always in the range 1-2.

So, we have equally spaced numbers inbetween each pair of power of two (positive and negative), and the spacing between the numbers doubles with each increase in the exponent, E.

Think about this, it means that there is a constant spacing of representable numbers between each of these numbers [(1,2),(2,4),(4,8), etc]. This also applies to the negative powers of two, so:

0.5 - 1
0.25 - 0.5
0.125 - 0.25
0.0625 - 0.125
etc.

And in each range, there are the same quantity of numbers. This means that if you take a number in the range (0.25,0.5) and add it to a number in the range (0.5,1), then you have a 50% chance that the number cannot be exactly represented.

If you sum two floating point numbers for which the exponents differ by D, then the chances of the sum being exactly representable are 2^-D.

If you then want to represent the range 0-1, then you'll have to be very careful about which floats you use (i.e. force the last N bits of the fraction to be zero, where N is a function of E).

If you go down this route, then you'll end up with far more floats at the top of the range than the bottom.

The alternative is to decide how close to zero you want to be able to get. Lets say you want to get down to 0.0001.

0.0001 = (1+M) * 2^E

log₂(0.0001) = -13.28771...

So we'll use -14 as our minimum exponent.

And then to get up to 1, we just leave the exponent as -1.

So now we have 13 ranges, each with twice as many values as the lower one which we can sum without having to worry about precision.

This also means though, that the top range has 213 more values we can use. This obviously isn't okay.

So, after picking a float, then round it to the nearest allowable value - in this case, by round I just mean set the last 13 bits to zero, and just put it all in a function, and apply it to your numbers immediately after you get them out of rand.

Something like this:

from ctypes import *

def roundf(x,bitsToRound):

    i = cast(pointer(c_float(x)), POINTER(c_int32)).contents.value

    bits = bin(i)

    bits = bits[:-bitsToRound] + "0"*bitsToRound

    i = int(bits,2)

    y = cast(pointer(c_int32(i)), POINTER(c_float)).contents.value

    return y

(images from wikipedia)