Purpose of function is unknown, what does it do?

Question

Could someone please help me understand what this function is doing if the input is a complex number a+bi and real = a, imag = b

I have no idea what it could be computing but maybe I am missing something obvious?

double function(double real, double imag)
{
  double y;
  double a;
  double b;
  a = fabs(real);
  b = fabs(imag);
  if (a < b) 
  {
    a /= b;
    y = b * sqrt(a * a + 1.0);
  } 
  else if (a > b) 
  {
    b /= a;
    y = a * sqrt(b * b + 1.0);
  } 
  else if (b == NAN) 
  {
    y = b;
  } 
  else 
  {
    y = a * sqrt(2);
  }

  return y;
}

Answer 1

The code is a defective attempt to compute the magnitude (absolute value) of the complex number passed to it without incurring needless overflow.

Consider the complex number a + b i, where a and b are the values assigned to a and b in the first few lines of the function. Its magnitude is √(a²+b²). However, if a or b is large, the floating-point calculation a*a might overflow the finite range of the floating-point format and produce infinity (∞) even though the magnitude is within the range. As a simple example, let a be 2¹⁰⁰⁰ and b be 0. Then the magnitude is √(2²⁰⁰⁰+0) = 2¹⁰⁰⁰, but computing sqrt(a*a + b*b) yields infinity. (Since a*a overflowed and produced ∞, and the rest of the calculations then produce ∞ too.)

The code attempts to solve that by dividing the smaller of a and b by the larger and using a calculation that is mathematically equivalent but that does not overflow. For example, if a < b is true, it executes:

a /= b;
y = b * sqrt(a * a + 1.0);

Then a /= b produces a value less than 1, so all calculation prior to the last are safely within the floating-point finite range: a * a is less than 1, a * a + 1.0 is less than 2, and sqrt(a * a + 1.0) is less than 1.42. When we multiply by b, the final result might overflow to ∞. There are two reasons this might happen: The magnitude of a + b i might exceed the floating-point finite range, so the overflow is correct. Or rounding in the prior calculations might have caused sqrt(a * a + 1.0) to be slightly larger than the mathematical result and sufficient to cause b * sqrt(a * a + 1.0) to overflow when actual value of the magnitude is within the range. As this is not our focus, I will not analyze this case further in this answer.

Aside from that rounding issue, the first two cases are fine. However, this code is incorrect:

else if (b == NAN)

Per IEEE-754 2008 5.11 and IEEE-754 1985 5.7, a NaN is not less than, equal to, or greater than any operand, including itself. It is unordered. This means b == NAN must return false if IEEE-754 is used. C 2018 does not require IEEE-754 be used, but footnote 22 (at 5.2.4.2.2 4) says that, if IEC 60559:1989 (effectively IEEE-754) is not supported, the terms “quiet NaN” and “signaling NaN” in the C standard are intended to apply to encodings with similar behavior. And 7.12 5 tells us that NAN expands to a float representing a quiet NaN. Thus, in b == NAN, NAN should behave as an IEEE-754 NaN, and so b == NAN should yield 0, for false.

Therefore, this code governed by else if (b == NAN) is never executed:

y = b;

Instead, execution falls through to the else, which executes:

y = a * sqrt(2);

If a is a NaN, the result is a NaN, as desired. However, if a is a number and b is a NaN, this produces a number as a result when the desired result would be a NaN. Thus, the code is broken.