🚀 go-pugleaf

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 -------------------------