Generating equally-spaced floats with least rounding error

S

Steven D'Aprano

I'm trying to generate a sequence of equally-spaced numbers between a lower
and upper limit. Given arbitrary limits, what is the best way to generate a
list of equally spaced floats with the least rounding error for each point?

For example, suppose I want to divide the range 0 to 2.1 into 7 equal
intervals, then the end-points of each interval are:

(0.0)---(0.3)---(0.6)---(0.9)---(1.2)---(1.5)---(1.8)---(2.1)

and I'd like to return the values in the brackets. Using Decimal or Fraction
is not an option, I must use floats. If the exact value isn't representable
as a float, I'm okay with returning the nearest possible float.

The width of each interval is:

width = (2.1 - 0.0)/7

The relevant points can be calculated in either of two methods:

#1 add multiples of the width to the starting value, 0.

#2 subtract multiples of width from the ending value, 2.1.

(Repeated addition or subtraction should be avoided, as it increases the
amount of rounding error.)

Mathematically the two are equivalent, but due to rounding, they may not be.
Here's a run using Python 3.2:
[0.0 + i*width for i in range(8)]
[0.0, 0.3, 0.6, 0.8999999999999999, 1.2, 1.5, 1.7999999999999998, 2.1]

The 4th and 7th values have rounding errors, the rest are exact.

[2.1 - (7-i)*width for i in range(8)]
[0.0, 0.30000000000000027, 0.6000000000000001, 0.9000000000000001,
1.2000000000000002, 1.5, 1.8, 2.1]

The 2nd, 3rd, 4th and 5th values have rounding errors. Note that the 7th
value is exact here, but not above.

Is there a way to pick between methods #1 and #2 (or some combination of the
two) without human intervention so as to minimise the rounding error? Or is
there some other way to generate equally-spaced floats? My efforts at
googling have not been helpful.
 
C

Chris Angelico

#1 add multiples of the width to the starting value, 0.

#2 subtract multiples of width from the ending value, 2.1.

#3, go half and half - generate the lower half by adding to the lower
bound, and the upper half by subtracting from the upper bound. Not
sure if it'll help or not but it might be worth a shot.

ChrisA
 
V

Vlastimil Brom

2011/9/24 Chris Angelico said:
#3, go half and half - generate the lower half by adding to the lower
bound, and the upper half by subtracting from the upper bound. Not
sure if it'll help or not but it might be worth a shot.

ChrisA
--
Just a naive way:
#4 compute the values in both directions and use the average of the
results (of, course, given, there isn't a better way to distinguish
the values showing rounding errors automatically).

(I guess dealing manually with values obtained by .as_integer_ratio()
doesn't count as pure float operation...)

vbr
 
M

Mel

Steven said:
I'm trying to generate a sequence of equally-spaced numbers between a
lower and upper limit. Given arbitrary limits, what is the best way to
generate a list of equally spaced floats with the least rounding error for
each point?

For example, suppose I want to divide the range 0 to 2.1 into 7 equal
intervals, then the end-points of each interval are:

(0.0)---(0.3)---(0.6)---(0.9)---(1.2)---(1.5)---(1.8)---(2.1)

and I'd like to return the values in the brackets. Using Decimal or
Fraction is not an option, I must use floats. If the exact value isn't
representable as a float, I'm okay with returning the nearest possible
float.

The width of each interval is:

width = (2.1 - 0.0)/7

The relevant points can be calculated in either of two methods:

#1 add multiples of the width to the starting value, 0.

#2 subtract multiples of width from the ending value, 2.1.

(Repeated addition or subtraction should be avoided, as it increases the
amount of rounding error.)

Mathematically the two are equivalent, but due to rounding, they may not
be. Here's a run using Python 3.2:
[0.0 + i*width for i in range(8)]
[0.0, 0.3, 0.6, 0.8999999999999999, 1.2, 1.5, 1.7999999999999998, 2.1]

The 4th and 7th values have rounding errors, the rest are exact.

[2.1 - (7-i)*width for i in range(8)]
[0.0, 0.30000000000000027, 0.6000000000000001, 0.9000000000000001,
1.2000000000000002, 1.5, 1.8, 2.1]

The 2nd, 3rd, 4th and 5th values have rounding errors. Note that the 7th
value is exact here, but not above.

Is there a way to pick between methods #1 and #2 (or some combination of
the two) without human intervention so as to minimise the rounding error?
Or is there some other way to generate equally-spaced floats? My efforts
at googling have not been helpful.

