Re: Weirdness with numeric_limits in new special functions

[Ed Smith-Rowland <3dw4rd at verizon dot net>]

Paolo Carlini wrote:
> 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.
The public functions in std::tr1: call these with floating point types only.

I think the errtol *is* a mistake because __eps / 8 will degenerate to 
0.  The errtol should wind up being rather larger than __eps.

It might be a good idea to block calling with integral _Tp somehow.

Actually, let me walk through the code and make sure what's going on.

Ed