This is the mail archive of the
libstdc++@gcc.gnu.org
mailing list for the libstdc++ project.
Aw: Re: TR1 Special Math
- From: "Florian Goth" <CaptainSifff at gmx dot de>
- To: "Ed Smith-Rowland" <3dw4rd at verizon dot net>
- Cc: libstdc++ <libstdc++ at gcc dot gnu dot org>
- Date: Fri, 8 May 2015 18:04:51 +0200
- Subject: Aw: Re: TR1 Special Math
- Authentication-results: sourceware.org; auth=none
- References: <CAH6eHdSK80f08Fr+2gAiX_+QpmfHOgv+-bTtM0DxtVPmRwTykw at mail dot gmail dot com>, <554CC2B4 dot 80401 at verizon dot net>
- Sensitivity: Normal
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
>
>