This is a deep question without a simple answer that suits all cases. The wikipedia article is a brief introduction: https://en.wikipedia.org/wiki/Numerical_analysis. Especially the sections on error propagation and on numerical stability. 

See also https://stackoverflow.com/questions/13940316/floating-point-comparison-revisited.

In short the appropriate precision depends on the number of floating point operations done before comparison, the numerical stability of the algorithm being used, and whether the values being compared are truly on a continuum or can be considered rational such that the values to be compared can be "clamped" to the nearest rational values to stop errors from propagating. 

This means that in your comparison code the precision should vary by the sizes of the values being compared AND by case. The number you use is way too big for most values. The logic should also differ for ==, >, >=, <=, and <.

My beginner's stab at this which has worked so far is in: https://github.com/gogins/csound-ac/blob/master/CsoundAC/ChordSpaceBase.hpp.

On Mon, Feb 6, 2023, 19:34 Hlöðver Sigurðsson <hlolli@gmail.com> wrote:

Currently in the csound codebase we have this

#define RELATN(OPNAME,OP)                                \
  int32_t OPNAME(CSOUND *csound, RELAT *p)               \
  {   IGN(csound); *p->rbool = (*p->a OP *p->b) ? 1 : 0; \
       return OK; }

RELATN(eq,==)

And it seems to work fine until the numbers I'm comparing are within a deep recursion.

I tried to see how these two numbers are getting different from one another with 100 decimal precision

int32_t eq(CSOUND *csound, RELAT *p) {
  IGN(csound);
  *p->rbool = FL(*p->a) == FL(*p->b) ? 1 : 0;
  printf("*p->a(%.100f) == *p->b(%.100f) => *p->rbool(%d)\n", *p->a, *p->b, *p->rbool);
  return OK;
}

and it prints literally
*p->a(0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000) == *p->b(0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000) => *p->rbool(0)

so *p->rbool = 0

I read from some sources only that one shouldn't really compare two floating point values this way.

I tested with replacing it with

int compare_float(MYFLT f1, MYFLT f2) {
   MYFLT precision = 0.00001;
   if (((f1 - precision) < f2) &&
       ((f1 + precision) > f2)) {
     return 1;
   } else {
    return 0;
   }
 }

which fixes this. But now I worry about the precision required to make this an acceptable solution. I wonder why this hasn't come up before. Could it be that my two zeroes there are in fact different? It's strange that this works just fine until it doesn't.

Yes I stumbled on similar stackoverflow with heated discussion on this topic. I find it strange that after all these years, this hasn't been an issue. You mention "precision depends on the number of floating point operations done before comparison" I think this may explain the cause of me seeing this.

> This means that in your comparison code the precision should vary by the sizes of the values being compared AND by case. 

by this logic you would say that the equality operation as it is currently in Csound is incorrectly implemented?

On Tue, 7 Feb 2023 at 03:11, Michael Gogins <michael.gogins@gmail.com> wrote:

This is a deep question without a simple answer that suits all cases. The wikipedia article is a brief introduction: https://en.wikipedia.org/wiki/Numerical_analysis. Especially the sections on error propagation and on numerical stability. 

See also https://stackoverflow.com/questions/13940316/floating-point-comparison-revisited.

In short the appropriate precision depends on the number of floating point operations done before comparison, the numerical stability of the algorithm being used, and whether the values being compared are truly on a continuum or can be considered rational such that the values to be compared can be "clamped" to the nearest rational values to stop errors from propagating. 

This means that in your comparison code the precision should vary by the sizes of the values being compared AND by case. The number you use is way too big for most values. The logic should also differ for ==, >, >=, <=, and <.

My beginner's stab at this which has worked so far is in: https://github.com/gogins/csound-ac/blob/master/CsoundAC/ChordSpaceBase.hpp.

On Mon, Feb 6, 2023, 19:34 Hlöðver Sigurðsson <hlolli@gmail.com> wrote:

Currently in the csound codebase we have this

#define RELATN(OPNAME,OP)                                \
  int32_t OPNAME(CSOUND *csound, RELAT *p)               \
  {   IGN(csound); *p->rbool = (*p->a OP *p->b) ? 1 : 0; \
       return OK; }

RELATN(eq,==)

And it seems to work fine until the numbers I'm comparing are within a deep recursion.

I tried to see how these two numbers are getting different from one another with 100 decimal precision

int32_t eq(CSOUND *csound, RELAT *p) {
  IGN(csound);
  *p->rbool = FL(*p->a) == FL(*p->b) ? 1 : 0;
  printf("*p->a(%.100f) == *p->b(%.100f) => *p->rbool(%d)\n", *p->a, *p->b, *p->rbool);
  return OK;
}

and it prints literally
*p->a(0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000) == *p->b(0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000) => *p->rbool(0)

so *p->rbool = 0

I read from some sources only that one shouldn't really compare two floating point values this way.

I tested with replacing it with

int compare_float(MYFLT f1, MYFLT f2) {
   MYFLT precision = 0.00001;
   if (((f1 - precision) < f2) &&
       ((f1 + precision) > f2)) {
     return 1;
   } else {
    return 0;
   }
 }

which fixes this. But now I worry about the precision required to make this an acceptable solution. I wonder why this hasn't come up before. Could it be that my two zeroes there are in fact different? It's strange that this works just fine until it doesn't.