Proposal: fix Enum Double instance

John Lato jwlato at gmail.com
Mon Jul 8 07:04:35 CEST 2013 

(apologies, meant this to go to libraries@ also)


On Mon, Jul 8, 2013 at 1:03 PM, John Lato <jwlato at gmail.com> wrote:

> Well, yes.  But at least we'd have a library function enumFromStepTo that
> works.  I suggest that anyone who cares about floating-point precision is
> going to be looking for a function with "Step" in the name for this anyway,
> so they're likely to find the most useful function first.
>
> Further I suggest that anyone who uses enumFromThenTo deserves what they
> get, particularly if they're using it on a type that doesn't admit a torsor
> (and consequently the user is relying upon some implicit compiler-defined
> torsor).
>
>
> On Mon, Jul 8, 2013 at 12:39 PM, Edward Kmett <ekmett at gmail.com> wrote:
>
>> Rather that would be the extension required if you wanted to implement
>> this as an overloadable method in the class.
>>
>> Implementing it separately and then wiring up enumFromThenTo in terms of
>> it still leaves you with the annoying initial massive cancellation of your
>> significand, so you wind up taking very precise steps of the wrong size. ;)
>>
>> -Edward
>>
>> On Sun, Jul 7, 2013 at 11:36 PM, Edward Kmett <ekmett at gmail.com> wrote:
>>
>>> Sadly this technically requires a language extension over Haskell
>>> 98/2010: ConstrainedClassMethods
>>>
>>>
>>> http://www.haskell.org/ghc/docs/6.12.2/html/users_guide/type-class-extensions.html#class-method-types
>>>
>>> GHC seems to be somewhat lax about requiring that extension on code in
>>> practice these days, but it does have implications when it comes to other
>>> approaches to implementing a compiler for Haskell.
>>>
>>> -Edward
>>>
>>>
>>> On Sun, Jul 7, 2013 at 11:18 PM, John Lato <jwlato at gmail.com> wrote:
>>>
>>>> My understanding was that the proposal required enumFromStepTo as a
>>>> separate function (i.e. not a class method) with a Num constraint on the
>>>> types.  So the function wouldn't be available for Ordering or derived Enum
>>>> instances unless they also had a Num instance.
>>>>
>>>> John
>>>>
>>>>
>>>> On Mon, Jul 8, 2013 at 10:46 AM, Edward Kmett <ekmett at gmail.com> wrote:
>>>>
>>>>> There is a pretty big potential problem with this.
>>>>>
>>>>> Not every type that is enumerable admits a torsor where you can talk
>>>>> about forward differences.
>>>>>
>>>>> What would the enumFromStepTo function mean for, say, Ordering or for
>>>>> most derived Enum instances?
>>>>>
>>>>> -Edward
>>>>>
>>>>> On Sun, Jul 7, 2013 at 12:17 PM, Edward A Kmett <ekmett at gmail.com>wrote:
>>>>>
>>>>>> +1 from me.
>>>>>>
>>>>>> The massive cancellation destroying significant figures of
>>>>>> significand in enumFromThenTo has always driven me nuts.
>>>>>>
>>>>>> Sent from my iPhone
>>>>>>
>>>>>> On Jul 5, 2013, at 8:50 AM, Twan van Laarhoven <twanvl at gmail.com>
>>>>>> wrote:
>>>>>>
>>>>>> > Hello All,
>>>>>> >
>>>>>> >
>>>>>> > I would like to make a counterproposal in the Enum Double debate:
>>>>>> Instead of deprecating or removing the instances, how about just fixing
>>>>>> them?
>>>>>> >
>>>>>> > A perfect instance for Enum Double is not possible, because
>>>>>> arithmetic is inexact. But you can actually get awfully close. I.e. instead
>>>>>> of allowing the final value to be at most step/2 past the end, we can allow
>>>>>> it to be at most about 2e-16*step past the end. In many practical
>>>>>> applications this is close enough to not be a problem.
>>>>>> >
>>>>>> >
>>>>>> >
>>>>>> > Now for some more detail on the design:
>>>>>> >
>>>>>> > * First of all, Haskell's enumFromThenTo is stupid, because
>>>>>> calculating step=then-from destroys numerical accuracy. So I will focus on
>>>>>> implementing a function enumFromStepTo. You can still implement
>>>>>> enumFromThenTo on top of it. But it might also make sense to expose
>>>>>> enumFromStepTo.
>>>>>> >
>>>>>> > * It is possible to write an exact instance for Rational. It is IMO
>>>>>> unacceptable that Haskell currently does not use this correct instance. I
>>>>>> will use this Rational instance as a baseline for comparison.
>>>>>> >
>>>>>> >    instance EnumFromStepTo Rational where
>>>>>> >      enumFromStepTo f s t
>>>>>> >        | s >= 0 && f <= t = f : enumFromStepTo (f + s) s t
>>>>>> >        | s <  0 && f >= t = f : enumFromStepTo (f + s) s t
>>>>>> >        | otherwise = []
>>>>>> >
>>>>>> > * Writing a loop naively, by recursing with from' = from+step will
>>>>>> result in accumulating the error of the addition. On the other hand,
>>>>>> Doubles can represent integer numbers exactly up to 2^53, so it is better
>>>>>> to use values from+i*step.
>>>>>> >
>>>>>> > * Assume for simplicity that step>0. The stopping condition will
>>>>>> then be of the form  from+i*step-to > 0. When this quantity is 0 then the
>>>>>> final value should be included, otherwise it shouldn't be.
>>>>>> >
>>>>>> > * Now for an error analysis. Assume that we start with
>>>>>> >   f,s,t :: Rational
>>>>>> >  and call
>>>>>> >   let f' = fromRational f
>>>>>> >   let s' = fromRational s
>>>>>> >   let t' = fromRational t
>>>>>> >   enumFromThenTo f' s' t'
>>>>>> >
>>>>>> > The fromRational function rounds a rational number to the nearest
>>>>>> representable Double. The relative error in f' is therefore bounded by abs
>>>>>> (f-f') ≤ ε * abs f, where ε is the half the maximum distance between two
>>>>>> adjacent doubles, which is about 1.11e-16. similarly for s' and t'. the
>>>>>> result of an addition or multiplication operation will also be rounded to
>>>>>> the nearest representable Double.
>>>>>> >
>>>>>> > So, the total error in our stopping condition is bounded by
>>>>>> > err(f+i*s-t)
>>>>>> >   ≤ ε * abs f             -- representing f
>>>>>> >   + ε * i * abs s         -- representing s, amplified by i
>>>>>> >   + ε * i * abs s         -- error in calculating (*), i * s
>>>>>> >   + ε * abs (f + i*s)     -- error in calculating (+), f + (i*s)
>>>>>> >   + ε * abs t             -- representing t
>>>>>> >   + ε * abs (f + i*s - t) -- error in calculating (-)
>>>>>> >
>>>>>> > Note that the last term of the bound is the thing we are bounding.
>>>>>> So we can handle that by picking a slightly larger epsilon, eps = ε/(1-ε).
>>>>>> >
>>>>>> > This leads to the following implementation of enumFromStepTo:
>>>>>> >
>>>>>> >    enumFromStepToEps eps !f !s !t = go 0
>>>>>> >      where
>>>>>> >      go i
>>>>>> >        | s >= 0 && x <= t = x : go (i+1)
>>>>>> >        | s <  0 && x >= t = x : go (i+1)
>>>>>> >        | abs (x - t) < eps * bound = [x]
>>>>>> >        | otherwise = []
>>>>>> >        where
>>>>>> >        x = f + i * s
>>>>>> >        bound = abs f + 2 * abs (i * s) + abs t + abs x
>>>>>> >
>>>>>> >  with eps = 1.12e-16 :: Double
>>>>>> >  or   eps = 5.97e-8 :: Float
>>>>>> >
>>>>>> > I have done extensive experiments with this implementation and many
>>>>>> variants, and I believe it will work in all cases.
>>>>>> >
>>>>>> > Now for enumFromTo, i.e. step=1 we could make a special case, since
>>>>>> we know that step is exactly equal to 1, so there is no representation
>>>>>> error. We could get rid of the multiplication, and replace i*s by 0 in the
>>>>>> error calculation.
>>>>>> >
>>>>>> > Another issue that remains is enumFromThenTo. If we take
>>>>>> >  step = then - from
>>>>>> > then the relative error in step is bounded by
>>>>>> >  |step' - step| ≤ ε * (|then| + |from| + |then - from|).
>>>>>> > This error gets multiplied by i, so it can unfortunately become
>>>>>> quite large.
>>>>>> >
>>>>>> > I think it would be best if we use the more general
>>>>>> >
>>>>>> >    enumFromStepToEps' !eps !f !s !sErr !t = go 0
>>>>>> >      where
>>>>>> >      go i
>>>>>> >        | s >= 0 && x <= t = x : go (i+1)
>>>>>> >        | s <  0 && x >= t = x : go (i+1)
>>>>>> >        | abs (x - t) < eps * bound = [x]
>>>>>> >        | otherwise = []
>>>>>> >        where
>>>>>> >        x = f + i * s
>>>>>> >        bound = abs f + abs (i * s) + abs (i * sErr) + abs t + abs x
>>>>>> >
>>>>>> >   enumFromStepTo f s t = enumFromStepToEps' eps f s s t
>>>>>> >   enumFromTo     f   t = enumFromStepToEps' eps f 1 0 t
>>>>>> >   enumFromThenTo f h t = enumFromStepToEps' eps f (h - f)
>>>>>> >                              (abs h + abs f + abs (h - f)) t
>>>>>> >
>>>>>> >
>>>>>> > I have also looked at what other language implementations do. So
>>>>>> far I have found that:
>>>>>> > * Octave uses a similar method, with the bound
>>>>>> >     3 * double_epsilon * max (abs x) (abs t),
>>>>>> >   which is less tight than my bound.
>>>>>> >   see
>>>>>> http://hg.savannah.gnu.org/hgweb/octave/file/787de2f144d9/liboctave/array/Range.cc#l525
>>>>>> >
>>>>>> > * numpy just uses an array with length ceil((to-from)/step)
>>>>>> >   see
>>>>>> https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/ctors.c#L2742
>>>>>> >   It therefore suffers from numerical inaccuracies:
>>>>>> >   $ python
>>>>>> >   >>> from numpy import arange
>>>>>> >   >>> notone = 9/7. - 1/7. - 1/7.
>>>>>> >   >>> notone
>>>>>> >   1.0000000000000002
>>>>>> >   >>> len(arange(1,1))
>>>>>> >   0
>>>>>> >   >>> len(arange(1,notone))
>>>>>> >   1
>>>>>> >
>>>>>> >
>>>>>> > In summary:
>>>>>> > * change Enum Rational to the always correct implementation
>>>>>> > * change Enum Double and Enum Float to the almost-always correct
>>>>>> implementation propose above.
>>>>>> > * (optional) expose the enumFromStepTo function.
>>>>>> >
>>>>>> >
>>>>>> >
>>>>>> > Twan
>>>>>> >
>>>>>> > _______________________________________________
>>>>>> > Libraries mailing list
>>>>>> > Libraries at haskell.org
>>>>>> > http://www.haskell.org/mailman/listinfo/libraries
>>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> Libraries mailing list
>>>>> Libraries at haskell.org
>>>>> http://www.haskell.org/mailman/listinfo/libraries
>>>>>
>>>>>
>>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.haskell.org/pipermail/libraries/attachments/20130708/8e937668/attachment-0001.htm>

More information about the Libraries mailing list