[racket] Inconsistency of `in-range` and numerical issues

On 02/24/2015 01:11 PM, Konrad Hinsen wrote:
> On 24/02/2015 16:41, Laurent wrote:
>
>> I've discovered a rather troubling behaviour when using `in-range` with
>> floating point numbers, which I think is worth knowing in case you
>> hadn't consider the issue before:
>>
>> On my machine, I get the following:
>>
>> (length (for/list ([i (in-range .1 .7 .1)]) i)) ; 6
>> (length (for/list ([i (in-range .1 .8 .1)]) i)) ; 8 (!)
>>
>> But:
>> (length (for/list ([i (in-range 1/10 7/10 1/10)]) i)) ; 6
>> (length (for/list ([i (in-range 1/10 8/10 1/10)]) i)) ; 7
>>
>>
>> Would it be a good idea to safe-guard these kinds of cases directly in
>> `in-range`?
>
> The problem is an old one that already troubled programmers in the age
> of Fortran. I suspect there is no reasonable safe-guard, with
> "reasonable" meaning that it does what people expect in all situations.
>
> The only way to stay out of trouble is to avoid loops defined by an
> accumulating floating-point value. This means either don't use floats
> (write the loop over integers or rationals and convert to floats when
> using the loop index), or don't use accumulation (define your range by
> two points and the number of subdivisions, rather than the width of
> subintervals).

I should have chimed in to support this two days ago, but this is 
exactly the right answer. Here's what Konrad means by his first 
alternative (write the loop over integers or rationals):

   (length
    (for*/list ([i  (in-range 1 8 1)]
                [i  (in-value (* i 0.1))])
      i))

The second alternative is a little harder to get right because of 
fencepost errors [1]. Fortunately, Racket has a library function for it. 
Unfortunately, it's buried in `plot/utils`. Here it is in action:

 > (require (only-in plot/utils linear-seq))

 > (linear-seq 0.0 1.0 4)
'(0.0 0.3333333333333333 0.6666666666666666 1.0)

 > (linear-seq 0.0 1.0 4 #:start? #f)
'(0.14285714285714285 0.42857142857142855 0.7142857142857142 1.0)

 > (linear-seq 0.0 1.0 4 #:start? #f #:end? #f)
'(0.125 0.375 0.625 0.875)

I should really move this function into `math/base`.

If you must use a flonum step size, do something like this:

   (define (flonum-range start end step)
     (define n (exact-ceiling (/ (- end start) step)))
     (for/list ([i  (in-range 0 n)])
       (+ start (* i step))))

To get points with 0.5 ulp error each, which is the best you can do, add 
(require math/flonum) to your program and change the loop body to (flfma 
step (fl i) start).

Arguably, `in-range` should do something like the above when given 
floating-point step lengths. I don't know how feasible that is, though.

Neil ⊥

[1] http://en.wikipedia.org/wiki/Off-by-one_error#Fencepost_error