Generating equally-spaced floats with least rounding error

[Terry Reedy]

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."

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.

 

[Chris Angelico]

False

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

This may be because one oughtn't to use == with floats anyway.
Fraction(4728779608739021, 2251799813685248)

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]

     ((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!

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
.... return [((n-i)*start + i*stop)/n for i in range(n+1)]
........ 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]

 

[Steven D'Aprano]

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.

 

[Terry Reedy]

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, 2251799813685248). 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