On Feb 16, 9:09 pm, Jerry Coffin <jcof...@[EMAIL PROTECTED]
> wrote:
> In article <154dedac-ee33-44d1-bb89-
> b12891145...@[EMAIL PROTECTED]
>, james.ka...@[EMAIL PROTECTED]
says...
> [ ... ]
> > You really do have to calculate the espilon as a function of the
> > two values involved, something like:
> That's what I was using frexp and ldexp to do. I didn't take all the
> corner cases into account, but the idea was definitely to adjust the
> magnitude of the precision to fit the magnitudes of the numbers.
I realized that. I was just wondering why you didn't take the
standard approach. Or what I believe is the standard
approach---I'm really not that up in numeric programming. As I
pointed out, the standard approach has problems with overflow
and underflow, but depending on the application, that will often
not be an issue.
[...]
> > Except, of course, that still has problems with very big values
> > (for which a + b might overflow) and very small values (for
> > which a - b might underflow, and result in zero). And I'm not
> > really sure what you should do for values which are near the
> > minimum representable: Do they compare approximately equal to
> > zero? Do they compare equal even if their signs are not the
> > same?
> See below -- I'm pretty sure my new code avoids problems with
> either underflow or overflow.
> I'm not sure the problems with numbers extremely close to zero are
> entirely real. The question being asked is whether two numbers are equal
> to a specified precision. The answer to that is a simple yes/no. That
> question may not always be the appropriate one for every possible
> calculation, but that doesn't mean there's anything wrong with the code
> as it stands.
[...]
> This is a basic question of what should be specified, and whether that
> specification fits all possible purposes.
[and later...]
> I wouldn't attempt to claim that it fits every possible
> situation,[...]
I think that that is really the crux of the problem. More
im****tant than the actual code is understanding exactly what it
is doing, and being able to analyse when it is appropriate, and
when not.
> The code in the FAQ might fit
> some specification, but only a very limited one (i.e. that all the
> numbers involved must be quite close to 1.0 or -1.0).
As I said, I haven't seen it. If it is something like "return
fabs(a-b) < prec", then I agree that it's almost never
appropriate.
Except, of course, when it is. In process control software,
I've had cases where the requirement was to keep the actual
value within a specified epsilon of the target value. In such
cases, of course, the "tolerance" is usually something between
1% and 10%, and the actual values will be of similar magnitude.
On the other hand, you rarely test for equality in such cases,
since the action you take will depend on whether the actual
value is less than or greater than the target value. (Is
"target value" the correct English word here: "valeur de
consigne" in French, "Sollwert" in German?)
> I wouldn't attempt to claim that it fits every possible
> situation, but I think the following code is (at least
> starting to get close to) correct. The specification I'm
> following is that it is supposed to determine whether two
> numbers are equal within a specified relative difference.
> Here's the code, along with a few test cases:
> #include <math.h>
> bool isNearlyEqual(double a, double b, double prec =3D 1e-6) {
> int mag1, mag2;
> double num1 =3D frexp(a, &mag1);
> double num2 =3D frexp(b, &mag2);
> if (abs(mag1-mag2) > 1)
> return false;
> num2 =3D ldexp(num2, mag2-mag1);
> return fabs(num2-num1) < prec;
> }
I'm still curious as to why you are trying to use frexp and
ldexp, rather than a more classical approach. Whatever that is:
I've been searching the Web for something, but none of my usual
numeric resources pages seem to turn up anything---should I
conclude that people competent in numeric processing don't have
such a function in their toolbox, and that thus, we are probably
wrong it trying to define one.
(FWIW: The reason I don't think that ldexp and frexp would be in
a classical solution, is that the classical solutions
doubtlessly date from the days of Fortran.)
[...]
> As an aside, I'd note that I don't really like the 'if
> (abs(mag2-mag1)> 1)' part, but I can't see anything a lot
> better at the moment.
I think it's a common characteristic for numerical functions to
special case limit cases. On the other hand, you definitely
need to test border values here.
> In theory, I'd like it to take the tolerance into account --
> e.g. if somebody specifies a tolerance of 100, it should allow
> numbers that are different by a factor of 100 to yield true.
> This would add some complexity, however, and I can't think of
> any real uses to justify it.
Agreed. I suspect as well that you'll get some errors when the
limit is around 2---your test will return false when when one
number is 1.9999 times the other, or some such. I'd say that
your missing an essential first line in your code:
assert( prec > 0.0 && prec < 1.0 ) ;
That seems like a reasonable constraint. (The first condition
should actually be something like "prec >=3D pow( 2.0,
-std::numeric_limits< double >::digits )" if the function is to
make sense. But maybe it is more reasonable to use prec >=3D 0.0,
with the do***entation saying that if prec is too small, the
function just becomes a very expensive way to write =3D=3D.)
--
James Kanze (GABI Software) email:james.kanze@[EMAIL PROTECTED]
en informatique orient=E9e objet/
Beratung in objektorientierter Datenverarbeitung
9 place S=E9mard, 78210 St.-Cyr-l'=C9cole, France, +33 (0)1 30 23 00 34
|
