public static double log1p(final double x) {

		if (x == -1) {
			return x / 0.0; // -Infinity
		}

		if (x > 0 && 1 / x == 0) { // x = Infinity
			return x;
		}

		if (x > 1e-6 || x < -1e-6) {
			double xpa = 1.0 + x;
			double xpb = -(xpa - 1.0 - x);

			double hiPrec[] = new double[2];

			final double lores = log(xpa, hiPrec);
			if (Double.isInfinite(lores)) { // don't allow this to be converted
											// to NaN
				return lores;
			}

			/* Do a taylor series expansion around xpa */
			/* f(x+y) = f(x) + f'(x)*y + f''(x)/2 y^2 */
			double fx1 = xpb / xpa;

			double epsilon = 0.5 * fx1 + 1.0;
			epsilon = epsilon * fx1;

			return epsilon + hiPrec[1] + hiPrec[0];
		}

		/* Value is small |x| < 1e6, do a Taylor series centered on 1.0 */
		double y = x * F_1_3 - F_1_2;
		y = y * x + 1.0;
		y = y * x;

		return y;
	}