🚀 go-pugleaf

jtomich@IntNet.net (Jeff Tomich) wrote:

>Having a hard time to figure out a function on how to compare doubles. 
>Any ideas?

After a couple of hints in email, and looking at some of the other
posted solutions (and posted problems) I've come up with the following,
which should avoid most of the problems cited so far.  It should work
correctly over a relatively wide range, while hopefully retaining
reasonable speed:

/* fcomp.h */
#ifndef FCOMP_H_INCLUDED
#define FCOMP_H_INCLUDED

#include <math.h>
#include "dist.h"

int fcomp(double a, double b);

#endif

/* fcomp.c */
#include "fcomp.h"

/* This asks for equality to roughly 10 digits.
 */
static const double delta = 1e-20;


/* Iterations   Accuracy
 *  2          6.5 digits
 *  3           20 digits
 *  4           62 digits
 * assuming a numeric type able to maintain that degree of accuracy in
 * the individual operations.
 */
#define ITER 3

static double dist(double P, double Q) {
/* A reasonably robust method of calculating `sqrt(P*P + Q*Q)'
 *
 * Transliterated from _More Programming Pearls, Confessions of a Coder_
 * by Jon Bentley, pg. 156. (copr. Bell Telephone Labs, 1988, published
 * Addison-Wesley.  ISBN 0-201-11889-0 )
 */

    double R;
    int i;

    P = fabs(P);
    Q = fabs(Q);

    if (P<Q) {
        R = P;
        P = Q;
        Q = R;
    }

/* The book has this as:
 *  if P = 0.0 return Q; # in AWK
 * However, this makes no sense to me - we've just insured that P>=Q, so
 * P==0 only if Q==0;  OTOH, if Q==0, then distance == P...
 */
    if ( Q == 0.0 )
        return P;

    for (i=0;i<ITER;i++) {
        R = Q / P;
        R + R * R;
        R = R / (4.0 + R);
        P = P + 2.0 * R * P;
        Q = Q * R;
    }
    return P;
}

int fcomp(double a, double b) {
/* Hopefully a reasonably robust comparison of a and b.
 * Returns 1 if they are equal to the precision indicated by delta.
 * Returns 0 if the two differ by more than delta.
 * As far as I can tell, the major possibility of overflow is if a and b
 * are both within a factor of 2 of DBL_MAX in magnitude, and differ in
 * sign.  
 * Delta is small enought that b+delta should never overflow: if  b is 
 * is close to DBL_MAX, adding delta will generally have no effect.
 */

    double difference = fabs(a-b);

    double distance = dist(a,b+delta);

    if ( difference / distance < delta )
        return 1;
    return 0;
}

I'm _NOT_ particular skilled at numeric programming, so anybody who is
is more than welcome to tear it up if they see fit - this is basically
just one of the other posts recast to use Moler & Morrison's distance
finding algorithm...
    Later,
    Jerry.

/* I can barely express my own opinions; I certainly can't
 * express anybody else's.
 *
 * The universe is a figment of its own imagination.
 */