<p>I pondered the various ideas I received in comments and also ran a few experiments based on that feedback. In the end I decided that a refined heuristic search was the best way forward. I have now managed to reduce the maximum error for <code>atanf_poly()</code> to 1.01036 ulps, with just three arguments exceeding my stated goal of a 1 ulp error bound:</p>

<pre><code>ulp = -1.00829 @ |a| =  9.80738342e-001 0x1.f62356p-1 (3f7b11ab)
ulp = -1.01036 @ |a| =  9.87551928e-001 0x1.f9a068p-1 (3f7cd034)
ulp =  1.00050 @ |a| =  9.99375939e-001 0x1.ffae34p-1 (3f7fd71a)
</code></pre>

<p>Based on the manner of generating the improved approximation there is no guarantee that this is a best approximation; no scientific breakthrough here. As the ulp error of the current solution is not yet perfectly balanced, and since continuing the search continues to deliver better approximations (albeit at exponentially increasing time intervals) my <em>guess</em> is that a 1 ulp error bound is achievable, but at the same time we seem to be very close to the best machine-optimized approximation already.</p>

<p>The better quality of the new approximation is the result of a refined search process. I observed that all of the largest ulp errors in the polynomial occur close to unity, say in [0.75,1.0] to be conservative. This allows for a fast scan for interesting coefficient sets whose maximum error is smaller than some bound, say 1.08 ulps. I can then test in detail and exhaustively all coefficient sets within a heuristically chosen hyper-cone anchored at that point. This second step searches for minimum ulp error as the primary goal, and maximum percentage of correctly rounded results as a secondary objective. By using this two-step process across all four cores of my CPU I was able to significantly speed up the search process: I have been able to check about 2<sup>21</sup> coefficient sets so far. </p>

<p>Based on the range of each coefficient across all "close" solutions I now estimate that the total useful search space for this approximation problem is >= 2<sup>24</sup> coefficient sets rather than the more optimistic number of 2<sup>20</sup> I threw out before. This seems like a feasible problem to solve for someone who is either very patient or has lots of computational horse-power at their disposal.</p>

<p>My updated code is as follows:</p>

<pre><code>// max ulp err = 1.01036
float atanf_poly (float a)
{
    float r, s;
    s = a * a;
    r =              0x1.7ed22cp-9f;
    r = fmaf (r, s, -0x1.0c2c2ep-6f);
    r = fmaf (r, s,  0x1.61fdf6p-5f);
    r = fmaf (r, s, -0x1.3556b4p-4f);
    r = fmaf (r, s,  0x1.b4e12ep-4f);
    r = fmaf (r, s, -0x1.230ae0p-3f);
    r = fmaf (r, s,  0x1.9978eep-3f);
    r = fmaf (r, s, -0x1.5554dap-2f);
    r = r * s;
    r = fmaf (r, a, a);
    return r;
}

// max ulp err = 1.51871
float my_atanf (float a)
{
    float r, t;
    t = fabsf (a);
    r = t;
    if (t &gt; 1.0f) {
        r = 1.0f / r;
    }
    r = atanf_poly (r);
    if (t &gt; 1.0f) {
        r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, -r); // pi/2 - r
    }
    r = copysignf (r, a);
    return r;
}
</code></pre>

<hr>

<p><strong>Update</strong> (after revisiting the issue two-and-a-half years later)</p>

<p>Using T. Myklebust's <a href="https://arxiv.org/abs/1508.03211" rel="nofollow noreferrer">draft publication</a> as a starting point, I found the arctangent approximation on [-1,1] that has the smallest error to have a maximum error of 0.94528 ulp.</p>

