Is there a C/C++ function to safely handle division by zero?

Question

We have a situation we want to do a sort of weighted average of two values w1 & w2, based on how far two other values v1 & v2 are away from zero... for example:

• If v1 is zero, it doesn't get weighted at all so we return w2
• If v2 is zero, it doesn't get weighted at all so we return w1
• If both values are equally far from zero, we do a mean average and return (w1 + w2 )/2

I've inherited code like:

float calcWeightedAverage(v1,v2,w1,w2)
{
  v1=fabs(v1);
  v2=fabs(v2);
  return (v1/(v1+v2))*w1 + (v2/(v1+v2)*w2);
}

For a bit of background, v1 & v2 represent how far two different knobs are turned, the weighting of their individual resultant effects only depends how much they are turned, not in which direction.

Clearly, this has a problem when v1==v2==0, since we end up with return (0/0)*w1 + (0/0)*w2 and you can't do 0/0. Putting a special test in for v1==v2==0 sounds horrible mathematically, even if it wasn't bad practice with floating-point numbers.

So I wondered if

• there was a standard library function to handle this
• there's a neater mathematical representation

Answer 1

If the denominator is zero, how do you want it to default? You can do something like this:

static inline float divide_default(float numerator, float denominator, float default) {
    return (denominator == 0) ? default : (numerator / denominator);
}

float calcWeightedAverage(v1, v2, w1, w2)
{
  v1 = fabs(v1);
  v2 = fabs(v2);
  return w1 * divide_default(v1, v1 + v2, 0.0) + w2 * divide_default(v2, v1 + v2, 0.0);
}

Note that the function definition and use of static inline should really let the compiler know that it can inline.

Answer 2

This is the best I could come up with quickly

float calcWeightedAverage(float v1,float v2,float w1,float w2)
{
    float a1 = 0.0;
    float a2 = 0.0;

    if (v1 != 0)
    { 
        a1 = v1/(v1+v2) * w1;
    }

    if (v2 != 0)
    { 
        a2 = v2/(v1+v2) * w2;
    }

    return a1 + a2;
}

Answer 3

I don't see what would be wrong with just doing this:

float calcWeightedAverage( float v1, float v2, float w1, float w2 ) {
    static const float eps = FLT_MIN; //Or some other suitably small value.
    v1 = fabs( v1 );
    v2 = fabs( v2 );

    if( v1 + v2 < eps )
        return (w1+w2)/2.0f;
    else
        return (v1/(v1+v2))*w1 + (v2/(v1+v2)*w2);
}

Sure, no "fancy" stuff to figure out your division, but why make it harder than it has to be?

Answer 4

I'd do like this:

float calcWeightedAverage(double v1, double v2, double w1, double w2)
{
  v1 = fabs(v1);
  v2 = fabs(v2);
  /* if both values are equally far from 0 */
  if (fabs(v1 - v2) < 0.000000001) return (w1 + w2) / 2;
  return (v1*w1 + v2*w2) / (v1 + v2);
}

Answer 5

Personally I don't see anything wrong with an explicit check for divide by zero. We all do them, so it could be argued that not having it is uglier.

However, it is possible to turn off the IEEE divide by zero exceptions. How you do this depends on your platform. I know on windows it has to be done process-wide, so you can inadvertantly mess with other threads (and they with you) by doing it if you aren't careful.

However, if you do that your result value will be NaN, not 0. I highly dooubt that's what you want. If you are going to have to put a special check in there anyway with different logic when you get NaN, you might as well just check for 0 in the denominator up front.

Answer 6

Actually divide by zero is almost never going to be your problem. Floating point numbers are almost never zero. However, they can get very close to zero which causes numerical overflow. The best way to handle this is:

const float epsilon = //some value. or you could go with a #DEFINE epsilon if you prefer

bool nearEpsilon(float _value){
    retnr fabs( _value ) < epsilon;
}

so that you can circumvent the problems and return

numeric_limits<float>::max;

taken from #include <limits>.

Answer 7

This should work

#include <float.h>

