restore signgam and errno when calling {l,t}gamma

  if (x.dpart() == T(0))

    return v;

  else {

    int errno_saved = errno;

    T h(T(1) / (1ull << (std::numeric_limits<T>::digits / 3)));

    return dual<T>(v,

                   x.dpart() * (tgamma(x.rpart()+h) - tgamma(x.rpart()-h))/(2*h));

    T w((tgamma(x.rpart()+h) - tgamma(x.rpart()-h))/(2*h));

    errno = errno_saved;

    return dual<T>(v, x.dpart() * w);

  }

}

/// Log of absolute value of gamma function.  Approximation of the dual part.

template<class T> dual<T> lgamma(const dual<T> & x) {

  using std::lgamma;

  T h(T(1) / (1ull << (std::numeric_limits<T>::digits / 3)));

  T v = lgamma(x.rpart());

  if (x.dpart() == T(0))

    return v;

  else {

    int signgam_saved = signgam;

    T w = (lgamma(x.rpart()+h) - lgamma(x.rpart()-h))/(2*h);

    int errno_saved = errno;

    T h(T(1) / (1ull << (std::numeric_limits<T>::digits / 3)));

    T w((lgamma(x.rpart()+h) - lgamma(x.rpart()-h))/(2*h));

    signgam = signgam_saved;

    errno = errno_saved;

    return dual<T>(v, x.dpart() * w);

  }

}