Proposal: fix Enum Double instance
Edward Kmett
ekmett at gmail.com
Mon Jul 8 06:39:11 CEST 2013
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/20130707/65c2d612/attachment.htm>
More information about the Libraries
mailing list