GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-21 12:24:25 +08:00
Code Issues Packages Projects Releases Wiki Activity

Added zeta function.

Browse Source
This commit is contained in:
Till Hoffmann 2016-04-01 13:32:29 +01:00
parent af4ef540bf
commit dd5d390daf
9 changed files with 256 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
Eigen/src/Core/GenericPacketMath.h
View File

@ -76,6 +76,7 @@ struct default_packet_traits
HasTanh = 0,
HasLGamma = 0,
HasDiGamma = 0,
HasZeta = 0,
HasErf = 0,
HasErfc = 0,
HasIGamma = 0,
@ -451,6 +452,10 @@ Packet plgamma(const Packet& a) { using numext::lgamma; return lgamma(a); }
template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet pdigamma(const Packet& a) { using numext::digamma; return digamma(a); }
/** \internal \returns the zeta function of two arguments (coeff-wise) */
template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet pzeta(const Packet& x, const Packet& q) { using numext::zeta; return zeta(x, q); }
/** \internal \returns the erf(\a a) (coeff-wise) */
template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet perf(const Packet& a) { using numext::erf; return erf(a); }

1
Eigen/src/Core/GlobalFunctions.h
View File

@ -51,6 +51,7 @@ namespace Eigen
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(tanh,scalar_tanh_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(lgamma,scalar_lgamma_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(digamma,scalar_digamma_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(zeta,scalar_zeta_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(erf,scalar_erf_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(erfc,scalar_erfc_op)
EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(exp,scalar_exp_op)

189
Eigen/src/Core/SpecialFunctions.h
View File

@ -722,6 +722,189 @@ struct igamma_impl {
#endif // EIGEN_HAS_C99_MATH
/****************************************************************************
* Implementation of Riemann zeta function of two arguments *
****************************************************************************/
template <typename Scalar>
struct zeta_retval {
typedef Scalar type;
};
#ifndef EIGEN_HAS_C99_MATH
template <typename Scalar>
struct zeta_impl {
EIGEN_DEVICE_FUNC
static Scalar run(Scalar x) {
EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
THIS_TYPE_IS_NOT_SUPPORTED);
return Scalar(0);
}
};
#else
template <typename Scalar>
struct zeta_impl {
EIGEN_DEVICE_FUNC
static Scalar run(Scalar x, Scalar q) {
/* zeta.c
*
* Riemann zeta function of two arguments
*
*
*
* SYNOPSIS:
*
* double x, q, y, zeta();
*
* y = zeta( x, q );
*
*
*
* DESCRIPTION:
*
*
*
* inf.
* - -x
* zeta(x,q) = > (k+q)
* -
* k=0
*
* where x > 1 and q is not a negative integer or zero.
* The Euler-Maclaurin summation formula is used to obtain
* the expansion
*
* n
* - -x
* zeta(x,q) = > (k+q)
* -
* k=1
*
* 1-x inf. B x(x+1)...(x+2j)
* (n+q) 1 - 2j
* + --------- - ------- + > --------------------
* x-1 x - x+2j+1
* 2(n+q) j=1 (2j)! (n+q)
*
* where the B2j are Bernoulli numbers. Note that (see zetac.c)
* zeta(x,1) = zetac(x) + 1.
*
*
*
* ACCURACY:
*
*
*
* REFERENCE:
*
* Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals,
* Series, and Products, p. 1073; Academic Press, 1980.
*
*/
int i;
/*double a, b, k, s, t, w;*/
Scalar p, r, a, b, k, s, t, w;
const double A[] = {
12.0,
-720.0,
30240.0,
-1209600.0,
47900160.0,
-1.8924375803183791606e9, /*1.307674368e12/691*/
7.47242496e10,
-2.950130727918164224e12, /*1.067062284288e16/3617*/
1.1646782814350067249e14, /*5.109094217170944e18/43867*/
-4.5979787224074726105e15, /*8.028576626982912e20/174611*/
1.8152105401943546773e17, /*1.5511210043330985984e23/854513*/
-7.1661652561756670113e18 /*1.6938241367317436694528e27/236364091*/
};
const Scalar maxnum = NumTraits<Scalar>::infinity();
const Scalar zero = 0.0, half = 0.5, one = 1.0;
const Scalar machep = igamma_helper<Scalar>::machep();
if( x == one )
return maxnum; //goto retinf;
if( x < one )
{
// domerr:
// mtherr( "zeta", DOMAIN );
return zero;
}
if( q <= zero )
{
if(q == numext::floor(q))
{
// mtherr( "zeta", SING );
// retinf:
return maxnum;
}
p = x;
r = numext::floor(p);
// if( x != floor(x) )
// goto domerr; /* because q^-x not defined */
if (p != r)
return zero;
}
/* Euler-Maclaurin summation formula */
/*
if( x < 25.0 )
*/
{
/* Permit negative q but continue sum until n+q > +9 .
* This case should be handled by a reflection formula.
* If q<0 and x is an integer, there is a relation to
* the polygamma function.
*/
s = numext::pow( q, -x );
a = q;
i = 0;
b = zero;
while( (i < 9) || (a <= Scalar(9.0)) )
{
i += 1;
a += one;
b = numext::pow( a, -x );
s += b;
if( numext::abs(b/s) < machep )
return s; // goto done;
}
w = a;
s += b*w/(x-one);
s -= half * b;
a = one;
k = zero;
for( i=0; i<12; i++ )
{
a *= x + k;
b /= w;
t = a*b/A[i];
s = s + t;
t = numext::abs(t/s);
if( t < machep )
return s; // goto done;
k += one;
a *= x + k;
b /= w;
k += one;
}
// done:
return(s);
}
}
};
#endif // EIGEN_HAS_C99_MATH
} // end namespace internal
namespace numext {
@ -738,6 +921,12 @@ EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(digamma, Scalar)
return EIGEN_MATHFUNC_IMPL(digamma, Scalar)::run(x);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(zeta, Scalar)
zeta(const Scalar& x, const Scalar& q) {
return EIGEN_MATHFUNC_IMPL(zeta, Scalar)::run(x, q);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erf, Scalar)
erf(const Scalar& x) {

14
Eigen/src/Core/arch/CUDA/MathFunctions.h
View File

@ -92,6 +92,20 @@ double2 pdigamma<double2>(const double2& a)
return make_double2(digamma(a.x), digamma(a.y));
}
template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
float4 pzeta<float4>(const float4& a)
{
using numext::zeta;
return make_float4(zeta(a.x), zeta(a.y), zeta(a.z), zeta(a.w));
}
template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
double2 pzeta<double2>(const double2& a)
{
using numext::zeta;
return make_double2(zeta(a.x), zeta(a.y));
}
template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
float4 perf<float4>(const float4& a)
{

1
Eigen/src/Core/arch/CUDA/PacketMath.h
View File

@ -40,6 +40,7 @@ template<> struct packet_traits<float> : default_packet_traits
HasRsqrt = 1,
HasLGamma = 1,
HasDiGamma = 1,
HasZeta = 1,
HasErf = 1,
HasErfc = 1,
HasIgamma = 1,

22
Eigen/src/Core/functors/UnaryFunctors.h
View File

@ -449,6 +449,28 @@ struct functor_traits<scalar_digamma_op<Scalar> >
};
};
/** \internal
* \brief Template functor to compute the Riemann Zeta function of two arguments.
* \sa class CwiseUnaryOp, Cwise::zeta()
*/
template<typename Scalar> struct scalar_zeta_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_zeta_op)
EIGEN_DEVICE_FUNC inline const Scalar operator() (const Scalar& x, const Scalar& q) const {
using numext::zeta; return zeta(x, q);
}
typedef typename packet_traits<Scalar>::type Packet;
EIGEN_DEVICE_FUNC inline Packet packetOp(const Packet& x, const Packet& q) const { return internal::pzeta(x, q); }
};
template<typename Scalar>
struct functor_traits<scalar_zeta_op<Scalar> >
{
enum {
// Guesstimate
Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
PacketAccess = packet_traits<Scalar>::HasZeta
};
};
/** \internal
* \brief Template functor to compute the Gauss error function of a
* scalar

9
Eigen/src/plugins/ArrayCwiseUnaryOps.h
View File

@ -23,6 +23,7 @@ typedef CwiseUnaryOp<internal::scalar_sinh_op<Scalar>, const Derived> SinhReturn
typedef CwiseUnaryOp<internal::scalar_cosh_op<Scalar>, const Derived> CoshReturnType;
typedef CwiseUnaryOp<internal::scalar_lgamma_op<Scalar>, const Derived> LgammaReturnType;
typedef CwiseUnaryOp<internal::scalar_digamma_op<Scalar>, const Derived> DigammaReturnType;
typedef CwiseUnaryOp<internal::scalar_zeta_op<Scalar>, const Derived> ZetaReturnType;
typedef CwiseUnaryOp<internal::scalar_erf_op<Scalar>, const Derived> ErfReturnType;
typedef CwiseUnaryOp<internal::scalar_erfc_op<Scalar>, const Derived> ErfcReturnType;
typedef CwiseUnaryOp<internal::scalar_pow_op<Scalar>, const Derived> PowReturnType;
@ -329,6 +330,14 @@ digamma() const
return DigammaReturnType(derived());
}
/** \returns an expression of the coefficient-wise zeta function.
*/
inline const ZetaReturnType
zeta() const
{
return ZetaReturnType(derived());
}
/** \returns an expression of the coefficient-wise Gauss error
* function of *this.
*

9
test/array.cpp
View File

@ -323,6 +323,15 @@ template<typename ArrayType> void array_real(const ArrayType& m)
VERIFY_IS_EQUAL(numext::digamma(Scalar(-1)),
std::numeric_limits<RealScalar>::infinity());
// Check the zeta function against scipy.special.zeta
VERIFY_IS_APPROX(numext::zeta(Scalar(1.5), Scalar(2)), RealScalar(1.61237534869));
VERIFY_IS_APPROX(numext::zeta(Scalar(4), Scalar(1.5)), RealScalar(0.234848505667));
VERIFY_IS_APPROX(numext::zeta(Scalar(10.5), Scalar(3)), RealScalar(1.03086757337e-5));
VERIFY_IS_APPROX(numext::zeta(Scalar(10000.5), Scalar(1.0001)), RealScalar(0.367879440865));
VERIFY_IS_APPROX(numext::zeta(Scalar(3), Scalar(-2.5)), RealScalar(0.054102025820864097));
VERIFY_IS_EQUAL(numext::zeta(Scalar(1), Scalar(1.2345)), // The second scalar does not matter
std::numeric_limits<RealScalar>::infinity());
{
// Test various propreties of igamma & igammac. These are normalized
// gamma integrals where

6
unsupported/Eigen/CXX11/src/Tensor/TensorBase.h
View File

@ -133,6 +133,12 @@ class TensorBase<Derived, ReadOnlyAccessors>
return unaryExpr(internal::scalar_digamma_op<Scalar>());
}
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE const TensorCwiseUnaryOp<internal::scalar_zeta_op<Scalar>, const Derived>
zeta() const {
return unaryExpr(internal::scalar_zeta_op<Scalar>());
}
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE const TensorCwiseUnaryOp<internal::scalar_erf_op<Scalar>, const Derived>
erf() const {