Weirdness with numeric_limits in new special functions

[Paolo Carlini <pcarlini at suse dot de>]

Hi all, hi Ed,

as part of my ongoing work on razionalizing the headers, I'm looking 
into some code in special functions which I don't understand:

   template<typename _Tp>
   _Tp
   __ellint_rd(const _Tp __x, const _Tp __y, const _Tp __z)
   {
     const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
     const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
     const _Tp __min = std::numeric_limits<_Tp>::min();
     const _Tp __max = std::numeric_limits<_Tp>::max();
     const _Tp __lolim = _Tp(2) / std::pow(__max, _Tp(2) / _Tp(3));
     const _Tp __uplim = std::pow(_Tp(0.1L) * __errtol / __min, _Tp(2) 
/ _Tp(3));
   ...

In this case, and others I suspect, something is very strange: if _Tp is 
a floating point type, for example __eps will be very small and the next 
std::pow for __errtol doesn't appear to make much sense, frankly; worse 
for __min, __max, __lolim, ... On the other hand, if _Tp is integral, 
epsilon() degenerates to zero, all the fractions also become zero, as 
far as I can see, most of the computation is either undefined or trivial.

I'd like to have a clarification, which are the right types for those 
functions (which ones exactly, could you please help auditing?), double 
check that we are not calling with the wrong types, hopefully, help a 
bit adjusting the code...

Thanks,
Paolo.