When I've done this with ints (usually in small embedded systems) I've
always preferred

low_limit + (total_width * i) / intervals

since it does the rounding on the biggest numbers where proportional error
will be least, and it's guaranteed to hit the low_limit and high_limit
exactly (as exactly as they can be represented, anyway.)

Mel.
 
M

Mark Dickinson

I'm trying to generate a sequence of equally-spaced numbers between a lower
and upper limit. Given arbitrary limits, what is the best way to generatea
list of equally spaced floats with the least rounding error for each point?

For example, suppose I want to divide the range 0 to 2.1 into 7 equal
intervals, then the end-points of each interval are:

(0.0)---(0.3)---(0.6)---(0.9)---(1.2)---(1.5)---(1.8)---(2.1)

and I'd like to return the values in the brackets. Using Decimal or Fraction
is not an option, I must use floats. If the exact value isn't representable
as a float, I'm okay with returning the nearest possible float.

Can you explain why you're constrained not to use Fraction? Speed?

Using Fraction for intermediate calculations actually works perfectly
for this, since conversions from float to Fraction are exact, while
conversions from Fraction to float are correctly rounded. So if
you're using Python, you're not too bothered about efficiency, and you
want provably correctly-rounded results, why not use Fraction?
from fractions import Fraction
start, stop, n = 0.0, 2.1, 7
[float(Fraction(start) + i * (Fraction(stop) - Fraction(start)) / n) for i in range(n+1)]
[0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1]
 
S

Steven D'Aprano

Mark said:
Can you explain why you're constrained not to use Fraction? Speed?

Speed is important, but secondary to correctness. To be honest, I never even
thought of using fractions as an intermediate result -- I was thinking of
generating lists of Fractions. I think I can use this approach.

But out of curiosity, how would you do it using nothing but floats? Is there
a way?
 
M

Mark Dickinson

Speed is important, but secondary to correctness. To be honest, I never even
thought of using fractions as an intermediate result -- I was thinking of
generating lists of Fractions. I think I can use this approach.

Speed could be improved greatly by using integer arithmetic (with some
help from the float.as_integer_ratio method) instead of using Fraction
directly, though it would require writing quite a lot more code;
Fraction's many and slow gcd calls for normalization are a problem
here, and since all denominators will be powers of 2 they're largely
unnecessary.
But out of curiosity, how would you do it using nothing but floats? Is there
a way?

Hmm. Beyond writing my own multiple-precision integer library using
floats as the base type, I can't think of anything obvious. :)
 
S

Steven D'Aprano

