Aw: Re: TR1 Special Math

[Florian Goth]

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
> 
>