TR1 Special Math

[Florian Goth]

Hi Jonathan,
yes my copyright assignment has been completed in the meantime.
Cheers,
Florian.


> Gesendet: Freitag, 08. Mai 2015 um 18:04 Uhr
> Von: "Florian Goth" <CaptainSifff@gmx.de>
> An: "Ed Smith-Rowland" <3dw4rd@verizon.net>
> Cc: libstdc++ <libstdc++@gcc.gnu.org>
> Betreff: Aw: Re: TR1 Special Math
>
> Hi all,
> I think this might be the time for merging the attached, regenerated patch for the Zeta function, as well as my contribution of the PolyLog?
> Cheers,
> Florian.
> 
> 
> See also the original post:
> https://gcc.gnu.org/ml/libstdc++/2014-08/msg00033.html
> 
> index 19def70..43a8b18 100644
> --- a/libstdc++-v3/include/tr1/riemann_zeta.tcc
> +++ b/libstdc++-v3/include/tr1/riemann_zeta.tcc
> @@ -155,8 +155,8 @@ namespace tr1
>  
>        const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
>        //  Max e exponent before overflow.
> -      const _Tp __max_bincoeff = std::numeric_limits<_Tp>::max_exponent10
> -                               * std::log(_Tp(10)) - _Tp(1);
> +      const _Tp __max_bincoeff = std::exp(std::numeric_limits<_Tp>::max_exponent10
> +                               * std::log(_Tp(10)) - _Tp(1));
>  
>        //  This series works until the binomial coefficient blows up
>        //  so use reflection.
> @@ -182,33 +182,25 @@ namespace tr1
>              }
>          }
>  
> -      _Tp __num = _Tp(0.5L);
> +      _Tp __num = _Tp(0.5L*0.5L);
>        const unsigned int __maxit = 10000;
> -      for (unsigned int __i = 0; __i < __maxit; ++__i)
> +      __zeta = _Tp(0.5L);//zeroth order contribution already calculated
> +      for (unsigned int __i = 1; __i < __maxit; ++__i)
>          {
>            bool __punt = false;
> -          _Tp __sgn = _Tp(1);
> -          _Tp __term = _Tp(0);
> -          for (unsigned int __j = 0; __j <= __i; ++__j)
> +          _Tp __term = _Tp(1.0);//again the zeroth order
> +          _Tp __bincoeff = _Tp(1);
> +          for (unsigned int __j = 1; __j <= __i; ++__j)
>              {
> -#if _GLIBCXX_USE_C99_MATH_TR1
> -              _Tp __bincoeff =  std::tr1::lgamma(_Tp(1 + __i))
> -                              - std::tr1::lgamma(_Tp(1 + __j))
> -                              - std::tr1::lgamma(_Tp(1 + __i - __j));
> -#else
> -              _Tp __bincoeff =  __log_gamma(_Tp(1 + __i))
> -                              - __log_gamma(_Tp(1 + __j))
> -                              - __log_gamma(_Tp(1 + __i - __j));
> -#endif
> -              if (__bincoeff > __max_bincoeff)
> -                {
> -                  //  This only gets hit for x << 0.
> -                  __punt = true;
> -                  break;
> -                }
> -              __bincoeff = std::exp(__bincoeff);
> -              __term += __sgn * __bincoeff * std::pow(_Tp(1 + __j), -__s);
> -              __sgn *= _Tp(-1);
> +              _Tp incr = _Tp(__i - __j + 1)/__j;
> +              __bincoeff *= -incr;
> +              if(std::fabs(__bincoeff) > __max_bincoeff )
> +              {
> +                //  This only gets hit for x << 0.
> +                __punt = true;
> +                break;
> +              }
> +              __term += __bincoeff * std::pow(_Tp(1 + __j), -__s);
>              }
>            if (__punt)
>              break;
> 
> 
> > Gesendet: Freitag, 08. Mai 2015 um 16:05 Uhr
> > Von: "Ed Smith-Rowland" <3dw4rd@verizon.net>
> > An: "Jonathan Wakely" <jwakely.gcc@gmail.com>, libstdc++ <libstdc++@gcc.gnu.org>
> > Betreff: Re: TR1 Special Math
> >
> > On 05/07/2015 12:06 PM, Jonathan Wakely wrote:
> > > Hi Ed,
> > >
> > > The C++ committee is considering the
> > > http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2015/n4437.pdf
> > > proposal to make C++17 include the contents of ISO 29124:2010 (the
> > > special math functions from TR1 that went into a separate standard,
> > > not into C++11).
> > >
> > > What is the status of our TR1 implementation? Is it complete? Good
> > > enough quality to move out of the tr1 sub-dir?
> > >
> > > Even if N4437 isn't accepted for C++17 we could move things around to
> > > turn the TR1 code into an iso29124 implementation, do you think that
> > > would make sense?
> > >
> > That would make absolute sense.
> > I actually have a tree where I've done that.
> > All the functions are in there (29124 removed the hypergeometric 
> > functions.  I'd like to keep those as extensions.
> > I have some bugfixes also.
> > 
> > I have a better version of the Carlson elliptic functions (which are 
> > used in the 29124 elliptic functions).
> > 
> > Ed
> > 
> > 
>