Best non-trigonometric floating point approximation of tanh(x) in 10 instructions or less

After doing much exploratory work, I came to the conclusion that approach 2 is the most promising direction. Since division is very fast on the asker's platform, rational approximations are attractive. The platform's support for FMA should be exploited aggressively. Below I am showing C code that implements a fast tanhf() in seven operations and achieves maximum absolute error of less than 3.2-3.

I used the Remez algorithm to compute the coefficients for the rational approximation and used a heuristic search to reduce these coefficients to as few bits as feasible, which may benefit some processor architectures that are able to incorporate floating-point data into an immediate field of commonly used floating-point instructions.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Fast computation of hyperbolic tangent. Rational approximation with clamping.
   Maximum absolute errror = 3.16989899e-3 @ +/-2.50056267
*/
float fast_tanhf_rat (float x)
{
    const float n0 = -8.70117188e-1f; // -0x1.bd8000p-1
    const float n1 = -8.92639160e-3f; // -0x1.248000p-7
    const float d0 =  2.72949219e+0f; //  0x1.5d6000p+1
    float x2 = x * x;
    float num = fmaf (n0, x2, n1);
    float den = x2 + d0;
    float quot = num / den;
    float res = fmaf (quot, x, x);
    res = fminf (fmaxf (res, -1.0f), 1.0f);
    return res;
}

int main (void)
{
    double ref, err, maxerr = 0;
    float arg, res, maxerrloc = INFINITY;
    maxerr = 0;
    arg = 0.0f;
    while (arg < 0x1.0p64f) {
        res = fast_tanhf_rat (arg);
        ref = tanh ((double)arg);
        err = fabs ((double)res - ref);
        if (err > maxerr) {
            maxerr = err;
            maxerrloc = arg;
        }
        arg = nextafterf (arg, INFINITY);
    }
    arg = -0.0f;
    while (arg > -0x1.0p64f) {
        res = fast_tanhf_rat (arg);
        ref = tanh ((double)arg);
        err = fabs ((double)res - ref);
        if (err > maxerr) {
            maxerr = err;
            maxerrloc = arg;
        }
        arg = nextafterf (arg, -INFINITY);
    }
    printf ("maximum absolute error = %15.8e @ %15.8e\n", maxerr, maxerrloc);
    return EXIT_SUCCESS;
}

Given that asker budgeted for up to ten operations, we can increase the degree of both numerator and denominator polynomials by one to achieve a fast tanhf() implementation comprising nine operations that has significantly lower maximum absolute error, less than 6.8e-5:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Fast computation of hyperbolic tangent. Rational approximation with clamping.
   Maximum absolute error = 6.74903450e-5 @ +/-5.14831972
 */
float fast_tanhf_rat2 (float x)
{
    const float n0 = -9.48013306e-001f; // -0x1.e56200p-1
    const float n1 = -2.61120605e+001f; // -0x1.a1cb00p+4
    const float n2 = -2.34234333e-003f; // -0x1.330400p-9
    const float d0 =  3.41279297e+001f; //  0x1.110600p+5
    const float d1 =  7.84062500e+001f; //  0x1.39a000p+6
    float x2 = x * x;
    float num = fmaf (fmaf (n0, x2, n1), x2, n2);
    float den = fmaf (x2 + d0, x2, d1);
    float quot = num / den;
    float res = fmaf (quot, x, x);
    res = fminf (fmaxf (res, -1.0f), 1.0f);
    return res;
}

int main (void)
{
    double ref, err, maxerr = 0;
    float arg, res, maxerrloc = INFINITY;
    maxerr = 0;
    arg = 0.0f;
    while (arg < 0x1.0p32f) {
        res = fast_tanhf_rat2 (arg);
        ref = tanh ((double)arg);
        err = fabs ((double)res - ref);
        if (err > maxerr) {
            maxerr = err;
            maxerrloc = arg;
        }
        arg = nextafterf (arg, INFINITY);
    }
    arg = -0.0f;
    while (arg > -0x1.0p32f) {
        res = fast_tanhf_rat2 (arg);
        ref = tanh ((double)arg);
        err = fabs ((double)res - ref);
        if (err > maxerr) {
            maxerr = err;
            maxerrloc = arg;
        }
        arg = nextafterf (arg, -INFINITY);
    }
    printf ("maximum absolute error = %15.8e @ %15.8e\n", maxerr, maxerrloc);
    return EXIT_SUCCESS;
}

Nic Schraudolph, author of the paper describing the exponential approximation that the previous version of this answer uses, suggests the following. It has error 0.5%.

Java implementation (for portable bit munging):

public class Tanh {
  private static final float m = (float)((1 << 23) / Math.log(2));
  private static final int b = Float.floatToRawIntBits(1);

  private static float tanh(float x) {
    int y = (int)(m * x);
    float exp_x = Float.intBitsToFloat(b + y);
    float exp_minus_x = Float.intBitsToFloat(b - y);
    return (exp_x - exp_minus_x) / (exp_x + exp_minus_x);
  }

  public static void main(String[] args) {
    double error = 0;
    int end = Float.floatToRawIntBits(10);
    for (int i = 0; i <= end; i++) {
      float x = Float.intBitsToFloat(i);
      error = Math.max(error, Math.abs(tanh(x) - Math.tanh(x)));
    }
    System.out.println(error);
  }
}