<pre><code>/* Based on: Tor Myklebust, "Computing accurate Horner form approximations 
   to special functions in finite precision arithmetic", arXiv:1508.03211,
   August 2015. maximum ulp err = 0.94528
*/
float atanf_poly (float a)
{
    float r, s;
    s = a * a;                        
    r =              0x1.6d2086p-9f;  //  2.78569828e-3
    r = fmaf (r, s, -0x1.03f2ecp-6f); // -1.58660226e-2
    r = fmaf (r, s,  0x1.5beebap-5f); //  4.24722321e-2
    r = fmaf (r, s, -0x1.33194ep-4f); // -7.49753043e-2
    r = fmaf (r, s,  0x1.b403a8p-4f); //  1.06448799e-1
    r = fmaf (r, s, -0x1.22f5c2p-3f); // -1.42070308e-1
    r = fmaf (r, s,  0x1.997748p-3f); //  1.99934542e-1
    r = fmaf (r, s, -0x1.5554d8p-2f); // -3.33331466e-1
    r = r * s;
    r = fmaf (r, a, a);
    return r;
}
</code></pre>
Best machine-optimized polynomial minimax approximation to arctangent on [-1,1]?
<p>For the simple and efficient implementation of fast math functions with reasonable accuracy, polynomial minimax approximations are often the method of choice. Minimax approximations are typically generated with a variant of the Remez algorithm. Various widely available tools such as Maple and Mathematica have built-in functionality for this. The generated coefficients are typically computed using high-precision arithmetic. It is well-known that simply rounding those coefficients to machine precision leads to suboptimal accuracy in the resulting implementation.</p>

<p>Instead, one searches for closely related sets of coefficients that are exactly representable as machine numbers to generate a machine-optimized approximation. Two relevant papers are:</p>

<p><em>Nicolas Brisebarre, Jean-Michel Muller, and Arnaud Tisserand, "Computing Machine-Efficient Polynomial Approximations", ACM Transactions on Mathematical Software, Vol. 32, No. 2, June 2006, pp. 236–256.</em></p>

<p><em>Nicolas Brisebarre and Sylvain Chevillard, "Efficient polynomial L∞-approximations", 18th IEEE Symposium on Computer Arithmetic (ARITH-18), Montpellier (France), June 2007, pp. 169-176.</em></p>

<p>An implementation of the LLL-algorithm from the latter paper is available as the <code>fpminimax()</code> command of the <a href="http://sollya.gforge.inria.fr/sollya-3.0/help.php?name=fpminimax">Sollya tool</a>. It is my understanding that all algorithms proposed for the generation of machine-optimized approximations are based on heuristics, and that it is therefore generally unknown what accuracy can be achieved by an optimal approximation. It is not clear to me whether the availability of FMA (fused multiply-add) for the evaluation of the approximation has an influence on the answer to that question. It seems to me naively that it should.</p>

<p>I am currently looking at a simple polynomial approximation for arctangent on [-1,1] that is evaluated in IEEE-754 single-precision arithmetic, using the Horner scheme and FMA. See function <code>atan_poly()</code> in the C99 code below. For lack of access to a Linux machine at the moment, I did not use Sollya to generate these coefficients, but used my own heuristic that could be loosely described as a mixture of steepest decent and simulated annealing (to avoid getting stuck on local minima). The maximum error of my machine-optimized polynomial is very close to 1 ulp, but ideally I would like the maximum ulp error to be below 1 ulp. </p>

<p>I am aware that I could change my computation to increase the accuracy, for example by using a leading coefficient represented to more than single-precision precision, but I would like to keep the code exactly as is (that is, as simple as possible) adjusting only the coefficients to deliver the most accurate result possible.</p>

<p>A "proven" optimal set of coefficients would be ideal, pointers to relevant literature are welcome. I did a literature search but could not find any paper that advances the state of the art meaningfully beyond Sollya's <code>fpminimax()</code>, and none that examine the role of FMA (if any) in this issue.</p>

<pre><code>// max ulp err = 1.03143
float atan_poly (float a)
{
    float r, s;
    s = a * a;
    r =              0x1.7ed1ccp-9f;
    r = fmaf (r, s, -0x1.0c2c08p-6f);
    r = fmaf (r, s,  0x1.61fdd0p-5f);
    r = fmaf (r, s, -0x1.3556b2p-4f);
    r = fmaf (r, s,  0x1.b4e128p-4f);
    r = fmaf (r, s, -0x1.230ad2p-3f);
    r = fmaf (r, s,  0x1.9978ecp-3f);
    r = fmaf (r, s, -0x1.5554dcp-2f);
    r = r * s;
    r = fmaf (r, a, a);
    return r;
}

// max ulp err = 1.52637
float my_atanf (float a)
{
    float r, t;
    t = fabsf (a);
    r = t;
    if (t &gt; 1.0f) {
        r = 1.0f / r;
    }
    r = atan_poly (r);
    if (t &gt; 1.0f) {
        r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, -r); // pi/2 - r
    }
    r = copysignf (r, a);
    return r;
}
</code></pre>