* @return log(1 + x)
     */
    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;
    }