Improved fast atan2 and hypot approximations

gistfile1.txt

double fast_atan(double x) {
	return x / (1 + 0.3211 * x * x);
}

double do_fast_atan2(double y, double x, bool abs_x_lt_abs_y) {
	if (abs_x_lt_abs_y) {
		double z = x / y;
		double atan_z = fast_atan(z);

		if (y > 0.0) {
			return -atan_z + M_PI / 2;
		}
		else {
			return -atan_z - M_PI / 2;
		}
	}

	double z = y / x;
	double atan_z = fast_atan(z);

	if (x > 0.0) {
		return atan_z;
	}
	else if (y >= 0.0) {
		return atan_z + M_PI;
	}
	return atan_z - M_PI;
}

double do_fast_atan2_rotated(double y, double x) {
	if (fabs(x) < fabs(y)) {
		double z = x / y;
		double atan_z = fast_atan(z);

		if (y > 0.0) {
			return -atan_z + M_PI / 2 - M_PI / 4;
		}
		else {
			return -atan_z - M_PI / 2 - M_PI / 4;
		}
	}

	double z = y / x;
	double atan_z = fast_atan(z);

	if (x > 0.0) {
		return atan_z - M_PI / 4;
	}
	return atan_z + M_PI - M_PI / 4;
}

// The maximum absolute error across all intervals is less than 0.00012.
double fast_atan2(double y, double x) {
	const double TAN_PI_OVER_8 = 0.41421356237;

	if (x == 0.0) {
		return (y > 0.0) ? M_PI / 2 : ((y < 0.0) ? -M_PI / 2 : 0.0);
	}

	double abs_x = fabs(x);
	double abs_y = fabs(y);

	bool abs_x_lt_abs_y = abs_x < abs_y;

	if (abs_x_lt_abs_y ? abs_x > TAN_PI_OVER_8 * abs_y : abs_y > TAN_PI_OVER_8 * abs_x) {
		double x_ = x - y;
		double y_ = x + y;
		return do_fast_atan2_rotated(y_, x_);
	}

	return do_fast_atan2(y, x, abs_x_lt_abs_y);
}

double fast_hypot(double x, double y) {
    x = std::abs(x);
    y = std::abs(y);
    auto t = std::min(x, y);
    x = std::max(x, y);
    y = t;
    t = t / x;

    const double THRESHOLD = M_SQRT2 - 1.0;
    if (t > THRESHOLD)
    {
        auto x_ = x + y;
        auto y_ = x - y;

        t = y_ / x_;
        return x_ * (1 + t * t / 2) * M_SQRT1_2;
    }
    return x * (1 + t * t / 2);
}

Approximations of the atan2 and hypot functions are commonly found in real-time signal processing applications. Numerous algorithms are available to implement the atan2 and hypot functions. Most solutions are based on the approximation interval [−1, 1]. However, for arguments close to one, this is inefficient. It turns out that the approximation interval can simply be reduced to [-tan(Pi/8), tan(Pi/8)] by rotating the vector formed by X and Y by 45 degrees.