atan2() is not implemented: it raises unavailable_feature_error

[Jordan08]

gaol::atan2(const interval& Y, const interval& X) is declared in gaol/gaol_interval.h, and expressions can call it (interval("atan2(1, 1)"), see gaol/gaol_interval_parser.ypp). But it is not implemented: it raises unavailable_feature_error ("atan2 not yet implemented") and returns the empty set (gaol_interval.cpp).

Found while writing tests/elementary.cpp, which leaves atan2() out for this reason.

A possible implementation: atan() of the quotients at the corners of [X]×[Y], quadrant by quadrant, with the cut along the negative x axis, as IBEX implements it for its intervals.

[Jordan08]

Implemented by 1bcda57, as the atan2 of IEEE 1788-2015 (Table 9.1): defined on the plane but (0, 0), with values in (−π, π].

• Algorithm. Not the arctangent of the quotients, whose division adds a rounding: mathlib has uatan2(), correctly rounded. In each quadrant the angle is monotonic in y and in x, so its least and greatest values over the box Y × X are at two corners, which the signs of the bounds give; uatan2() is taken there and moved one double outward. With CRlibm, which has no atan2, and with the math library of the system, the atan2() of the libm is taken, moved outward as the hyperbolic functions are.
• The cut. The angle jumps from π to −π across the half-line y = 0, x < 0. A box with points on that half-line, whose angle is π, and points below it, whose angles are next to −π, gives [-pi, pi], the hull of its angles: atan2([-1, 0], [-2, -1]) and atan2([-1, 1], [-2, -1]) are [-pi, pi], while atan2([0, 1], [-2, -1]) is [3pi/4, pi].
• Special cases. atan2([0], [0]) is empty; a corner on an axis, on a diagonal of the right half-plane or at an infinity gives 0, ±π/4, ±π/2 or π as the tightest bounds of GAOL's constants: atan2([0], [-2, -1]) is the tightest enclosure of π, atan2([1, +oo], [1, +oo]) is [0, pi/2], atan2([0], [-1, 1]) is [0, pi].
• In expressions. interval("4*atan2(1, 1)") is the tightest enclosure of π. The expressions the tests used for an exception going through the parser (1+pow(2, atan2(1,1))) are now computed.

Accuracy (doc/accuracy.md): within one double of the tightest bounds; tightest where a bound is 0, ±π/4, ±π/2 or ±π.

Tests.

• tests/elementary.cpp: 87 points of the four quadrants, next to the axes and with ratios beyond the doubles, and 324 boxes in every position about the axes and the cut, against values computed by mpmath with 2000 bits (the hull being taken over the corners and the points where the box meets the axes, independently of GAOL's choice of corners): all within one double. 25 boxes with infinite bounds or on an axis: the tightest bounds.
• tests/rounding_direction.cpp: the same results whatever the rounding direction of the caller, left upward afterwards.
• doc/compare: 29 cases added (263 to 291). GAOL gives the result of IEEE 1788-2015 in 24 and one double more in 5; libieeep1788 gives it in the 29; filib++ has no atan2; Solaris Studio continues the angle across the cut (atan2([-1, 1], [-2, -1]) is [3pi/4, 5pi/4]) and gives [-pi, pi] for atan2([0], [0]).
• The 110 jobs of the continuous integration pass (Linux, Windows).

Time: 100 ns for a box within a quadrant on an Intel i7-1185G7 (GCC 9.4), against 68 ns for atan().

manual/gaol.tex describes it since 42a8114 (the PDF is not rebuilt: #13). Closing.