public static double tanh(double x) {
		boolean negate = false;

		if (x != x) {
			return x;
		}

		// tanh[z] = sinh[z] / cosh[z]
		// = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
		// = (exp(2x) - 1) / (exp(2x) + 1)

		// for magnitude > 20, sinh[z] == cosh[z] in double precision

		if (x > 20.0) {
			return 1.0;
		}

		if (x < -20) {
			return -1.0;
		}

		if (x == 0) {
			return x;
		}

		if (x < 0.0) {
			x = -x;
			negate = true;
		}

		double result;
		if (x >= 0.5) {
			double hiPrec[] = new double[2];
			// tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
			exp(x * 2.0, 0.0, hiPrec);

			double ya = hiPrec[0] + hiPrec[1];
			double yb = -(ya - hiPrec[0] - hiPrec[1]);

			/* Numerator */
			double na = -1.0 + ya;
			double nb = -(na + 1.0 - ya);
			double temp = na + yb;
			nb += -(temp - na - yb);
			na = temp;

			/* Denominator */
			double da = 1.0 + ya;
			double db = -(da - 1.0 - ya);
			temp = da + yb;
			db += -(temp - da - yb);
			da = temp;

			temp = da * HEX_40000000;
			double daa = da + temp - temp;
			double dab = da - daa;

			// ratio = na/da
			double ratio = na / da;
			temp = ratio * HEX_40000000;
			double ratioa = ratio + temp - temp;
			double ratiob = ratio - ratioa;

			// Correct for rounding in division
			ratiob += (na - daa * ratioa - daa * ratiob - dab * ratioa - dab
					* ratiob)
					/ da;

			// Account for nb
			ratiob += nb / da;
			// Account for db
			ratiob += -db * na / da / da;

			result = ratioa + ratiob;
		} else {
			double hiPrec[] = new double[2];
			// tanh(x) = expm1(2x) / (expm1(2x) + 2)
			expm1(x * 2.0, hiPrec);

			double ya = hiPrec[0] + hiPrec[1];
			double yb = -(ya - hiPrec[0] - hiPrec[1]);

			/* Numerator */
			double na = ya;
			double nb = yb;

			/* Denominator */
			double da = 2.0 + ya;
			double db = -(da - 2.0 - ya);
			double temp = da + yb;
			db += -(temp - da - yb);
			da = temp;

			temp = da * HEX_40000000;
			double daa = da + temp - temp;
			double dab = da - daa;

			// ratio = na/da
			double ratio = na / da;
			temp = ratio * HEX_40000000;
			double ratioa = ratio + temp - temp;
			double ratiob = ratio - ratioa;

			// Correct for rounding in division
			ratiob += (na - daa * ratioa - daa * ratiob - dab * ratioa - dab
					* ratiob)
					/ da;

			// Account for nb
			ratiob += nb / da;
			// Account for db
			ratiob += -db * na / da / da;

			result = ratioa + ratiob;
		}

		if (negate) {
			result = -result;
		}

		return result;
	}