hi, I think this is the well-known issue of using real numbers in decimal representation on a machine that thinks binary, but I don't know what to do with it, and some of you maybe do. I want to shift+stretch a list of doubles into a given interval. example: | x1 = [2, 3, 4, 5, 10] | y1 = normInterval x1 0 1 | => y1 = [0.0,0.125,0.25,0.375,1.0] The function that does this looks something like this: | normInterval :: [Double] -> Double -> Double -> [Double] | normInterval ps lower upper = map (\ x -> (x - oldLower) * stretch + lower) ps | where | oldLower = head ps | oldUpper = last ps | stretch = (upper - lower) / (oldUpper - oldLower) However, with bigger numbers I get rounding errors: | x2 = [0.0,1.9569471624266144e-3,5.870841487279843e-3,1.5655577299412915e-2,3.913894324853229e-2,9.393346379647749e-2,0.2191780821917808,0.5009784735812133,1.1272015655577299,2.504892367906066] | y2 = normInterval x2 0 1 | => y2 = [0.0,7.8125e-4,2.3437500000000003e-3,6.25e-3,1.5625000000000003e-2,3.7500000000000006e-2,8.750000000000001e-2,0.2,0.45000000000000007,0.9999999999999999] The solution that pops to my mind is very simple: | normInterval :: [Double] -> Double -> Double -> [Double] | normInterval ps lower upper = repair (map (\ x -> (x - oldLower) * stretch + lower) ps) | where | oldLower = head ps | oldUpper = last ps | stretch = (upper - lower) / (oldUpper - oldLower) | | -- fix rounding error: | repair [i] = [upper] | repair (h:t) = h : repair t It works, but it's ugly. Is there a better way to do this? Thanks, Matthias _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe signature.asc (196 bytes) Download Attachment |
Your solution works, but is slightly wasteful with (repair) traversing
the whole list again. Here is a slightly more efficient expression: -- Precondition: The first parameter (xs) is sorted (ascending) : -- assert (all (zipWith (<=) (xs, tail xs))) -- low' < high' -- low < high normInterval :: [Double] -> Double -> Double -> [Double] normInterval xs low high = let low' = head xs; high' = last xs; scale = (high-low)/(high'-low') middle = init (tail xs) in (low : [scale*(x-low')+low) | x <-middle])++[high] -- Chris _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe |
In reply to this post by Matthias Fischmann
Well, if you are relying on exact results from floating point
arithmetic you're in trouble no matter what you do. I would just ignore the slight error and when finally printing the results do some rounding. Trying to fudge things is just going to bite you somewhere else. (BTW, I much prefer the name "floating point" to "real". I think the latter should be reserved for when you actually have real numbers.) -- Lennart Matthias Fischmann wrote: > > hi, > > I think this is the well-known issue of using real numbers in decimal > representation on a machine that thinks binary, but I don't know what > to do with it, and some of you maybe do. > > I want to shift+stretch a list of doubles into a given interval. > example: > > | x1 = [2, 3, 4, 5, 10] > | y1 = normInterval x1 0 1 > | => y1 = [0.0,0.125,0.25,0.375,1.0] > > The function that does this looks something like this: > > | normInterval :: [Double] -> Double -> Double -> [Double] > | normInterval ps lower upper = map (\ x -> (x - oldLower) * stretch + lower) ps > | where > | oldLower = head ps > | oldUpper = last ps > | stretch = (upper - lower) / (oldUpper - oldLower) > > However, with bigger numbers I get rounding errors: > > | x2 = [0.0,1.9569471624266144e-3,5.870841487279843e-3,1.5655577299412915e-2,3.913894324853229e-2,9.393346379647749e-2,0.2191780821917808,0.5009784735812133,1.1272015655577299,2.504892367906066] > | y2 = normInterval x2 0 1 > | => y2 = [0.0,7.8125e-4,2.3437500000000003e-3,6.25e-3,1.5625000000000003e-2,3.7500000000000006e-2,8.750000000000001e-2,0.2,0.45000000000000007,0.9999999999999999] > > The solution that pops to my mind is very simple: > > | normInterval :: [Double] -> Double -> Double -> [Double] > | normInterval ps lower upper = repair (map (\ x -> (x - oldLower) * stretch + lower) ps) > | where > | oldLower = head ps > | oldUpper = last ps > | stretch = (upper - lower) / (oldUpper - oldLower) > | > | -- fix rounding error: > | repair [i] = [upper] > | repair (h:t) = h : repair t > > It works, but it's ugly. Is there a better way to do this? > > > > Thanks, > Matthias > > > ------------------------------------------------------------------------ > > _______________________________________________ > Haskell-Cafe mailing list > [hidden email] > http://www.haskell.org/mailman/listinfo/haskell-cafe _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe |
In reply to this post by haskell-2
On Sun, Feb 26, 2006 at 01:00:54PM +0000, Chris Kuklewicz wrote: > To: Matthias Fischmann <[hidden email]> > Cc: [hidden email] > From: Chris Kuklewicz <[hidden email]> > Date: Sun, 26 Feb 2006 13:00:54 +0000 > Subject: Re: [Haskell-cafe] rounding errors with real numbers. > > Your solution works, but is slightly wasteful with (repair) traversing > the whole list again. Here is a slightly more efficient expression: > > -- Precondition: The first parameter (xs) is sorted (ascending) : > -- assert (all (zipWith (<=) (xs, tail xs))) > -- low' < high' > -- low < high > normInterval :: [Double] -> Double -> Double -> [Double] > normInterval xs low high = let low' = head xs; high' = last xs; > scale = (high-low)/(high'-low') > middle = init (tail xs) > in (low : [scale*(x-low')+low) | x <-middle])++[high] hi chris, i like this code better, too. slightly better still, because of fewer typos: normInterval :: [Double] -> Double -> Double -> [Double] normInterval ps lower upper = assert check ([lower] ++ [ stretch * (p - oldLower) + lower | p <- middle ] ++ [upper]) where check = lower < upper && oldLower < oldUpper && (and (zipWith (<=) ps (tail ps))) oldLower = head ps oldUpper = last ps stretch = (upper - lower) / (oldUpper - oldLower) middle = init (tail ps) but then i'll just take this as the best solution. thanks, matthias _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe signature.asc (196 bytes) Download Attachment |
In reply to this post by Lennart Augustsson
> Well, if you are relying on exact results from floating point
> arithmetic you're in trouble no matter what you do. As long as you don't do anything irrational (exp, sin, sqrt, etc.), you should be able to get away with using Rational. Number constants with decimals are not automatically constructors for floating point numbers in Haskell; they are exact (fromRational) until you make them Doubles or some other floating point value. If you have to use Doubles for other reasons (performance/memory, interfacing with other code, etc.) this won't help you... Just my 2c. Jared. -- http://www.updike.org/~jared/ reverse ")-:" _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe |
In reply to this post by Matthias Fischmann
Matthias Fischmann wrote:
>| -- fix rounding error: >| repair [i] = [upper] >| repair (h:t) = h : repair t Just to point out that this only fixes the last element of the list, so inputs like [1,2,10.8,10.8] would not be handled properly if you require the same input values to map to the same output values (I assume such inputs don't arise in the context you're using but in a general context the above wouldn't be a solution). Another thing is that when using floating point numbers, is there really any difference between 1.0 and 0.9999999 anyway? It's usually not recommended to ever test floats for equality since, depending on the architecture, the "same" float can end up being represented differently depending on what optimizations are happening eg an implementation could conceivably be making use of two different fp units if values are passed between different concurrently executing threads in a multi-processor or distributed processing environment... Regards, Brian. _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe |
In reply to this post by Matthias Fischmann
On Sun, 26 Feb 2006, Matthias Fischmann wrote: > I think this is the well-known issue of using real numbers in decimal > representation on a machine that thinks binary, but I don't know what > to do with it, and some of you maybe do. > > I want to shift+stretch a list of doubles into a given interval. > example: > > | x1 = [2, 3, 4, 5, 10] > | y1 = normInterval x1 0 1 > | => y1 = [0.0,0.125,0.25,0.375,1.0] > > The function that does this looks something like this: > > | normInterval :: [Double] -> Double -> Double -> [Double] > | normInterval ps lower upper = map (\ x -> (x - oldLower) * stretch + lower) ps Is there --------------------------------------------------------------^ a cancellation problem? Maybe you should use a kind of convex combination, that is (x-oldLower)*a + (oldUpper-x)*b _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe |
On Mon, Feb 27, 2006 at 03:11:35PM +0100, Henning Thielemann wrote: > To: Matthias Fischmann <[hidden email]> > cc: [hidden email] > From: Henning Thielemann <[hidden email]> > Date: Mon, 27 Feb 2006 15:11:35 +0100 (MET) > Subject: Re: [Haskell-cafe] rounding errors with real numbers. > > > On Sun, 26 Feb 2006, Matthias Fischmann wrote: > > >I think this is the well-known issue of using real numbers in decimal > >representation on a machine that thinks binary, but I don't know what > >to do with it, and some of you maybe do. > > > >I want to shift+stretch a list of doubles into a given interval. > >example: > > > >| x1 = [2, 3, 4, 5, 10] > >| y1 = normInterval x1 0 1 > >| => y1 = [0.0,0.125,0.25,0.375,1.0] > > > >The function that does this looks something like this: > > > >| normInterval :: [Double] -> Double -> Double -> [Double] > >| normInterval ps lower upper = map (\ x -> (x - oldLower) * stretch + > >lower) ps > > Is there --------------------------------------------------------------^ > a cancellation problem? > Maybe you should use a kind of convex combination, that is > > (x-oldLower)*a + (oldUpper-x)*b i don't quite understand this either. is 'x' an old element in my input list and your expression is the corresponding new element? then how does the resulting curve relate to the original one? does this keep the ratios between distances between elements in the list intact? (this is the property that i am interested in.) but it sounds intriguing. perhaps i should play with this a little and find out myself. thanks, matthias _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe signature.asc (196 bytes) Download Attachment |
On Mon, 27 Feb 2006, Matthias Fischmann wrote: > On Mon, Feb 27, 2006 at 03:11:35PM +0100, Henning Thielemann wrote: >> >> Is there a cancellation problem? > > what's a cancellation problem? zu deutsch: Auslöschung cancellation happens for instance here: 1 + 1e-50 - 1 == 0 >> Maybe you should use a kind of convex combination, that is >> >> (x-oldLower)*a + (oldUpper-x)*b > > i don't quite understand this either. You have to adjust 'a' and 'b' in order to obtain 0 on x==oldLower and 1 on x==oldUpper. The expression still depends linearly on x (plus an absolute term). Maybe it reduces the cancellations if the lower and the upper bound differ much in magnitude. _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe |
In reply to this post by Henning Thielemann
Henning Thielemann wrote:
> Maybe you should use a kind of convex combination, that is > > (x-oldLower)*a + (oldUpper-x)*b Maybe lower*(1-z) + upper*z where z = (x-oldLower) / (oldUpper-oldLower). I think this will always map oldLower and oldUpper to lower and upper exactly, but I'm not sure it's monotonic. -- Ben _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe |
In reply to this post by Henning Thielemann
> cancellation happens for instance here: 1 + 1e-50 - 1 == 0 the function again (in the wasteful original form, for clarity). where do you think cancellation may take place? isn't what you call canellation a generic rounding error? ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ normInterval :: [Double] -> Double -> Double -> [Double] normInterval ps lower upper = repair (map (\ x -> (x - oldLower) * stretch + lower) ps) where oldLower = head ps oldUpper = last ps stretch = (upper - lower) / (oldUpper - oldLower) -- fix rounding error: repair [i] = [upper] repair (h:t) = h : repair t ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ > >>(x-oldLower)*a + (oldUpper-x)*b i got this into my head though. neat. thanks! i will rewrite the function right now. cheers, matthias _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe signature.asc (196 bytes) Download Attachment |
On Thu, 2 Mar 2006, Matthias Fischmann wrote: > > cancellation happens for instance here: 1 + 1e-50 - 1 == 0 > > the function again (in the wasteful original form, for clarity). > where do you think cancellation may take place? isn't what you call > canellation a generic rounding error? Cancellation is a special kind of rounding error. Rounding errors appear everywhere, in (*), sin, exp and so on, but cancellations only arise on differences. They are especially bad, because as the example above shows, even if all numbers are given in double precision in the computation a+b-a, no digit of the result is correct, that is 100% rounding error! The danger of cancellation is everywhere where you subtract numbers of similar magnitude. In your example: If lower=-100000, upper=1, x approximately oldUpper then you get the effect at the '+' in ((x - oldLower) * stretch + lower). _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe |
On Fri, Mar 03, 2006 at 01:29:06PM +0100, Henning Thielemann wrote: > To: Matthias Fischmann <[hidden email]> > cc: [hidden email] > From: Henning Thielemann <[hidden email]> > Date: Fri, 3 Mar 2006 13:29:06 +0100 (MET) > Subject: Re: [Haskell-cafe] rounding errors with real numbers. > > > On Thu, 2 Mar 2006, Matthias Fischmann wrote: > > > > cancellation happens for instance here: 1 + 1e-50 - 1 == 0 > > > > the function again (in the wasteful original form, for clarity). > > where do you think cancellation may take place? isn't what you call > > canellation a generic rounding error? > > Cancellation is a special kind of rounding error. Rounding errors appear > everywhere, in (*), sin, exp and so on, but cancellations only arise on > differences. They are especially bad, because as the example above shows, > even if all numbers are given in double precision in the computation > a+b-a, no digit of the result is correct, that is 100% rounding error! The > danger of cancellation is everywhere where you subtract numbers of similar > magnitude. rounding error somewhere deep in the e-minus-somethings. how is the error getting worse than that, for which numbers? m. _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe signature.asc (196 bytes) Download Attachment |
On Fri, 3 Mar 2006, Matthias Fischmann wrote: > 1 + epsilon - 1 == epsilon, which is zero except for a very small > rounding error somewhere deep in the e-minus-somethings. how is the > error getting worse than that, for which numbers? I meant the relative error. epsilon should be the result, but the computer says 0, so the absolute error is epsilon and the relative error is 1, that is the number of reliable digits in the mantissa of the result is 0. _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe |
ok, now i am intimidated enough to give up. (-:. thanks for trying, though. m. On Fri, Mar 03, 2006 at 03:19:44PM +0100, Henning Thielemann wrote: > To: Matthias Fischmann <[hidden email]> > cc: [hidden email] > From: Henning Thielemann <[hidden email]> > Date: Fri, 3 Mar 2006 15:19:44 +0100 (MET) > Subject: Re: [Haskell-cafe] rounding errors with real numbers. > > > On Fri, 3 Mar 2006, Matthias Fischmann wrote: > > > 1 + epsilon - 1 == epsilon, which is zero except for a very small > > rounding error somewhere deep in the e-minus-somethings. how is the > > error getting worse than that, for which numbers? > > I meant the relative error. epsilon should be the result, but the computer > says 0, so the absolute error is epsilon and the relative error is 1, that > is the number of reliable digits in the mantissa of the result is 0. _______________________________________________ Haskell-Cafe mailing list [hidden email] http://www.haskell.org/mailman/listinfo/haskell-cafe signature.asc (196 bytes) Download Attachment |
Free forum by Nabble | Edit this page |