Skip to content

Instantly share code, notes, and snippets.

@aliakseis
Last active December 30, 2022 10:28
Show Gist options
  • Save aliakseis/bdc9a4db2412b697ddcf4515bfd5853f to your computer and use it in GitHub Desktop.
Save aliakseis/bdc9a4db2412b697ddcf4515bfd5853f to your computer and use it in GitHub Desktop.
fast atan2
inline double fast_atan(double x) {
// http://www-labs.iro.umontreal.ca/~mignotte/IFT2425/Documents/EfficientApproximationArctgFunction.pdf.
//constexpr auto scaling_constant = 0.28086;
//constexpr auto N2 = 0.273;
constexpr auto scaling_constant = 0.274;
constexpr auto N2 = 0.258;
return ((M_PI_4 + N2) / 2. - (N2 / 2.) * std::abs(x)) * x
+ x / (2. + (scaling_constant * 2.) * x * x);
}
// The maximum absolute error across all intervals is less than 0.0011.
inline double fast_atan2(double y, double x) {
if (x != 0.0) {
if (fabs(x) > fabs(y)) {
auto z = y / x;
if (x > 0.0) {
// atan2(y,x) = atan(y/x) if x > 0
return fast_atan(z);
}
else if (y >= 0.0) {
// atan2(y,x) = atan(y/x) + PI if x < 0, y >= 0
return fast_atan(z) + M_PI;
}
else {
// atan2(y,x) = atan(y/x) - PI if x < 0, y < 0
return fast_atan(z) - M_PI;
}
}
else {
// Use property atan(y/x) = PI/2 - atan(x/y) if |y/x| > 1.
auto z = x / y;
if (y > 0.0) {
// atan2(y,x) = PI/2 - atan(x/y) if |y/x| > 1, y > 0
return -fast_atan(z) + M_PI / 2;
}
else {
// atan2(y,x) = -PI/2 - atan(x/y) if |y/x| > 1, y < 0
return -fast_atan(z) - M_PI / 2;
}
}
}
else
{
if (y > 0.0) {
// x = 0, y > 0
return M_PI / 2;
}
else if (y < 0.0) {
// x = 0, y < 0
return -M_PI / 2;
}
}
return 0.0; // x,y = 0. Could return NaN instead.
}
@aliakseis
Copy link
Author

The maximum absolute error across all intervals is less than 0.0011.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment