function y = dlog1p(x)
    % derivate of ln(1+x) wrt to x where x is a vector, and perhaps is small
    y = zeros(size(x));
    ii = abs(x) > eps(1)^.33;
    y(ii==1) = 1./(1+x(ii==1));
    y(ii==0) = 1 - x(ii==0);
  end