How to compare two complex numbers?

Question

In C, complex numbers are float or double and have same problem as canonical types:

#include <stdio.h>
#include <complex.h>

int main(void)
{
    double complex a = 0 + I * 0;
    double complex b = 1 + I * 1;

    for (int i = 0; i < 10; i++) {
        a += .1 + I * .1;
    }

    if (a == b) {
        puts("Ok");
    }
    else {
        printf("Fail: %f + i%f != %f + i%f\n", creal(a), cimag(a), creal(b), cimag(b));
    }

    return 0;
}

The result:

$ clang main.c
$ ./a.out 
Fail: 1.000000 + i1.000000 != 1.000000 + i1.000000

I try this syntax:

a - b < DBL_EPSILON + I * DBL_EPSILON

But the compiler hate it:

main.c:24:15: error: invalid operands to binary expression ('_Complex double' and '_Complex double')
    if (a - b < DBL_EPSILON + I * DBL_EPSILON) {
        ~~~~~ ^ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

This last works fine but it’s a little fastidious:

fabs(creal(a) - creal(b)) < DBL_EPSILON && fabs(cimag(a) - cimag(b)) < DBL_EPSILON

Answer 1

Comparing 2 complex floating point numbers is much like comparing 2 real floating point numbers.

Comparing for exact equivalences often is insufficient as the numbers involved contain small computational errors.

So rather than if (a == b) code needs to be if (nearlyequal(a,b))

---

The usual approach is double diff = cabs(a - b) and then comparing diff to some small constant value like DBL_EPSILON.

This fails when a,b are large numbers as their difference may many orders of magnitude larger than DBL_EPSILON, even though a,b differ only by their least significant bit.

This fails for small numbers too as the difference between a,b may be relatively great, but many orders of magnitude smaller than DBL_EPSILON and so return true when the value are relatively quite different.

Complex numbers literally add another dimensional problem to the issue as the real and imaginary components themselves may be greatly different. Thus the best answer for nearlyequal(a,b) is highly dependent on the code's goals.

---

For simplicity, let us use the magnitude of the difference as compared to the average magnitude of a,b. A control constant ULP_N approximates the number of binary digits of least significance that a,b are allowed to differ.

#define ULP_N 4

bool nearlyequal(complex double a, complex double b) {
  double diff = cabs(a - b);
  double mag = (cabs(a) + cabs(b))/2;
  return diff <= (mag * DBL_EPSILON * (1ull << ULP_N));
}