float calcWeightedAverage(v1,v2,w1,w2)
{
  v1=fabs(v1);
  v2=fabs(v2);
  return (v1/(v1+v2+FLT_EPSILON))*w1 + (v2/(v1+v2+FLT_EPSILON)*w2);
}

edit: I saw there may be problems with some precision so instead of using FLT_EPSILON use DBL_EPSILON for accurate results (I guess you will return a float value).

Answer 8

So with a weighted average, you need to look at the special case where both are zero. In that case you want to treat it as 0.5 * w1 + 0.5 * w2, right? How about this?

float calcWeightedAverage(float v1,float v2,float w1,float w2)
{
  v1=fabs(v1);
  v2=fabs(v2);
  if (v1 == v2) {
    v1 = 0.5;
  } else {
    v1 = v1 / (v1 + v2); // v1 is between 0 and 1
  }
  v2 = 1 - v1; // avoid addition and division because they should add to 1      

  return v1 * w1 + v2 * w2;
}

Answer 9

You chould test for fabs(v1)+fabs(v2)==0 (this seems to be the fastest given that you've already computed them), and return whatever value makes sense in this case (w1+w2/2?). Otherwise, keep the code as-is.

However, I suspect the algorithm itself is broken if v1==v2==0 is possible. This kind of numerical instability when the knobs are "near 0" hardly seems desirable.

If the behavior actually is right and you want to avoid special-cases, you could add the minimum positive floating point value of the given type to v1 and v2 after taking their absolute values. (Note that DBL_MIN and friends are not the correct value because they're the minimum normalized values; you need the minimum of all positive values, including subnormals.) This will have no effect unless they're already extremely small; the additions will just yield v1 and v2 in the usual case.

Answer 10

You're trying to implement this mathematical function:

F(x, y) = (W1 * |x| + W2 * |y|) / (|x| + |y|)

This function is discontinuous at the point x = 0, y = 0. Unfortunately, as R. stated in a comment, the discontinuity is not removable - there is no sensible value to use at this point.

This is because the "sensible value" changes depending on the path you take to get to x = 0, y = 0. For example, consider following the path F(0, r) from r = R1 to r = 0 (this is equivalent to having the X knob at zero, and smoothly adjusting the Y knob down from R1 to 0). The value of F(x, y) will be constant at W2 until you get to the discontinuity.

Now consider following F(r, 0) (keeping the Y knob at zero and adjusting the X knob smoothly down to zero) - the output will be constant at W1 until you get to the discontinuity.

Now consider following F(r, r) (keeping both knobs at the same value, and adjusting them down simulatneously to zero). The output here will be constant at W1 + W2 / 2 until you go to the discontinuity.

This implies that any value between W1 and W2 is equally valid as the output at x = 0, y = 0. There's no sensible way to choose between them. (And further, always choosing 0 as the output is completely wrong - the output is otherwise bounded to be on the interval W1..W2 (ie, for any path you approach the discontinuity along, the limit of F() is always within that interval), and 0 might not even lie in this interval!)

---

You can "fix" the problem by adjusting the function slightly - add a constant (eg 1.0) to both v1 and v2 after the fabs(). This will make it so that the minimum contribution of each knob can't be zero - just "close to zero" (the constant defines how close).

It may be tempting to define this constant as "a very small number", but that will just cause the output to change wildly as the knobs are manipulated close to their zero points, which is probably undesirable.

Answer 11

The problem with using an explicit check for zero is that you can end up with discontinuities in behaviour unless you are careful as outlined in cafs response ( and if its in the core of your algorithm the if can be expensive - but dont care about that until you measure...)

I tend to use something that just smooths out the weighting near zero instead.

float calcWeightedAverage(v1,v2,w1,w2)
{
  eps = 1e-7; // Or whatever you like...
  v1=fabs(v1)+eps;
  v2=fabs(v2)+eps;
  return (v1/(v1+v2))*w1 + (v2/(v1+v2)*w2);
}

Your function is now smooth, with no asymptotes or division by zero, and so long as one of v1 or v2 is above 1e-7 by a significant amount it will be indistinguishable from a "real" weighted average.