Velocity Reviews > Generating equally-spaced floats with least rounding error

# Generating equally-spaced floats with least rounding error

Terry Reedy
Guest
Posts: n/a

 09-25-2011
On 9/25/2011 1:21 AM, Steven D'Aprano wrote:
> Terry Reedy wrote:
>
>> On 9/24/2011 9:53 AM, Steven D'Aprano wrote:

> [...]
>>>>>> [0.0 + i*width for i in range(]
>>> [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."

I do hope you did not stop with my lead-in sentence, and read to the
end, where I gave you most of the answer you were looking for, without
using the fractions module.

--
Terry Jan Reedy

Chris Angelico
Guest
Posts: n/a

 09-25-2011
On Sun, Sep 25, 2011 at 3:31 PM, Steven D'Aprano
<(E-Mail Removed)> wrote:
>>>> 2.1 == F(21, 10)

> False
>

Ahh, but that may just mean that Fraction doesn't implement an == that
compares with floats.

>>> 2.1==float(F(21,10))

True

This may be because one oughtn't to use == with floats anyway.

>>> F(2.1)==F(21,10)

False
>>> F(2.1)

Fraction(4728779608739021, 225179981368524

That denominator is a power of 2 (2**51 as it happens), which is
unsurprising for something created from an IEEE floating point value.
Rounding F(21,10) off to fit into float produces this same value.

ChrisA

Arnaud Delobelle
Guest
Posts: n/a

 09-25-2011
On 24 September 2011 23:10, Terry Reedy <(E-Mail Removed)> wrote:

> On 9/24/2011 10:18 AM, Arnaud Delobelle wrote:
>> Â* Â* Â*((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(]]

> ['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!

Of course *0.0 adds nothing in this case. I was suggesting a general
solution to the problem of partitioning the interval (a, b) into n
subintervals of equal length. As for 7*2.1/7 it is exactly the same
as 21/10 to 18 decimal places as the output of your own solution shows
below.

[...]
> 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']

Can you give an implementation for any interval (a, b) where a, b are
floats and any partition size n where n is a positive int?

My suggestion was a simple answer which tries to avoid the obvious
pitfalls that the naive approaches may fall into, as shown

>>> def bary(start, stop, n):

.... return [((n-i)*start + i*stop)/n for i in range(n+1)]
....
>>> def width(start, stop, n):

.... width = stop - start
.... return [start + i*width/n for i in range(n+1)]
....

Here is where the 'width' approach will fail:

>>> width(-1.0, 1e-20, 4)

[-1.0, -0.75, -0.5, -0.25, 0.0]
>>> bary(-1.0, 1e-20, 4)

[-1.0, -0.75, -0.5, -0.25, 1e-20]

--
Arnaud

Steven D'Aprano
Guest
Posts: n/a

 09-25-2011
Terry Reedy wrote:

> I do hope you did not stop with my lead-in sentence, and read to the
> end, where I gave you most of the answer you were looking for, without
> using the fractions module.

Yes, I read your entire post, thank you, and for my purposes I'm happy with
the fractions module.

--
Steven

Terry Reedy
Guest
Posts: n/a

 09-26-2011
On 9/25/2011 3:36 AM, Arnaud Delobelle wrote:
> On 24 September 2011 23:10, Terry Reedy<(E-Mail Removed)> wrote:

>> 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']

>
> Can you give an implementation for any interval (a, b) where a, b are
> floats and any partition size n where n is a positive int?

I was leaving that as an exercise for Stephen
but since he is satisfied with using Franctions ...
I accepted the challenge and (I believe) did it.

I first separated float-fractions conversions from generating equally
spaced fractions. While this does not make much if any difference if and
when one converts back to floats, testing with (21,10), for instance, as
input is much easier than with (2.1).as_integer_ratio() ==
(4728779608739021, 225179981368524. In some applications, small
integer ratios are wanted anyway instead of floats. The main
complication in the code is getting all outputs to have the smallest
possible common denominator across the entire series.

#! python3 -- frange.py
# Copyright 2011 Terry Jan Reedy: use with acknowledgement
'''Generate n+1 equally spaced fractions or float given two endpoints
and the number n of intervals'''

from fractions import gcd

def floatrange(a, b, n):
'''Yield floats a, b, and n-1 equally spaced floats in between.'''
for num,dem in fracrange(a.as_integer_ratio(),
b.as_integer_ratio(), n):
yield num/dem

def fracrange(frac1, frac2, n):
'''Yield fractions frac1, frac2 and n-1 equally spaced fractions in
between.
Fractions are represented as (numerator, denominator > 0) pairs.
For output, use the smallest common denominator of the inputs
that makes the numerator range an even multiple of n.
'''
n1, d1 = frac1
n2, d2 = frac2
dem = d1 * d2 // gcd(d1, d2)
start = n1 * (dem // d1)
stop = n2 * (dem // d2)
rang = stop - start
q, r = divmod(rang, n)
if r:
gcd_r_n = gcd(r, n)
m = n // gcd_r_n
dem *= m
start *= m
stop *= m
step = rang // gcd_r_n # rang * m // n
else:
step = q # if r==0: gcd(r,n)==n, m==1, rang//n == q
for num in range(start, stop+1, step):
yield num,dem

for (i,j), x in zip(fracrange((0,1), (21,10), 7), floatrange(0.0, 2.1, 7)):
print((i,j), '%20.18f' % (i/j ), '%20.18f' % x )
print()
for i,j in fracrange((1,10), (22,10), 7): print(i,j)
print()
for i,j in fracrange((1,5), (1,1), 6): print(i,j)

# output
(0, 10) 0.000000000000000000 0.000000000000000000
(3, 10) 0.299999999999999989 0.299999999999999989
(6, 10) 0.599999999999999978 0.599999999999999978
(9, 10) 0.900000000000000022 0.900000000000000022
(12, 10) 1.199999999999999956 1.199999999999999956
(15, 10) 1.500000000000000000 1.500000000000000000
(18, 10) 1.800000000000000044 1.800000000000000044
(21, 10) 2.100000000000000089 2.100000000000000089

1 10
4 10
7 10
10 10
13 10
16 10
19 10
22 10

3 15
5 15
7 15
9 15
11 15
13 15
15 15

--
Terry Jan Reedy