From mboxrd@z Thu Jan 1 00:00:00 1970 Received: (from majordomo@localhost) by pauillac.inria.fr (8.7.6/8.7.3) id VAA06889; Mon, 10 Nov 2003 21:36:35 +0100 (MET) X-Authentication-Warning: pauillac.inria.fr: majordomo set sender to owner-caml-list@pauillac.inria.fr using -f Received: from nez-perce.inria.fr (nez-perce.inria.fr [192.93.2.78]) by pauillac.inria.fr (8.7.6/8.7.3) with ESMTP id VAA06804 for ; Mon, 10 Nov 2003 21:36:34 +0100 (MET) Received: from moby.atcorp.com (moby.atcorp.com [204.72.172.2]) by nez-perce.inria.fr (8.11.1/8.11.1) with ESMTP id hAAKaX103239 for ; Mon, 10 Nov 2003 21:36:33 +0100 (MET) Received: from atcorp.com (seahorse.atcorp.com [204.72.172.13]) by moby.atcorp.com (8.11.6/8.11.2) with ESMTP id hAAKcvZ30254; Mon, 10 Nov 2003 14:38:57 -0600 Message-ID: <3FAFF6AD.4090009@atcorp.com> Date: Mon, 10 Nov 2003 14:35:57 -0600 From: Eric Dahlman User-Agent: Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US; rv:1.5) Gecko/20031013 Thunderbird/0.3 X-Accept-Language: en-us, en MIME-Version: 1.0 To: Brian Hurt CC: Ocaml Mailing List Subject: Re: [Caml-list] Rounding mode References: In-Reply-To: Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit X-Loop: caml-list@inria.fr X-Spam: no; 0.00; caml-list:01 rounding:01 raffalli:01 rounding:01 implemented:01 simplistic:01 multiplying:01 simplistic:01 solvable:01 analysts:99 figuring:01 differing:01 doubles:01 floats:01 tightly:01 Sender: owner-caml-list@pauillac.inria.fr Precedence: bulk Brian Hurt wrote: >On Mon, 10 Nov 2003, Eric Dahlman wrote: > > > >>Christophe Raffalli wrote: >> >> >> >>>I saw a previous discussion about rounding mode for OCaml (in 2000). >>> >>>As anyone implemented the functions to change the rounding mode for >>>floating point from OCaml (possibly as a patch to add primitives and >>>save a C function call) ? >>> >>>This is necessary to implement interval arithmetic ... >>> >>> >>Somewhat off topic but why is this necessary from a numerical math type >>of perspective. I am honestly curious as I don't see how this would >>interact with the calculation in a meaningful way. >> >> >> > >Simplistic discussion of why changing the rounding modes is usefull for >numeric programming. If you know numerical analysis, stop reading now >:-). > > I read it anyway ;-) [Snip a description of how floating point works.] >Consider, for example, multiplying two numbers- x and y. But you don't >have x and y exactly, you have (x + e) and (y + e) for some error amount >e. You multiply them, and you get (x + e)*(y + e) = x*y + x*e + y*e + >e*e. Now, e*e is vanishingly small, we can ignore it. And, assuming x >and y are both about 1.0, our answer is x*y + ~2*e- in other words, we've >about doubled the amount of error. Even starting at 1e-15, doubling it >every operation and you eat up all your precision surprisingly fast. > > >How can we tell how much error we have? Well, you can either do a whole >bunch of graduate-level mathematics, or you can run the same algorithm >with different rounding modes. If you always round towards positive >infinity, you'll always be adding an error term- i.e. it'll be (x + e) and >(y + e) with a result of x*y + 2e. If you always round to negative >infinity, you'll be subtracting the error term- (x - e) and (y - e) with a >result x*Y - 2e. Comparing the results gives you both a good idea of >where the answer really is, and how big your errors are. > > Now I don't know that it is this simple. I understand that by always rounding in one direction you are in effect calculating a bound on the error in that direction but that will only work in the most simplistic of cases. For instance if half of your values had their sign changed before being summed the sign of half the associated error terms would also change. The result would not give you any information about the effective lower/upper bound you thought you were calculating. It seems to me that this approach is only applicable to problems which would be trivially solvable with graduate level mathematics and of absolutely no value in any interesting cases. >Note, all is not necessarily lost. If instead we have (x + e) and (y - e) >a little bit of quick algebra shows that the result is actually *more* >accurate than either of the original inputs. The errors tend to cancel >instead of accumulate. Numerical analysts call algorithms whose errors >tend to cancel "stable", while algorithsm whose errors tend to accumulate >are "unstable". A good textbook on numerical methods will generally only >discuss stable algorithms- use 'em. > > This is incorrect because you are thinking about real numbers while working with floating point numbers, this is especially dangerous when you are operating near or below the epsilon for the particular floating point representation like we are now. If you were to look at the set of numbers which are precisely represented by a given floating point representation you would see that it is concentrated near 0 and that it tends to thin out as you move away from there. The distance between any two adjacent floating point numbers is highly variable and not uniform or "nice" by almost any measure. I really am having a hard time figuring out what you can say about the tiny e terms in this case as we move between regions with differing "floating point resolutions". Just to illustrate this the different "resolutions" the smallest positive double float which can be represented is 2.225 x 10^(-308) while the epsilon for doubles in general is 1.11 x 10^(-16) and can be thought of at the coarsest resolution. That is a big difference and is because floats are packed in tightly around 0. As for a numerically stable algorithm being able to produce a result which is more accurate than the original inputs that is total hogwash. You cannot ever have more accuracy than was present in your original input, ever. Nihil ex nihilo. >The effect can also be sort of acheived by using a random rounding method. >For each operation, the CPU picks one of it's rounding methods at random >to use (well, pseudo-random anyways). Now, while this will in general >increase the accuracy of floating point programs, it also means that every >time you run a numerical program, even over the same data, you will get a >slightly different answer. Newbies wondering why their program keeps >changing it's answer will probably rival newbies wondering why their stack >overflows pretty quickly. Which is why I'm hesitant to advocate doing >this to all programs. I still would like to be able to select the >rounding mode- either at compile time or at run time, which ever is >easier. Actually, run time is probably easier to implement and more >usefull (I can switch rounding modes for different parts of the program). > > As an oldby I would be disconcerted by this as well. All you are doing is trading one sort of numerical noise for another which you guess should be better behaved but in the absence of a neat proof to the contrary I don't believe is. There are just *SO* many sources of huger error terms in numerical calculations that the tiny rounding errors are really lost in the noise. I even tried to come up with a case where it would matter but that got so farcical that I gave up on it. What I was really wondering about when I first asked this question was if the other person had exactly that kind of example where they were using some clever trick to usefully misapply floating point calculations in their problem domain. -Eric ------------------- To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/ Beginner's list: http://groups.yahoo.com/group/ocaml_beginners