Mark said:
Using Fraction for intermediate calculations actually works perfectly
for this, since conversions from float to Fraction are exact, while
conversions from Fraction to float are correctly rounded. So if
you're using Python, you're not too bothered about efficiency, and you
want provably correctly-rounded results, why not use Fraction?
from fractions import Fraction
start, stop, n = 0.0, 2.1, 7
[float(Fraction(start) + i * (Fraction(stop) - Fraction(start)) / n)
[for i in range(n+1)]
[0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1]


Ah, I knew it was too easy!
from fractions import Fraction as F
start, stop, n = 1, 3.1, 7
[float(F(start) + i*(F(stop)-F(start))/n) for i in range(n+1)] [1.0, 1.3, 1.6, 1.9000000000000001, 2.2, 2.5, 2.8000000000000003, 3.1]

start, stop, n = -1, 1.1, 7
[float(F(start) + i*(F(stop)-F(start))/n) for i in range(n+1)]
[-1.0, -0.7, -0.39999999999999997, -0.09999999999999996,
0.20000000000000004, 0.5000000000000001, 0.8, 1.1]
 
C

Chris Angelico

Ah, I knew it was too easy!

Try using Fraction for the start and stop too:
from fractions import Fraction as F
start,stop,n = F(0),F(21,10),7
[float(start+i*(stop-start)/n) for i in range(n+1)] [0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1]
[float(start+i*(stop-start)/n) for i in range(n+1)]
[-1.0, -0.7, -0.4, -0.1, 0.2, 0.5, 0.8, 1.1]

(Tested in Py 3.2 for Win)

ChrisA
 
M

Mark Dickinson

Ah, I knew it was too easy!
from fractions import Fraction as F
start, stop, n = 1, 3.1, 7
[float(F(start) + i*(F(stop)-F(start))/n) for i in range(n+1)]

[1.0, 1.3, 1.6, 1.9000000000000001, 2.2, 2.5, 2.8000000000000003, 3.1]

I believe that's still giving correctly-rounded results. Note that
the stop value of 3.1 isn't exactly 3.1 here: it's

3.100000000000000088817841970012523233890533447265625

So the 4th value above is the closest float to 4/7 * 1.0 + 3/7 *
3.100000000000000088817841970012523233890533447265625.
 
A

Arnaud Delobelle

Mark said:
Using Fraction for intermediate calculations actually works perfectly
for this, since conversions from float to Fraction are exact, while
conversions from Fraction to float are correctly rounded.  So if
you're using Python, you're not too bothered about efficiency, and you
want provably correctly-rounded results, why not use Fraction?
from fractions import Fraction
start, stop, n = 0.0, 2.1, 7
[float(Fraction(start) + i * (Fraction(stop) - Fraction(start)) / n)
[for i in range(n+1)]
[0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1]


Ah, I knew it was too easy!
from fractions import Fraction as F
start, stop, n = 1, 3.1, 7
[float(F(start) + i*(F(stop)-F(start))/n) for i in range(n+1)]
[1.0, 1.3, 1.6, 1.9000000000000001, 2.2, 2.5, 2.8000000000000003, 3.1]
start, stop, n = 1, 3.1, 7
[((n-i)*start + i*stop)/n for i in range(n+1)]
[1.0, 1.3, 1.5999999999999999, 1.9000000000000001, 2.2, 2.5,
2.8000000000000003, 3.1]
start, stop, n = -1, 1.1, 7
[float(F(start) + i*(F(stop)-F(start))/n) for i in range(n+1)]
[-1.0, -0.7, -0.39999999999999997, -0.09999999999999996,
0.20000000000000004, 0.5000000000000001, 0.8, 1.1]
start, stop, n = -1, 1.1, 7
[((n-i)*start + i*stop)/n for i in range(n+1)]
[-1.0, -0.7000000000000001, -0.39999999999999997,
-0.09999999999999996, 0.20000000000000004, 0.5, 0.8, 1.1]

On these examples, using fractions is no better than what I suggested
in my previous post.
 
S

Stefan Krah

Arnaud Delobelle said:
start, stop, n = -1, 1.1, 7
[float(F(start) + i*(F(stop)-F(start))/n) for i in range(n+1)]
[-1.0, -0.7, -0.39999999999999997, -0.09999999999999996,
0.20000000000000004, 0.5000000000000001, 0.8, 1.1]
start, stop, n = -1, 1.1, 7
[((n-i)*start + i*stop)/n for i in range(n+1)]
[-1.0, -0.7000000000000001, -0.39999999999999997,
-0.09999999999999996, 0.20000000000000004, 0.5, 0.8, 1.1]

On these examples, using fractions is no better than what I suggested
in my previous post.

Why not use Decimal if one needs exact endpoints? Something like:

def drange(start, stop, step=1):
if (step == 0):
raise ValueError("step must be != 0")
with localcontext() as ctx:
ctx.traps[Inexact] = True
x = start
cmp = -1 if step > 0 else 1
while x.compare(stop) == cmp:
yield float(x)
x += step

list(drange(Decimal(1), Decimal("3.1"), Decimal("0.3"))) [1.0, 1.3, 1.6, 1.9, 2.2, 2.5, 2.8]
list(drange(Decimal(-1), Decimal("1.1"), Decimal("0.3"))) [-1.0, -0.7, -0.4, -0.1, 0.2, 0.5, 0.8]
list(drange(Decimal(-1), Decimal("1.1"), Decimal("0.1823612873"))) [-1.0, -0.8176387127, -0.6352774254, -0.4529161381, -0.2705548508, -0.0881935635, 0.0941677238, 0.2765290111, 0.4588902984, 0.6412515857, 0.823612873, 1.0059741603]
list(drange(Decimal(-1), Decimal("1.1"), Decimal(1)/3))
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "<stdin>", line 10, in drange
File "/usr/lib/python3.2/decimal.py", line 1178, in __add__
ans = ans._fix(context)
File "/usr/lib/python3.2/decimal.py", line 1652, in _fix
context._raise_error(Inexact)
File "/usr/lib/python3.2/decimal.py", line 3836, in _raise_error
raise error(explanation)
decimal.Inexact: None



Stefan Krah
 
T

Terry Reedy

I'm trying to generate a sequence of equally-spaced numbers between a lower
and upper limit. Given arbitrary limits, what is the best way to generate a
list of equally spaced floats with the least rounding error for each point?

For example, suppose I want to divide the range 0 to 2.1 into 7 equal
intervals, then the end-points of each interval are:

(0.0)---(0.3)---(0.6)---(0.9)---(1.2)---(1.5)---(1.8)---(2.1)

and I'd like to return the values in the brackets. Using Decimal or Fraction
is not an option, I must use floats. If the exact value isn't representable
as a float, I'm okay with returning the nearest possible float.

The width of each interval is:

width = (2.1 - 0.0)/7

Calculating width is the fundamental problem. .3 cannot be exactly
represented in binary, and higher multiples with multiply the error.
The relevant points can be calculated in either of two methods:

#1 add multiples of the width to the starting value, 0.

#2 subtract multiples of width from the ending value, 2.1.

(Repeated addition or subtraction should be avoided, as it increases the
amount of rounding error.)

Mathematically the two are equivalent, but due to rounding, they may not be.
Here's a run using Python 3.2:
[0.0 + i*width for i in range(8)]
[0.0, 0.3, 0.6, 0.8999999999999999, 1.2, 1.5, 1.7999999999999998, 2.1]

The 4th and 7th values have rounding errors, the rest are exact

No they are not. Their errors are just smaller and not visible with 16
digits.
width = (2.1 - 0.0)/7
['%20.18f' % x for x in [0.0 + i*width for i in range(8)]]
['0.000000000000000000', '0.299999999999999989', '0.599999999999999978',
'0.899999999999999911', '1.199999999999999956', '1.500000000000000000',
'1.799999999999999822', '2.100000000000000089']

((n-i)*a + i*b)/n for i in range(n+1)
['%20.18f' % x for x in [((7-i)*0.0 + i*2.1)/7 for i in range(8)]]
['0.000000000000000000', '0.299999999999999989', '0.599999999999999978',
'0.900000000000000133', '1.199999999999999956', '1.500000000000000000',
'1.800000000000000266', '2.100000000000000089']

In the two places where this disagrees with the previous result, I
believe it is worse. The *0.0 adds nothing. You are still multiplying an
inexactly represented 2.1 by increasingly large numbers. Even 7*2.1/7 is
not exact!

My answer is the same as Mark's. Do exact calculation with fractions
(whether using Fraction or not) and convert to float. I believe the
least common multiple of a.as_integer_ratio[1], b.as_integer_ratio[1],
and n is the common denominator needed to convert the numerators to a
range() output. For the current problem, that is lcm(10,7) = 70.

The best you can do for this example is
['%20.18f' % (i/10 ) for i in range(0, 22, 3)]
['0.000000000000000000', '0.299999999999999989', '0.599999999999999978',
'0.900000000000000022', '1.199999999999999956', '1.500000000000000000',
'1.800000000000000044', '2.100000000000000089']

Notice that the last place errors are all less than 50, so printing to
16 places will make them appear 'exact'.

Without being fancy (to see than you can cancel 7), you get the same
output with
['%20.18f' % (i/70 ) for i in range(0, 148, 21)]
['0.000000000000000000', '0.299999999999999989', '0.599999999999999978',
'0.900000000000000022', '1.199999999999999956', '1.500000000000000000',
'1.800000000000000044', '2.100000000000000089']

In the two places where these disagree with the first (your original)
they are *better*, with absolute error in the last places of 22 versus
89 and 44 versus 178 for .9 and 1.8. The last place errors greater than
50 gave you the 'inexact' answers in your post, and indeed, they are not
the best possible.
 
S

Steven D'Aprano

Arnaud said:
On these examples, using fractions is no better than what I suggested
in my previous post.

I'm sorry, I can't see your previous post. Perhaps it isn't available on my
news provider?
 
S

Steven D'Aprano

Terry said:
On 9/24/2011 9:53 AM, Steven D'Aprano wrote: [...]
[0.0 + i*width for i in range(8)]
[0.0, 0.3, 0.6, 0.8999999999999999, 1.2, 1.5, 1.7999999999999998, 2.1]

The 4th and 7th values have rounding errors, the rest are exact

No they are not. Their errors are just smaller and not visible with 16
digits.

Pardon. I meant that they were as close to exact as is possible in binary
floats. With the exception of 0.8999... and 1.7999... I don't believe any
other float can be closer to the exact value.

I did mention that "If the exact value isn't representable as a float, I'm
okay with returning the nearest possible float." :)
 
S

Steven D'Aprano

Chris said:
Try using Fraction for the start and stop too:

If you look again at the code I used, I did. I just did it inside the list
comp.

from fractions import Fraction as F
start,stop,n = F(0),F(21,10),7
[float(start+i*(stop-start)/n) for i in range(n+1)] [0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1]
[float(start+i*(stop-start)/n) for i in range(n+1)]
[-1.0, -0.7, -0.4, -0.1, 0.2, 0.5, 0.8, 1.1]


Something looks a tad suspicious here... the two list comps are identical,
but give completely different results, even though you don't re-assign
start and stop between calls. You're not editing your results are you?
<wink>

But seriously, you were freaking me out there for a bit. I couldn't see why
pulling the conversion to fraction outside of the list comp was giving
different results. And then it hit me...
False

You're testing different numbers from me. Try again with F(2.1) as the stop
value.
 
J

John Ladasky

I'm trying to generate a sequence of equally-spaced numbers between a lower
and upper limit. Given arbitrary limits, what is the best way to generatea
list of equally spaced floats with the least rounding error for each point?

For example, suppose I want to divide the range 0 to 2.1 into 7 equal
intervals, then the end-points of each interval are:

(0.0)---(0.3)---(0.6)---(0.9)---(1.2)---(1.5)---(1.8)---(2.1)

and I'd like to return the values in the brackets. Using Decimal or Fraction
is not an option, I must use floats. If the exact value isn't representable
as a float, I'm okay with returning the nearest possible float.

Use numpy.
from numpy import mgrid
seq = mgrid[0.0:2.1:8j]
seq

array([ 0. , 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1])
print repr(x)

0.0
0.29999999999999999
0.59999999999999998
0.89999999999999991
1.2
1.5
1.7999999999999998
2.1000000000000001
 
C

Chris Angelico

On Sun, Sep 25, 2011 at 3:31 PM, Steven D'Aprano
the code I used, I did. I just did it inside the list> comp.
You did, but you created a Fraction from a float; my version used
Fraction(21,10) instead of (effectively) Fraction(2.1). My description
was a little sloppy, but what I meant was to use Fraction for the
actual arguments to the "function" that this is implementing.
Something looks a tad suspicious here... the two list comps are identical,
but give completely different results, even though you don't re-assign
start and stop between calls. You're not editing your results are you?
<wink>

Whooops! I had a whole lot of other commands in the scrollback
(displaying intermediate results and such), and elided the rather
crucial reassignment of parameters along with them!

Fortunately, thanks to my reiplophobia, I still have the original
session up in IDLE. This is even the most recent thing in it.
sys.version '3.2 (r32:88445, Feb 20 2011, 21:29:02) [MSC v.1500 32 bit (Intel)]'
from fractions import Fraction as F
start,stop,n = F(0),F(21,10),7
[float(start+i*(stop-start)/n) for i in range(n+1)] [0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1]
start, stop, n = F(-1), F(11,10), 7
[float(start+i*(stop-start)/n) for i in range(n+1)]
[-1.0, -0.7, -0.4, -0.1, 0.2, 0.5, 0.8, 1.1]

There, now it's honest. Sorry about that!!

ChrisA
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Members online

No members online now.

Forum statistics

Threads
473,744
Messages
2,569,484
Members
44,903
Latest member
orderPeak8CBDGummies

Latest Threads

Top