Path: utzoo!utgpu!jarvis.csri.toronto.edu!mailrus!uwm.edu!cs.utexas.edu!uunet!munnari.oz.au!cs.mu.oz.au!ok From: ok@cs.mu.oz.au (Richard O'Keefe) Newsgroups: comp.lang.prolog Subject: Re: Floating Point Message-ID: <2550@munnari.oz.au> Date: 27 Oct 89 09:09:44 GMT References: <2491@munnari.oz.au> <36169@srcsip.UUCP> <2524@munnari.oz.au> <6688@latcs1.oz> Sender: news@cs.mu.oz.au Lines: 248 In article <6688@latcs1.oz>, mahler@latcs1.oz (Daniel Mahler) writes: > Having both kinds of real unification is naturally preferable if > you have the computational resources (remember the discussion on swapping > and threaded code last year). Computational resources? Oh big deal. Look, it took me all of half an hour to write the code which follows at the end of this message. (And half of that was some testing.) The code runs in Quintus Prolog, and with a little bit of editing should work in SICStus Prolog and NU Prolog as well. > But if I had to choose I'd go fuzzy. > The reason is that the least significant digit(s) of reals are wildly > inaccurate. Speak for yourself. There are two quite different questions: (1) how accurate are your data (2) how accurate is the machine's arithmetic. As for (1), only you know. As for (2), on machines that conform to the IEEE standard, the answer is "accurate to the last bit". In fact, Peter Tang has recently published in ACM TOMS a method of calculating exp(x) whose error is around 0.47 ULP, mind-bogglingly good. With IEEE double precision arithmetic, I can add, subtract, multiply, divide, and get remainders of integers up to 4 x 10**15; I mean by this that the answers I get are *exact* for integers in that range. > If you have done high school physics you know about the evils > of over representing accuracy by not stating error and including > too many significant digits. This has nothing to do with floating point arithmetic. When you are publishing your results, you ought not to claim (even implicitly) more accuracy than your error bounds warrant. But the technical point here is a moral one: "Thou shalt not tell lies to thy scientific peers". No *numeric* evils follow from giving the computer all the digits you have. > The precision is not lost as the heavily numerical functions will > presumably be either system provided or interfaced in from another > language That will perhaps serve as an excuse for making Prolog floating-point arithmetic SLOW, it will *not* serve as an excuse for making it BAD. Why do you want to stop people developing floating-point code in Prolog? I assure you that there are people who would like to calculate coefficients for Runga-Kutta schemes, or would like to develop new coding methods to calculate special functions more precisely, ... who would be quite happy to use Prolog if only it had floating-point arithmetic they could TRUST. (The results of their research would be applied in other languages.) > The problem is that floating point arithmetic is conceptually > incompatible with unification, only identical objects should unify > (be made one) and 3.141593 is not the same object as pi. But this is my point! The number 3.141593 is equal to itself and so should unify with itself, but it should not unify with any other number. > The only interpretation of floating unification I can think of > is equal to the known accuracy. The "known" accuracy is "full machine precision". We are not just talking about numbers that you typed in, remember. We are talking about numbers which have been derived from long computations as well as numbers that you typed in. And indeed, there are functions f(-) such that the relative error in f(x) is considerably less than the relative error in x; you may supply an argument which is accurate to only 4 decimal places but get out a result which is accurate to 16 places. > Most of the time you should > be doing | X - Y | < tol rather than unifying, Just so. But since that is so, you don't *NEED* fuzz in unification. An extremely important point about this sort of test is that the appropriate tolerance depends on the PROBLEM as well as the machine arithmetic; there is no one number (quite contrary to what SB Prolog assumes) which can be used for all problems. I have exchanged E-mail with Cody (chairman of the IEEE 854 committee). We didn't actually discuss fuzzy unification, but something similar does happen in existing languages such as Ada. He said very strongly indeed that "sharp comparisons" are important (that is, non-fuzzy comparison), and suggested that a Prolog standard should "demand that floating-point comparison be exact". He has been thinking about such matters for years, and it would be foolish to disregard his advice. As I have said that fuzzy comparisons can be useful but do not belong in unification, I ought to show how in a language like Quintus Prolog or SICStus Prolog or NU Prolog you can get the effect of fuzzy comparison when you need it. It is really very simple, if you have either - the equivalent of ldexp() and frexp() in your Prolog, or - the IEEE functions scalb() and logb() in your Prolog, or - a C interface. It would be quite easy to build these operations into a Prolog system. (I wouldn't actually use ldexp() and frexp(), but would fiddle the exponent fields directly.) The bulk of the cost of an implementation would be adding the new built-ins to the compiler's tables. Note well: I am offering here *more* than you get from fuzzy unification, I am offering fuzzy less than and fuzzy greater than as well. The code has a Quintus copyright on it, but it was thrown together in half an hour and has received no additional testing. It is not warranted as fit for any purpose other than discussion. It may be used freely provided the copyright notices are preserved and if you use it in a product you must acknowledge its Quintus origin. (Help stamp out GNU restrictions.) ------------------------- cut here for fuzzy.pl ------------------------- % Package: fuzzy % Author : Richard A. O'Keefe % Updated: Friday October 27th, 1989. % Defines: Knuth-style fuzzy comparison predicates in Quintus Prolog. % This code is Copyright (C) 1989, Quintus Computer Systems, Inc., % but it may be used freely provided the copyight notice is preserved and % Quintus is acknowledged. :- module(fuzzy, [ 'fuzzy <'/2, % Def (21) p 218 sec 4.2.2 'fuzzy ~'/2, % Def (22) p 218 sec 4.2.2 'fuzzy >'/2, % Def (23) p 218 sec 4.2.2 'fuzzy ='/2, % Def (24) p 218 sec 4.2.2 set_fuzz/1 ]). sccs_id('"%Z%%E% %M% %I%"'). /* For a fixed epsilon, the fuzzy comparison predicates are defined in section 4.2.2 (Accuracy of floating point arithmetic) of Knuth, The Art of Computer Programming, vol 2, Seminumerical Algorithms (2nd edition) on page 218. (21) 'fuzzy <'(U, V) <-> V-U > epsilon.b^max(e_u, e_v) (22) 'fuzzy ~'(U, V) <-> |V-U| =< epsilon.b^max(e_u, e_v) (23) 'fuzzy >'(U, V) <-> U-V > epsilon.b^max(e_u, e_v) (24) 'fuzzy ='(U, V) <-> |V-U| =< epsilon.b^min(e_u, e_v) where U = F.b^e_u and |F| < 1 and V = G.b^e_v and |G| < 1 This is a very crude first implementation. */ 'fuzzy <'(U, V) :- 'QPFUZZCMP'(U, V, 2). 'fuzzy ~'(U, V) :- 'QPFUZZCMP'(U, V, C), C < 2. 'fuzzy >'(U, V) :- 'QPFUZZCMP'(U, V, 3). 'fuzzy ='(U, V) :- 'QPFUZZCMP'(U, V, 1). set_fuzz(Epsilon) :- number(Epsilon), Epsilon > 0, 'QPSETFUZZ'(Epsilon). foreign_file('fuzzy.o', [ 'QPFUZZCMP', 'QPSETFUZZ' ]). foreign('QPFUZZCMP', 'QPFUZZCMP'(+float,+float,[-integer])). foreign('QPSETFUZZ', 'QPSETFUZZ'(+float)). :- load_foreign_files(['fuzzy.o'], ['-lm']), abolish(foreign_file, 2), abolish(foreign, 2). ------------------------- cut here for fuzzy.c ------------------------- /* File : fuzzy.c Author : Richard A. O'Keefe Updated: Friday October 27th, 1989. Support: for library(fuzzy) -- Knuth-style fuzzy comparison This code is Copyright (C) 1989, Quintus Computer Systems, Inc., but it may be used freely provided the copyight notice is preserved and Quintus is acknowledged. */ #ifndef lint static char SCCSid[] = "%Z%%E% %M% %I%"; #endif/*lint*/ #define PRETTY_CLOSE 0 #define VERY_CLOSE 1 #define LESS_THAN 2 #define GREATER_THAN 3 #ifdef __STDC__ extern double frexp(double, int*); extern double ldexp(double, int); #else extern double frexp(); extern double ldexp(); #endif static double epsilon = 1.0e-5; void QPSETFUZZ(x) double x; { epsilon = x; } int QPFUZZCMP(U, V) double U, V; { double t; int e_u, e_v; if (U < 0.0 && V < 0.0) { /* U and V are both negative. We want them positive */ /* The order of U and V is the same as the order of -V and -U */ /* so we have to exchange them as well as change sign *. t = -U, U = -V, V = t; } else if (U <= 0.0 || V <= 0.0) { /* either U and V have opposite signs, */ /* or at least one of them is +/- 0.0 */ return U < V ? LESS_THAN : U > V ? GREATER_THAN : VERY_CLOSE; } /* now U and V are both strictly positive */ /* we are likely to need both scale factors, so */ (void)frexp(U, &e_u); /* U = F*2**e_u where 0 < F < 1 */ (void)frexp(V, &e_v); /* V = G*2**e_v where 0 < G < 1 */ if (U >= V) { /* there are three possibilities: U :>:, :~:, or :=: V */ /* since U >= V, e_u >= e_v, so */ /* max(e_u, e_v) = e_u, min(e_u, e_v) = e_v */ /* Calculating U-V may underflow, but as long as underflow */ /* is flushed to 0, that's ok. */ t = U-V; return t > ldexp(epsilon, e_u) ? GREATER_THAN : t > ldexp(epsilon, e_v) ? PRETTY_CLOSE : VERY_CLOSE; } else { t = V-U; return t > ldexp(epsilon, e_v) ? LESS_THAN : t > ldexp(epsilon, e_u) ? PRETTY_CLOSE : VERY_CLOSE; } } ------------------------- end of code inclusion -------------------------