How should I tackle ARCTAN (and prevent deviding by zero) without the convenience of ARCTAN2(n,m) being available?

Question

I try to determine the angle from a point (n,m) to (0,0). Without arctan2 being available, I'm running into the problem that m can be 0, which leads to a possible division by zero.

What would be an elegant, correct solution to tackling this issue?

Answer 1

Define your version of arctan2. An example in C as a macro:

#define atan2(n,m)   (m)==0 ? M_PI_2 : atan((n)/(m))

Of course, you can elaborate for finding the quadrant depending on signs of n and m.

Answer 2

If atan2 is not available you have to check for the divide by zero condition and all the other special cases in your code. Easy as that. The wikipedia entry on atan2 has all the conditions you need.

If your target hardware supports divide by zero exceptions for floating point operations you have another option:

Install a low-level handler that checks the exception cause and if it happends to be a atan division fix the problem. This will make your atan2 faster if exceptions are rare, but it requires low level tinkering and is not portable.

Answer 3

I believe this is a correct implementation of atan2 using atan (doesn't handle infinities, though):

float my_atan2(float y, float x)
{
    if(x == 0) // might also want to use fabs(x) < 1e-6 or something like that
    {
        if(y > 0)
            return M_PI_2;
        else
            return -M_PI_2;
    }
    else if(x > 0)
    {
        return atan(y/x);
    }
    else 
    {
        // x < 0                                                                                
        if(y > 0)
            return M_PI + atan(y/x);
        else
            return -M_PI + atan(y/x);
    }
}

Test harness:

int main()
{
    for(int i = -360; i <= 360; i++)
    {
        float x = cos(i / 180.0 * M_PI);
        float y = sin(i / 180.0 * M_PI);


        float good = atan2(y, x);
        float mine = my_atan2(y, x);


        if(fabs(good - mine) > 1e-6)
        {
            printf("%d %f %f %f %f\n", i, x, y, good, mine);
        }
    }
}

Answer 4

don't use the conventional quadrants, use the branch points defined by the lines y = +/- x, and use the 1st two steps of a CORDIC-like algorithm (e.g. rotate coordinates by a known angle and keep track of how much you've rotated):

function atan2_substitute(x,y)
{
   double angle = 0;
   if (x < y)
   { angle = M_PI; x = -x; y = -y; }
   // this guarantees that the angle is between -135 and +45 degrees

   if (x < -y)
   {
     angle -= M_PI/2; tmp = x; x = -y; y = tmp;
   }
   // this guarantees that the angle is between -45 and +45

   angle += atan(y/x);

   if (angle > M_PI)
      angle -= 2*M_PI;
   // fails at 0,0; otherwise is accurate over the entire plane
}

the reason for doing it this way is that atan() may be more likely to be accurate for ratios y/x between -1 and +1, than for ratios greater than 1. (though a good atan() algorithm would recognize this and take the reciprocal)

Answer 5

Implement the standard arctan(n, m) using Taylor series and do the following before computing arctan:

if (m == 0) {
    if (n < 0) return Pi;
    return 0;
}

Another few tricks:

1) if |m| < |n|, swap m, n and then compute arctan. Finally subtract the result from Pi/2

2) if |m| is close to |n|, compute arctan of the half angle using the half angle formula arctan(x) = 2*arctan(x/(1+sqrt(1+x*x))), otherwise, for such values, arctan converges very slowly