diff --git a/math/mathcore/inc/Math/PdfFuncMathCore.h b/math/mathcore/inc/Math/PdfFuncMathCore.h
index 00cf35f96d96fa4d63a9e447dd6b0f98bddeb94c..b2267ef1523a00a587a862f5d7d860c7fc43c2bf 100644
--- a/math/mathcore/inc/Math/PdfFuncMathCore.h
+++ b/math/mathcore/inc/Math/PdfFuncMathCore.h
@@ -249,20 +249,24 @@ namespace Math {
/**
- Probability density function of the Landau distribution.
-
- \f[ p(x) = \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} e^{x s + s \log{s}} ds\f]
-
-
- Where s = (x-x0)/sigma. For detailed description see
+ Probability density function of the Landau distribution:
+ \f[ p(x) = \frac{1}{sigma} \phi (\lambda) \f]
+ with
+ \f[ \phi(\lambda) = \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} e^{\lambda s + s \log{s}} ds\f]
+ where \f$\lambda = (x-x0)/sigma\f$. For a detailed description see
+ K.S. K\"olbig and B. Schorr, A program package for the Landau distribution,
+ <A HREF="http://dx.doi.org/10.1016/0010-4655(84)90085-7">Computer Phys. Comm. 31 (1984) 97-111</A>
+ <A HREF="http://dx.doi.org/10.1016/j.cpc.2008.03.002">[Erratum-ibid. 178 (2008) 972]</A>.
+ The same algorithms as in
<A HREF="http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/g110/top.html">
- CERNLIB</A>. The same algorithms as in CERNLIB (DENLAN) is used
+ CERNLIB</A> (DENLAN) is used
@ingroup PdfFunc
*/
- double landau_pdf(double x, double s = 1, double x0 = 0.);
+ double landau_pdf(double x, double sigma = 1, double x0 = 0);
+
/**
@@ -362,11 +366,6 @@ namespace Math {
-
-
-
-
-
} // namespace Math
} // namespace ROOT
diff --git a/math/mathcore/inc/Math/ProbFuncMathCore.h b/math/mathcore/inc/Math/ProbFuncMathCore.h
index 47a2ce6cc32bea54af98c64681796aea8c74314d..d8b3e582ca0d45b6c23660089bf5127105928354 100644
--- a/math/mathcore/inc/Math/ProbFuncMathCore.h
+++ b/math/mathcore/inc/Math/ProbFuncMathCore.h
@@ -52,6 +52,8 @@ namespace Math {
* These names are currently kept for backward compatibility, but
* their usage is deprecated.
*
+ * These functions are defined in the header file <em>Math/ProbFunc.h<em> or in the global one
+ * including all statistical dunctions <em>Math/DistFunc.h<em>
*
*/
@@ -323,15 +325,20 @@ namespace Math {
Cumulative distribution function of the Landau
distribution (lower tail).
-
+
\f[ D(x) = \int_{-\infty}^{x} p(x) dx \f]
- Where p(x) is the Landau probability density function :
-
- \f[ p(x) = \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} e^{x s + s \log{s}} ds\f]
- Where s = (x-x0)/sigma. For detailed description see
+ where \f$p(x)\f$ is the Landau probability density function :
+ \f[ p(x) = \frac{1}{sigma} \phi (\lambda) \f]
+ with
+ \f[ \phi(\lambda) = \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} e^{\lambda s + s \log{s}} ds\f]
+ with \f$\lambda = (x-x0)/sigma\f$. For a detailed description see
+ K.S. K\"olbig and B. Schorr, A program package for the Landau distribution,
+ <A HREF="http://dx.doi.org/10.1016/0010-4655(84)90085-7">Computer Phys. Comm. 31 (1984) 97-111</A>
+ <A HREF="http://dx.doi.org/10.1016/j.cpc.2008.03.002">[Erratum-ibid. 178 (2008) 972]</A>.
+ The same algorithms as in
<A HREF="http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/g110/top.html">
- CERNLIB</A>. The same algorithms as in CERNLIB (DISLAN) is used.
+ CERNLIB</A> (DISLAN) is used.
@ingroup ProbFunc
@@ -648,6 +655,65 @@ namespace Math {
#endif
+ /** @defgroup TruncFunc Statistical functions from truncated distributions
+
+ @ingroup StatFunc
+
+ Statistical functions for the truncated distributions. Examples of such functions are the
+ first or the second momentum of the truncated distribution.
+ In the case of the Landau, first and second momentum functions are provided for the Landau
+ distribution truncated only on the right side.
+ These functions are defined in the header file <em>Math/ProbFunc.h<em> or in the global one
+ including all statistical dunctions <em>Math/StatFunc.h<em>
+
+ */
+
+ /**
+
+ First moment (mean) of the truncated Landau distribution.
+ \f[ \frac{1}{D (x)} int_{-infty}^{x} t p(t) d t \f]
+ where \f$p(x)\f$ is the Landau distribution:
+ and \f$D(x)\f$ its cumulative distribution function.
+
+ For detailed description see
+ K.S. K\"olbig and B. Schorr, A program package for the Landau distribution,
+ <A HREF="http://dx.doi.org/10.1016/0010-4655(84)90085-7">Computer Phys. Comm. 31 (1984) 97-111</A>
+ <A HREF="http://dx.doi.org/10.1016/j.cpc.2008.03.002">[Erratum-ibid. 178 (2008) 972]</A>.
+ The same algorithms as in
+ <A HREF="http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/g110/top.html">
+ CERNLIB</A> (XM1LAN) is used
+
+ @ingroup TruncFunc
+
+ */
+
+ double landau_xm1(double x, double sigma = 1, double x0 = 0);
+
+
+
+ /**
+
+ Second moment of the truncated Landau distribution.
+ \f[ \frac{1}{D (x)} int_{-infty}^{x} t p(t) d t \f]
+ where \f$p(x)\f$ is the Landau distribution:
+ and \f$D(x)\f$ its cumulative distribution function.
+
+ For detailed description see
+ K.S. K\"olbig and B. Schorr, A program package for the Landau distribution,
+ <A HREF="http://dx.doi.org/10.1016/0010-4655(84)90085-7">Computer Phys. Comm. 31 (1984) 97-111</A>
+ <A HREF="http://dx.doi.org/10.1016/j.cpc.2008.03.002">[Erratum-ibid. 178 (2008) 972]</A>.
+ The same algorithms as in
+ <A HREF="http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/g110/top.html">
+ CERNLIB</A> (XM1LAN) is used
+
+ @ingroup TruncFunc
+
+ */
+
+ double landau_xm2(double x, double sigma = 1, double x0 = 0);
+
+
+
} // namespace Math
} // namespace ROOT
diff --git a/math/mathcore/inc/Math/QuantFuncMathCore.h b/math/mathcore/inc/Math/QuantFuncMathCore.h
index e30d0c916dbc9aee189f9e33d99f14a959cf6c71..d92ae9cabcb53342510f559e45691d7eec01b6ce 100644
--- a/math/mathcore/inc/Math/QuantFuncMathCore.h
+++ b/math/mathcore/inc/Math/QuantFuncMathCore.h
@@ -57,19 +57,22 @@ namespace Math {
*
* \f[ D(x) = \int_{x}^{+\infty} p(x') dx' \f]
*
- * The implementation used is that of
- * <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Distributions.html">GSL</A>.
+ * These functions are defined in the header file <em>Math/ProbFunc.h<em> or in the global one
+ * including all statistical dunctions <em>Math/DistFunc.h<em>
+ *
*
* <strong>NOTE:</strong> In the old releases (< 5.14) the <em>_quantile</em> functions were called
* <em>_quant_inv</em> and the <em>_quantile_c</em> functions were called
* <em>_prob_inv</em>.
* These names are currently kept for backward compatibility, but
* their usage is deprecated.
+ *
*/
/** @name Quantile Functions from MathCore
- * The implementation used is that of
- * <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Distributions.html">GSL</A>.
+ * The implementation is provided in MathCore and for the majority of the function comes from
+ * <A HREF="http://www.netlib.org/cephes">Cephes</A>.
+
*/
//@{
@@ -521,6 +524,41 @@ namespace Math {
+
+ /**
+
+ Inverse (\f$D^{-1}(z)\f$) of the cumulative distribution
+ function of the lower tail of the Landau distribution
+ (#landau_cdf).
+
+ For detailed description see
+ K.S. K\"olbig and B. Schorr, A program package for the Landau distribution,
+ <A HREF="http://dx.doi.org/10.1016/0010-4655(84)90085-7">Computer Phys. Comm. 31 (1984) 97-111</A>
+ <A HREF="http://dx.doi.org/10.1016/j.cpc.2008.03.002">[Erratum-ibid. 178 (2008) 972]</A>.
+ The same algorithms as in
+ <A HREF="http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/g110/top.html">
+ CERNLIB</A> (RANLAN) is used.
+
+
+ @ingroup QuantFunc
+
+ */
+
+ double landau_quantile(double z, double sigma = 1);
+
+
+ /**
+
+ Inverse (\f$D^{-1}(z)\f$) of the cumulative distribution
+ function of the upper tail of the landau distribution
+ (#landau_cdf_c).
+ Implemented using #landau_quantile
+
+ */
+
+ double landau_quantile_c(double z, double sigma = 1);
+
+
#ifdef HAVE_OLD_STAT_FUNC
//@}
diff --git a/math/mathcore/inc/Math/SpecFuncMathCore.h b/math/mathcore/inc/Math/SpecFuncMathCore.h
index ec89b4a61f9990a1dfd97795f53f31c7f7003c2a..333ab86ed99e2fbcc58428214b3564620184a1f7 100644
--- a/math/mathcore/inc/Math/SpecFuncMathCore.h
+++ b/math/mathcore/inc/Math/SpecFuncMathCore.h
@@ -183,6 +183,55 @@ namespace Math {
double inc_beta( double x, double a, double b);
+
+
+
+ /**
+
+ Calculates the sine integral.
+
+ \f[ Si(x) = - \int_{0}^{x} \frac{\sin t}{t} dt \f]
+
+ For detailed description see
+ <A HREF="http://mathworld.wolfram.com/SineIntegral.html">
+ Mathworld</A>. The implementation used is that of
+ <A HREF="http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/c336/top.html">
+ CERNLIB</A>,
+ based on Y.L. Luke, The special functions and their approximations, v.II, (Academic Press, New York l969) 325-326.
+
+
+ @ingroup SpecFunc
+
+ */
+
+ double sinint(double x);
+
+
+
+
+ /**
+
+ Calculates the real part of the cosine integral \Re(Ci).
+
+ For x<0, the imaginary part is \pi i and has to be added by the user,
+ for x>0 the imaginary part of Ci(x) is 0.
+
+ \f[ Ci(x) = - \int_{x}^{\infty} \frac{\cos t}{t} dt = \gamma + \ln x + \int_{0}^{x} \frac{\cos t - 1}{t} dt\f]
+
+ For detailed description see
+ <A HREF="http://mathworld.wolfram.com/CosineIntegral.html">
+ Mathworld</A>. The implementation used is that of
+ <A HREF="http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/c336/top.html">
+ CERNLIB</A>,
+ based on Y.L. Luke, The special functions and their approximations, v.II, (Academic Press, New York l969) 325-326.
+
+
+ @ingroup SpecFunc
+
+ */
+
+ double cosint(double x);
+
diff --git a/math/mathcore/src/PdfFuncMathCore.cxx b/math/mathcore/src/PdfFuncMathCore.cxx
index 96c8ca3ad57e3cf168fdb0bb2ed49058f22c2b2a..cd594743202ad1101ec0635e12d5d7727b65d615 100644
--- a/math/mathcore/src/PdfFuncMathCore.cxx
+++ b/math/mathcore/src/PdfFuncMathCore.cxx
@@ -188,8 +188,9 @@ namespace Math {
denlan = u*u*(1+(a2[0]+a2[1]*u)*u);
}
return denlan/sigma;
-}
-
+
+ }
+
@@ -245,7 +246,8 @@ namespace Math {
}
-
+
+
} // namespace Math
} // namespace ROOT
diff --git a/math/mathcore/src/ProbFuncMathCore.cxx b/math/mathcore/src/ProbFuncMathCore.cxx
index b6bfbfeba784157b707d9cbf202b764190ac48b9..5de6629a7f49bac09bd4db85ef7a34e57dd0b732 100644
--- a/math/mathcore/src/ProbFuncMathCore.cxx
+++ b/math/mathcore/src/ProbFuncMathCore.cxx
@@ -273,26 +273,26 @@ namespace Math {
//distribution", Computer Phys.Comm., 31(1984), 97-111
static double p1[5] = {0.2514091491e+0,-0.6250580444e-1, 0.1458381230e-1, -0.2108817737e-2, 0.7411247290e-3};
- static double q1[5] = {1.0 ,-0.5571175625e-2, 0.6225310236e-1, -0.3137378427e-2, 0.1931496439e-2};
+ static double q1[5] = {1.0 ,-0.5571175625e-2, 0.6225310236e-1, -0.3137378427e-2, 0.1931496439e-2};
static double p2[4] = {0.2868328584e+0, 0.3564363231e+0, 0.1523518695e+0, 0.2251304883e-1};
- static double q2[4] = {1.0 , 0.6191136137e+0, 0.1720721448e+0, 0.2278594771e-1};
+ static double q2[4] = {1.0 , 0.6191136137e+0, 0.1720721448e+0, 0.2278594771e-1};
static double p3[4] = {0.2868329066e+0, 0.3003828436e+0, 0.9950951941e-1, 0.8733827185e-2};
- static double q3[4] = {1.0 , 0.4237190502e+0, 0.1095631512e+0, 0.8693851567e-2};
+ static double q3[4] = {1.0 , 0.4237190502e+0, 0.1095631512e+0, 0.8693851567e-2};
static double p4[4] = {0.1000351630e+1, 0.4503592498e+1, 0.1085883880e+2, 0.7536052269e+1};
- static double q4[4] = {1.0 , 0.5539969678e+1, 0.1933581111e+2, 0.2721321508e+2};
+ static double q4[4] = {1.0 , 0.5539969678e+1, 0.1933581111e+2, 0.2721321508e+2};
static double p5[4] = {0.1000006517e+1, 0.4909414111e+2, 0.8505544753e+2, 0.1532153455e+3};
- static double q5[4] = {1.0 , 0.5009928881e+2, 0.1399819104e+3, 0.4200002909e+3};
+ static double q5[4] = {1.0 , 0.5009928881e+2, 0.1399819104e+3, 0.4200002909e+3};
static double p6[4] = {0.1000000983e+1, 0.1329868456e+3, 0.9162149244e+3, -0.9605054274e+3};
- static double q6[4] = {1.0 , 0.1339887843e+3, 0.1055990413e+4, 0.5532224619e+3};
+ static double q6[4] = {1.0 , 0.1339887843e+3, 0.1055990413e+4, 0.5532224619e+3};
static double a1[4] = {0, -0.4583333333e+0, 0.6675347222e+0,-0.1641741416e+1};
- static double a2[4] = {0, 1.0 ,-0.4227843351e+0,-0.2043403138e+1};
+ static double a2[4] = {0, 1.0 ,-0.4227843351e+0,-0.2043403138e+1};
double v = (x - x0)/sigma;
double u;
@@ -332,6 +332,170 @@ namespace Math {
}
+
+ double landau_xm1(double x, double sigma, double x0) {
+ // implementation of first momentum of Landau distribution
+ // translated from Cernlib (XM1LAN function) by Benno List
+
+ static double p1[5] = {
+ -0.8949374280E+0, 0.4631783434E+0,-0.4053332915E-1,
+ 0.1580075560E-1,-0.3423874194E-2};
+ static double q1[5] = {
+ 1.0 , 0.1002930749E+0, 0.3575271633E-1,
+ -0.1915882099E-2, 0.4811072364E-4};
+ static double p2[5] = {
+ -0.8933384046E+0, 0.1161296496E+0, 0.1200082940E+0,
+ 0.2185699725E-1, 0.2128892058E-2};
+ static double q2[5] = {
+ 1.0 , 0.4935531886E+0, 0.1066347067E+0,
+ 0.1250161833E-1, 0.5494243254E-3};
+ static double p3[5] = {
+ -0.8933322067E+0, 0.2339544896E+0, 0.8257653222E-1,
+ 0.1411226998E-1, 0.2892240953E-3};
+ static double q3[5] = {
+ 1.0 , 0.3616538408E+0, 0.6628026743E-1,
+ 0.4839298984E-2, 0.5248310361E-4};
+ static double p4[4] = {
+ 0.9358419425E+0, 0.6716831438E+2,-0.6765069077E+3,
+ 0.9026661865E+3};
+ static double q4[4] = {
+ 1.0 , 0.7752562854E+2,-0.5637811998E+3,
+ -0.5513156752E+3};
+ static double p5[4] = {
+ 0.9489335583E+0, 0.5561246706E+3, 0.3208274617E+5,
+ -0.4889926524E+5};
+ static double q5[4] = {
+ 1.0 , 0.6028275940E+3, 0.3716962017E+5,
+ 0.3686272898E+5};
+ static double a0[6] = {
+ -0.4227843351E+0,-0.1544313298E+0, 0.4227843351E+0,
+ 0.3276496874E+1, 0.2043403138E+1,-0.8681296500E+1};
+ static double a1[4] = {0,
+ -0.4583333333E+0, 0.6675347222E+0,-0.1641741416E+1};
+ static double a2[5] = {0,
+ -0.1958333333E+1, 0.5563368056E+1,-0.2111352961E+2,
+ 0.1006946266E+3};
+
+ double v = (x-x0)/sigma;
+ double xm1lan;
+ if (v < -4.5) {
+ double u = std::exp(v+1);
+ xm1lan = v-u*(1+(a2[1]+(a2[2]+(a2[3]+a2[4]*u)*u)*u)*u)/
+ (1+(a1[1]+(a1[2]+a1[3]*u)*u)*u);
+ } else if(v < -2) {
+ xm1lan = (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
+ (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
+ } else if(v < 2) {
+ xm1lan = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*v)*v)*v)*v)/
+ (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*v)*v)*v)*v);
+ } else if(v < 10) {
+ xm1lan = (p3[0]+(p3[1]+(p3[2]+(p3[3]+p3[4]*v)*v)*v)*v)/
+ (q3[0]+(q3[1]+(q3[2]+(q3[3]+q3[4]*v)*v)*v)*v);
+ } else if(v < 40) {
+ double u = 1/v;
+ xm1lan = std::log(v)*(p4[0]+(p4[1]+(p4[2]+p4[3]*u)*u)*u)/
+ (q4[0]+(q4[1]+(q4[2]+q4[3]*u)*u)*u);
+ } else if(v < 200) {
+ double u = 1/v;
+ xm1lan= std::log(v)*(p5[0]+(p5[1]+(p5[2]+p5[3]*u)*u)*u)/
+ (q5[0]+(q5[1]+(q5[2]+q5[3]*u)*u)*u);
+ } else {
+ double u = v-v*std::log(v)/(v+1);
+ v = 1/(u-u*(u+ std::log(u)-v)/(u+1));
+ u = -std::log(v);
+ xm1lan = (u+a0[0]+(-u+a0[1]+(a0[2]*u+a0[3]+(a0[4]*u+a0[5])*v)*v)*v)/
+ (1-(1-(a0[2]+a0[4]*v)*v)*v);
+ }
+ return xm1lan*sigma + x0;
+
+ }
+
+
+
+ double landau_xm2(double x, double sigma, double x0) {
+ // implementation of second momentum of Landau distribution
+ // translated from Cernlib (XM2LAN function) by Benno List
+
+ static double p1[5] = {
+ 0.1169837582E+1,-0.4834874539E+0, 0.4383774644E+0,
+ 0.3287175228E-2, 0.1879129206E-1};
+ static double q1[5] = {
+ 1.0 , 0.1795154326E+0, 0.4612795899E-1,
+ 0.2183459337E-2, 0.7226623623E-4};
+ static double p2[5] = {
+ 0.1157939823E+1,-0.3842809495E+0, 0.3317532899E+0,
+ 0.3547606781E-1, 0.6725645279E-2};
+ static double q2[5] = {
+ 1.0 , 0.2916824021E+0, 0.5259853480E-1,
+ 0.3840011061E-2, 0.9950324173E-4};
+ static double p3[4] = {
+ 0.1178191282E+1, 0.1011623342E+2,-0.1285585291E+2,
+ 0.3641361437E+2};
+ static double q3[4] = {
+ 1.0 , 0.8614160194E+1, 0.3118929630E+2,
+ 0.1514351300E+0};
+ static double p4[5] = {
+ 0.1030763698E+1, 0.1216758660E+3, 0.1637431386E+4,
+ -0.2171466507E+4, 0.7010168358E+4};
+ static double q4[5] = {
+ 1.0 , 0.1022487911E+3, 0.1377646350E+4,
+ 0.3699184961E+4, 0.4251315610E+4};
+ static double p5[4] = {
+ 0.1010084827E+1, 0.3944224824E+3, 0.1773025353E+5,
+ -0.7075963938E+5};
+ static double q5[4] = {
+ 1.0 , 0.3605950254E+3, 0.1392784158E+5,
+ -0.1881680027E+5};
+ static double a0[7] = {
+ -0.2043403138E+1,-0.8455686702E+0,-0.3088626596E+0,
+ 0.5821346754E+1, 0.4227843351E+0, 0.6552993748E+1,
+ -0.1076714945E+2};
+ static double a1[4] = {0,
+ -0.4583333333E+0, 0.6675347222E+0,-0.1641741416E+1};
+ static double a2[4] = {
+ -0.1958333333E+1, 0.5563368056E+1,-0.2111352961E+2,
+ 0.1006946266E+3};
+ static double a3[4] = {
+ -1.0 , 0.4458333333E+1,-0.2116753472E+2,
+ 0.1163674359E+3};
+
+ double v = (x-x0)/sigma;
+ double xm2lan;
+ if (v < -4.5) {
+ double u = std::exp(v+1);
+ xm2lan = v*v-2*u*u*
+ (v/u+a2[0]*v+a3[0]+(a2[1]*v+a3[1]+(a2[2]*v+a3[2]+
+ (a2[3]*v+a3[3])*u)*u)*u)/
+ (1+(a1[1]+(a1[2]+a1[3]*u)*u)*u);
+ } else if (v < -2) {
+ xm2lan = (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
+ (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
+ } else if (v < 2) {
+ xm2lan = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*v)*v)*v)*v)/
+ (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*v)*v)*v)*v);
+ } else if (v < 5) {
+ double u = 1/v;
+ xm2lan = v*(p3[0]+(p3[1]+(p3[2]+p3[3]*u)*u)*u)/
+ (q3[0]+(q3[1]+(q3[2]+q3[3]*u)*u)*u);
+ } else if (v < 50) {
+ double u = 1/v;
+ xm2lan = v*(p4[0]+(p4[1]+(p4[2]+(p4[3]+p4[4]*u)*u)*u)*u)/
+ (q4[0]+(q4[1]+(q4[2]+(q4[3]+q4[4]*u)*u)*u)*u);
+ } else if (v < 200) {
+ double u = 1/v;
+ xm2lan = v*(p5[0]+(p5[1]+(p5[2]+p5[3]*u)*u)*u)/
+ (q5[0]+(q5[1]+(q5[2]+q5[3]*u)*u)*u);
+ } else {
+ double u = v-v*std::log(v)/(v+1);
+ v = 1/(u-u*(u+log(u)-v)/(u+1));
+ u = -std::log(v);
+ xm2lan = (1/v+u*u+a0[0]+a0[1]*u+(-u*u+a0[2]*u+a0[3]+
+ (a0[4]*u*u+a0[5]*u+a0[6])*v)*v)/(1-(1-a0[4]*v)*v);
+ }
+ if (x0 == 0) return xm2lan*sigma*sigma;
+ double xm1lan = ROOT::Math::landau_xm1(x, sigma, x0);
+ return xm2lan*sigma*sigma + (2*xm1lan-x0)*x0;
+ }
} // namespace Math
} // namespace ROOT
diff --git a/math/mathcore/src/QuantFuncMathCore.cxx b/math/mathcore/src/QuantFuncMathCore.cxx
index 62454c1e24368c8f05bfea286423b8e3c361ab2c..76c869d8a8327083052d8a2a5a42bdbeaaa25971 100644
--- a/math/mathcore/src/QuantFuncMathCore.cxx
+++ b/math/mathcore/src/QuantFuncMathCore.cxx
@@ -188,6 +188,216 @@ namespace Math {
}
+ double landau_quantile(double z, double sigma) {
+ // LANDAU quantile : algorithm from CERNLIB G110 ranlan
+ // with scale parameter sigma
+ // Converted by Rene Brun from CERNLIB routine ranlan(G110),
+ // Moved and adapted to QuantFuncMathCore by B. List 29.4.2010
+
+ static double f[982] = {
+ 0 , 0 , 0 ,0 ,0 ,-2.244733,
+ -2.204365,-2.168163,-2.135219,-2.104898,-2.076740,-2.050397,
+ -2.025605,-2.002150,-1.979866,-1.958612,-1.938275,-1.918760,
+ -1.899984,-1.881879,-1.864385,-1.847451,-1.831030,-1.815083,
+ -1.799574,-1.784473,-1.769751,-1.755383,-1.741346,-1.727620,
+ -1.714187,-1.701029,-1.688130,-1.675477,-1.663057,-1.650858,
+ -1.638868,-1.627078,-1.615477,-1.604058,-1.592811,-1.581729,
+ -1.570806,-1.560034,-1.549407,-1.538919,-1.528565,-1.518339,
+ -1.508237,-1.498254,-1.488386,-1.478628,-1.468976,-1.459428,
+ -1.449979,-1.440626,-1.431365,-1.422195,-1.413111,-1.404112,
+ -1.395194,-1.386356,-1.377594,-1.368906,-1.360291,-1.351746,
+ -1.343269,-1.334859,-1.326512,-1.318229,-1.310006,-1.301843,
+ -1.293737,-1.285688,-1.277693,-1.269752,-1.261863,-1.254024,
+ -1.246235,-1.238494,-1.230800,-1.223153,-1.215550,-1.207990,
+ -1.200474,-1.192999,-1.185566,-1.178172,-1.170817,-1.163500,
+ -1.156220,-1.148977,-1.141770,-1.134598,-1.127459,-1.120354,
+ -1.113282,-1.106242,-1.099233,-1.092255,
+ -1.085306,-1.078388,-1.071498,-1.064636,-1.057802,-1.050996,
+ -1.044215,-1.037461,-1.030733,-1.024029,-1.017350,-1.010695,
+ -1.004064, -.997456, -.990871, -.984308, -.977767, -.971247,
+ -.964749, -.958271, -.951813, -.945375, -.938957, -.932558,
+ -.926178, -.919816, -.913472, -.907146, -.900838, -.894547,
+ -.888272, -.882014, -.875773, -.869547, -.863337, -.857142,
+ -.850963, -.844798, -.838648, -.832512, -.826390, -.820282,
+ -.814187, -.808106, -.802038, -.795982, -.789940, -.783909,
+ -.777891, -.771884, -.765889, -.759906, -.753934, -.747973,
+ -.742023, -.736084, -.730155, -.724237, -.718328, -.712429,
+ -.706541, -.700661, -.694791, -.688931, -.683079, -.677236,
+ -.671402, -.665576, -.659759, -.653950, -.648149, -.642356,
+ -.636570, -.630793, -.625022, -.619259, -.613503, -.607754,
+ -.602012, -.596276, -.590548, -.584825, -.579109, -.573399,
+ -.567695, -.561997, -.556305, -.550618, -.544937, -.539262,
+ -.533592, -.527926, -.522266, -.516611, -.510961, -.505315,
+ -.499674, -.494037, -.488405, -.482777,
+ -.477153, -.471533, -.465917, -.460305, -.454697, -.449092,
+ -.443491, -.437893, -.432299, -.426707, -.421119, -.415534,
+ -.409951, -.404372, -.398795, -.393221, -.387649, -.382080,
+ -.376513, -.370949, -.365387, -.359826, -.354268, -.348712,
+ -.343157, -.337604, -.332053, -.326503, -.320955, -.315408,
+ -.309863, -.304318, -.298775, -.293233, -.287692, -.282152,
+ -.276613, -.271074, -.265536, -.259999, -.254462, -.248926,
+ -.243389, -.237854, -.232318, -.226783, -.221247, -.215712,
+ -.210176, -.204641, -.199105, -.193568, -.188032, -.182495,
+ -.176957, -.171419, -.165880, -.160341, -.154800, -.149259,
+ -.143717, -.138173, -.132629, -.127083, -.121537, -.115989,
+ -.110439, -.104889, -.099336, -.093782, -.088227, -.082670,
+ -.077111, -.071550, -.065987, -.060423, -.054856, -.049288,
+ -.043717, -.038144, -.032569, -.026991, -.021411, -.015828,
+ -.010243, -.004656, .000934, .006527, .012123, .017722,
+ .023323, .028928, .034535, .040146, .045759, .051376,
+ .056997, .062620, .068247, .073877,
+ .079511, .085149, .090790, .096435, .102083, .107736,
+ .113392, .119052, .124716, .130385, .136057, .141734,
+ .147414, .153100, .158789, .164483, .170181, .175884,
+ .181592, .187304, .193021, .198743, .204469, .210201,
+ .215937, .221678, .227425, .233177, .238933, .244696,
+ .250463, .256236, .262014, .267798, .273587, .279382,
+ .285183, .290989, .296801, .302619, .308443, .314273,
+ .320109, .325951, .331799, .337654, .343515, .349382,
+ .355255, .361135, .367022, .372915, .378815, .384721,
+ .390634, .396554, .402481, .408415, .414356, .420304,
+ .426260, .432222, .438192, .444169, .450153, .456145,
+ .462144, .468151, .474166, .480188, .486218, .492256,
+ .498302, .504356, .510418, .516488, .522566, .528653,
+ .534747, .540850, .546962, .553082, .559210, .565347,
+ .571493, .577648, .583811, .589983, .596164, .602355,
+ .608554, .614762, .620980, .627207, .633444, .639689,
+ .645945, .652210, .658484, .664768,
+ .671062, .677366, .683680, .690004, .696338, .702682,
+ .709036, .715400, .721775, .728160, .734556, .740963,
+ .747379, .753807, .760246, .766695, .773155, .779627,
+ .786109, .792603, .799107, .805624, .812151, .818690,
+ .825241, .831803, .838377, .844962, .851560, .858170,
+ .864791, .871425, .878071, .884729, .891399, .898082,
+ .904778, .911486, .918206, .924940, .931686, .938446,
+ .945218, .952003, .958802, .965614, .972439, .979278,
+ .986130, .992996, .999875, 1.006769, 1.013676, 1.020597,
+ 1.027533, 1.034482, 1.041446, 1.048424, 1.055417, 1.062424,
+ 1.069446, 1.076482, 1.083534, 1.090600, 1.097681, 1.104778,
+ 1.111889, 1.119016, 1.126159, 1.133316, 1.140490, 1.147679,
+ 1.154884, 1.162105, 1.169342, 1.176595, 1.183864, 1.191149,
+ 1.198451, 1.205770, 1.213105, 1.220457, 1.227826, 1.235211,
+ 1.242614, 1.250034, 1.257471, 1.264926, 1.272398, 1.279888,
+ 1.287395, 1.294921, 1.302464, 1.310026, 1.317605, 1.325203,
+ 1.332819, 1.340454, 1.348108, 1.355780,
+ 1.363472, 1.371182, 1.378912, 1.386660, 1.394429, 1.402216,
+ 1.410024, 1.417851, 1.425698, 1.433565, 1.441453, 1.449360,
+ 1.457288, 1.465237, 1.473206, 1.481196, 1.489208, 1.497240,
+ 1.505293, 1.513368, 1.521465, 1.529583, 1.537723, 1.545885,
+ 1.554068, 1.562275, 1.570503, 1.578754, 1.587028, 1.595325,
+ 1.603644, 1.611987, 1.620353, 1.628743, 1.637156, 1.645593,
+ 1.654053, 1.662538, 1.671047, 1.679581, 1.688139, 1.696721,
+ 1.705329, 1.713961, 1.722619, 1.731303, 1.740011, 1.748746,
+ 1.757506, 1.766293, 1.775106, 1.783945, 1.792810, 1.801703,
+ 1.810623, 1.819569, 1.828543, 1.837545, 1.846574, 1.855631,
+ 1.864717, 1.873830, 1.882972, 1.892143, 1.901343, 1.910572,
+ 1.919830, 1.929117, 1.938434, 1.947781, 1.957158, 1.966566,
+ 1.976004, 1.985473, 1.994972, 2.004503, 2.014065, 2.023659,
+ 2.033285, 2.042943, 2.052633, 2.062355, 2.072110, 2.081899,
+ 2.091720, 2.101575, 2.111464, 2.121386, 2.131343, 2.141334,
+ 2.151360, 2.161421, 2.171517, 2.181648, 2.191815, 2.202018,
+ 2.212257, 2.222533, 2.232845, 2.243195,
+ 2.253582, 2.264006, 2.274468, 2.284968, 2.295507, 2.306084,
+ 2.316701, 2.327356, 2.338051, 2.348786, 2.359562, 2.370377,
+ 2.381234, 2.392131, 2.403070, 2.414051, 2.425073, 2.436138,
+ 2.447246, 2.458397, 2.469591, 2.480828, 2.492110, 2.503436,
+ 2.514807, 2.526222, 2.537684, 2.549190, 2.560743, 2.572343,
+ 2.583989, 2.595682, 2.607423, 2.619212, 2.631050, 2.642936,
+ 2.654871, 2.666855, 2.678890, 2.690975, 2.703110, 2.715297,
+ 2.727535, 2.739825, 2.752168, 2.764563, 2.777012, 2.789514,
+ 2.802070, 2.814681, 2.827347, 2.840069, 2.852846, 2.865680,
+ 2.878570, 2.891518, 2.904524, 2.917588, 2.930712, 2.943894,
+ 2.957136, 2.970439, 2.983802, 2.997227, 3.010714, 3.024263,
+ 3.037875, 3.051551, 3.065290, 3.079095, 3.092965, 3.106900,
+ 3.120902, 3.134971, 3.149107, 3.163312, 3.177585, 3.191928,
+ 3.206340, 3.220824, 3.235378, 3.250005, 3.264704, 3.279477,
+ 3.294323, 3.309244, 3.324240, 3.339312, 3.354461, 3.369687,
+ 3.384992, 3.400375, 3.415838, 3.431381, 3.447005, 3.462711,
+ 3.478500, 3.494372, 3.510328, 3.526370,
+ 3.542497, 3.558711, 3.575012, 3.591402, 3.607881, 3.624450,
+ 3.641111, 3.657863, 3.674708, 3.691646, 3.708680, 3.725809,
+ 3.743034, 3.760357, 3.777779, 3.795300, 3.812921, 3.830645,
+ 3.848470, 3.866400, 3.884434, 3.902574, 3.920821, 3.939176,
+ 3.957640, 3.976215, 3.994901, 4.013699, 4.032612, 4.051639,
+ 4.070783, 4.090045, 4.109425, 4.128925, 4.148547, 4.168292,
+ 4.188160, 4.208154, 4.228275, 4.248524, 4.268903, 4.289413,
+ 4.310056, 4.330832, 4.351745, 4.372794, 4.393982, 4.415310,
+ 4.436781, 4.458395, 4.480154, 4.502060, 4.524114, 4.546319,
+ 4.568676, 4.591187, 4.613854, 4.636678, 4.659662, 4.682807,
+ 4.706116, 4.729590, 4.753231, 4.777041, 4.801024, 4.825179,
+ 4.849511, 4.874020, 4.898710, 4.923582, 4.948639, 4.973883,
+ 4.999316, 5.024942, 5.050761, 5.076778, 5.102993, 5.129411,
+ 5.156034, 5.182864, 5.209903, 5.237156, 5.264625, 5.292312,
+ 5.320220, 5.348354, 5.376714, 5.405306, 5.434131, 5.463193,
+ 5.492496, 5.522042, 5.551836, 5.581880, 5.612178, 5.642734,
+ 5.673552, 5.704634, 5.735986, 5.767610,
+ 5.799512, 5.831694, 5.864161, 5.896918, 5.929968, 5.963316,
+ 5.996967, 6.030925, 6.065194, 6.099780, 6.134687, 6.169921,
+ 6.205486, 6.241387, 6.277630, 6.314220, 6.351163, 6.388465,
+ 6.426130, 6.464166, 6.502578, 6.541371, 6.580553, 6.620130,
+ 6.660109, 6.700495, 6.741297, 6.782520, 6.824173, 6.866262,
+ 6.908795, 6.951780, 6.995225, 7.039137, 7.083525, 7.128398,
+ 7.173764, 7.219632, 7.266011, 7.312910, 7.360339, 7.408308,
+ 7.456827, 7.505905, 7.555554, 7.605785, 7.656608, 7.708035,
+ 7.760077, 7.812747, 7.866057, 7.920019, 7.974647, 8.029953,
+ 8.085952, 8.142657, 8.200083, 8.258245, 8.317158, 8.376837,
+ 8.437300, 8.498562, 8.560641, 8.623554, 8.687319, 8.751955,
+ 8.817481, 8.883916, 8.951282, 9.019600, 9.088889, 9.159174,
+ 9.230477, 9.302822, 9.376233, 9.450735, 9.526355, 9.603118,
+ 9.681054, 9.760191, 9.840558, 9.922186,10.005107,10.089353,
+ 10.174959,10.261958,10.350389,10.440287,10.531693,10.624646,
+ 10.719188,10.815362,10.913214,11.012789,11.114137,11.217307,
+ 11.322352,11.429325,11.538283,11.649285,
+ 11.762390,11.877664,11.995170,12.114979,12.237161,12.361791,
+ 12.488946,12.618708,12.751161,12.886394,13.024498,13.165570,
+ 13.309711,13.457026,13.607625,13.761625,13.919145,14.080314,
+ 14.245263,14.414134,14.587072,14.764233,14.945778,15.131877,
+ 15.322712,15.518470,15.719353,15.925570,16.137345,16.354912,
+ 16.578520,16.808433,17.044929,17.288305,17.538873,17.796967,
+ 18.062943,18.337176,18.620068,18.912049,19.213574,19.525133,
+ 19.847249,20.180480,20.525429,20.882738,21.253102,21.637266,
+ 22.036036,22.450278,22.880933,23.329017,23.795634,24.281981,
+ 24.789364,25.319207,25.873062,26.452634,27.059789,27.696581,
+ 28.365274,29.068370,29.808638,30.589157,31.413354,32.285060,
+ 33.208568,34.188705,35.230920,36.341388,37.527131,38.796172,
+ 40.157721,41.622399,43.202525,44.912465,46.769077,48.792279,
+ 51.005773,53.437996,56.123356,59.103894 };
+
+ if (sigma <= 0) return 0;
+ if (z <= 0) return -std::numeric_limits<double>::infinity();
+ if (z >= 1) return std::numeric_limits<double>::infinity();
+
+ double ranlan, u, v;
+ u = 1000*z;
+ int i = int(u);
+ u -= i;
+ if (i >= 70 && i < 800) {
+ ranlan = f[i-1] + u*(f[i] - f[i-1]);
+ } else if (i >= 7 && i <= 980) {
+ ranlan = f[i-1] + u*(f[i]-f[i-1]-0.25*(1-u)*(f[i+1]-f[i]-f[i-1]+f[i-2]));
+ } else if (i < 7) {
+ v = std::log(z);
+ u = 1/v;
+ ranlan = ((0.99858950+(3.45213058E1+1.70854528E1*u)*u)/
+ (1 +(3.41760202E1+4.01244582 *u)*u))*
+ (-std::log(-0.91893853-v)-1);
+ } else {
+ u = 1-z;
+ v = u*u;
+ if (z <= 0.999) {
+ ranlan = (1.00060006+2.63991156E2*u+4.37320068E3*v)/
+ ((1 +2.57368075E2*u+3.41448018E3*v)*u);
+ } else {
+ ranlan = (1.00001538+6.07514119E3*u+7.34266409E5*v)/
+ ((1 +6.06511919E3*u+6.94021044E5*v)*u);
+ }
+ }
+ return sigma*ranlan;
+ }
+
+ double landau_quantile_c(double z, double sigma) {
+ return landau_quantile(1.-z,sigma);
+ }
diff --git a/math/mathcore/src/SpecFuncMathCore.cxx b/math/mathcore/src/SpecFuncMathCore.cxx
index d00cd55dc78d5ca60928fcd888f7683267c61577..24a11e8d26fd2397226b3500c6b42c4a84d366d1 100644
--- a/math/mathcore/src/SpecFuncMathCore.cxx
+++ b/math/mathcore/src/SpecFuncMathCore.cxx
@@ -17,6 +17,11 @@
#include <cmath>
+#include <limits>
+
+#ifndef PI
+#define PI 3.14159265358979323846264338328 /* pi */
+#endif
// use cephes for functions which are also in C99
#define USE_CEPHES
@@ -111,6 +116,188 @@ double inc_beta( double x, double a, double b) {
return ROOT::Math::Cephes::incbet(a,b,x);
}
+// Sine integral
+// Translated from CERNLIB SININT (C336) by B. List 29.4.2010
+
+double sinint(double x) {
+
+ static const double z1 = 1, r8 = z1/8;
+
+ static const double pih = PI/2;
+
+ static const double s[16] = {
+ +1.95222097595307108, -0.68840423212571544,
+ +0.45518551322558484, -0.18045712368387785,
+ +0.04104221337585924, -0.00595861695558885,
+ +0.00060014274141443, -0.00004447083291075,
+ +0.00000253007823075, -0.00000011413075930,
+ +0.00000000418578394, -0.00000000012734706,
+ +0.00000000000326736, -0.00000000000007168,
+ +0.00000000000000136, -0.00000000000000002};
+
+ static const double p[29] = {
+ +0.96074783975203596, -0.03711389621239806,
+ +0.00194143988899190, -0.00017165988425147,
+ +0.00002112637753231, -0.00000327163256712,
+ +0.00000060069211615, -0.00000012586794403,
+ +0.00000002932563458, -0.00000000745695921,
+ +0.00000000204105478, -0.00000000059502230,
+ +0.00000000018322967, -0.00000000005920506,
+ +0.00000000001996517, -0.00000000000699511,
+ +0.00000000000253686, -0.00000000000094929,
+ +0.00000000000036552, -0.00000000000014449,
+ +0.00000000000005851, -0.00000000000002423,
+ +0.00000000000001025, -0.00000000000000442,
+ +0.00000000000000194, -0.00000000000000087,
+ +0.00000000000000039, -0.00000000000000018,
+ +0.00000000000000008};
+
+ static const double q[25] = {
+ +0.98604065696238260, -0.01347173820829521,
+ +0.00045329284116523, -0.00003067288651655,
+ +0.00000313199197601, -0.00000042110196496,
+ +0.00000006907244830, -0.00000001318321290,
+ +0.00000000283697433, -0.00000000067329234,
+ +0.00000000017339687, -0.00000000004786939,
+ +0.00000000001403235, -0.00000000000433496,
+ +0.00000000000140273, -0.00000000000047306,
+ +0.00000000000016558, -0.00000000000005994,
+ +0.00000000000002237, -0.00000000000000859,
+ +0.00000000000000338, -0.00000000000000136,
+ +0.00000000000000056, -0.00000000000000024,
+ +0.00000000000000010};
+
+ double h;
+ if (std::abs(x) <= 8) {
+ double y = r8*x;
+ h = 2*y*y-1;
+ double alfa = h+h;
+ double b0 = 0;
+ double b1 = 0;
+ double b2 = 0;
+ for (int i = 15; i >= 0; --i) {
+ b0 = s[i]+alfa*b1-b2;
+ b2 = b1;
+ b1 = b0;
+ }
+ h = y*(b0-b2);
+ } else {
+ double r = 1/x;
+ h = 128*r*r-1;
+ double alfa = h+h;
+ double b0 = 0;
+ double b1 = 0;
+ double b2 = 0;
+ for (int i = 28; i >= 0; --i) {
+ b0 = p[i]+alfa*b1-b2;
+ b2 = b1;
+ b1 = b0;
+ }
+ double pp = b0-h*b2;
+ b1 = 0;
+ b2 = 0;
+ for (int i = 24; i >= 0; --i) {
+ b0 = q[i]+alfa*b1-b2;
+ b2 = b1;
+ b1 = b0;
+ }
+ h = (x > 0 ? pih : -pih)-r*(r*pp*std::sin(x)+(b0-h*b2)*std::cos(x));
+ }
+ return h;
+}
+
+// Real part of the cosine integral
+// Translated from CERNLIB COSINT (C336) by B. List 29.4.2010
+
+double cosint(double x) {
+
+ static const double z1 = 1, r8 = z1/8, r32 = z1/32;
+
+ static const double ce = 0.57721566490153286;
+ static const double pih = PI/2;
+
+ static const double c[16] = {
+ +1.94054914648355493, +0.94134091328652134,
+ -0.57984503429299276, +0.30915720111592713,
+ -0.09161017922077134, +0.01644374075154625,
+ -0.00197130919521641, +0.00016925388508350,
+ -0.00001093932957311, +0.00000055223857484,
+ -0.00000002239949331, +0.00000000074653325,
+ -0.00000000002081833, +0.00000000000049312,
+ -0.00000000000001005, +0.00000000000000018};
+
+ static const double p[29] = {
+ +0.96074783975203596, -0.03711389621239806,
+ +0.00194143988899190, -0.00017165988425147,
+ +0.00002112637753231, -0.00000327163256712,
+ +0.00000060069211615, -0.00000012586794403,
+ +0.00000002932563458, -0.00000000745695921,
+ +0.00000000204105478, -0.00000000059502230,
+ +0.00000000018322967, -0.00000000005920506,
+ +0.00000000001996517, -0.00000000000699511,
+ +0.00000000000253686, -0.00000000000094929,
+ +0.00000000000036552, -0.00000000000014449,
+ +0.00000000000005851, -0.00000000000002423,
+ +0.00000000000001025, -0.00000000000000442,
+ +0.00000000000000194, -0.00000000000000087,
+ +0.00000000000000039, -0.00000000000000018,
+ +0.00000000000000008};
+
+ static const double q[25] = {
+ +0.98604065696238260, -0.01347173820829521,
+ +0.00045329284116523, -0.00003067288651655,
+ +0.00000313199197601, -0.00000042110196496,
+ +0.00000006907244830, -0.00000001318321290,
+ +0.00000000283697433, -0.00000000067329234,
+ +0.00000000017339687, -0.00000000004786939,
+ +0.00000000001403235, -0.00000000000433496,
+ +0.00000000000140273, -0.00000000000047306,
+ +0.00000000000016558, -0.00000000000005994,
+ +0.00000000000002237, -0.00000000000000859,
+ +0.00000000000000338, -0.00000000000000136,
+ +0.00000000000000056, -0.00000000000000024,
+ +0.00000000000000010};
+
+ double h = 0;
+ if(x == 0) {
+ h = - std::numeric_limits<double>::infinity();
+ } else if (std::abs(x) <= 8) {
+ h = r32*x*x-1;
+ double alfa = h+h;
+ double b0 = 0;
+ double b1 = 0;
+ double b2 = 0;
+ for (int i = 15; i >= 0; --i) {
+ b0 = c[i]+alfa*b1-b2;
+ b2 = b1;
+ b1 = b0;
+ }
+ h = ce+std::log(std::abs(x))-b0+h*b2;
+ } else {
+ double r = 1/x;
+ h = 128*r*r-1;
+ double alfa = h+h;
+ double b0 = 0;
+ double b1 = 0;
+ double b2 = 0;
+ for (int i = 28; i >= 0; --i) {
+ b0 = p[i]+alfa*b1-b2;
+ b2 = b1;
+ b1 = b0;
+ }
+ double pp = b0-h*b2;
+ b1 = 0;
+ b2 = 0;
+ for (int i = 24; i >= 0; --i) {
+ b0 = q[i]+alfa*b1-b2;
+ b2 = b1;
+ b1 = b0;
+ }
+ h = r*((b0-h*b2)*std::sin(x)-r*pp*cos(x));
+ }
+ return h;
+}
+
diff --git a/test/stressMathCore.cxx b/test/stressMathCore.cxx
index 2b0ff0649270617244cecfb6949824079a884697..baadb6665ad0312f7c9a839f3fbb664cf4ebffb0 100644
--- a/test/stressMathCore.cxx
+++ b/test/stressMathCore.cxx
@@ -56,7 +56,6 @@
#ifndef __CINT__
-
#include "Math/DistFuncMathCore.h"
#ifdef USE_MATHMORE
#include "Math/DistMathMore.h"
@@ -357,6 +356,7 @@ typedef StatFunction<F2,F1,2> Dist_gaussian;
typedef StatFunction<F3,F2,3> Dist_lognormal;
typedef StatFunction<F2,F1,2> Dist_tdistribution;
typedef StatFunction<F2,F1,2> Dist_exponential;
+typedef StatFunction<F2,F1,2> Dist_landau;
typedef StatFunction<F3,F2,3> Dist_uniform;
@@ -502,6 +502,21 @@ int testStatFunctions(int nfunc = 100 ) {
distc.ScaleTol2(100);
iret |= distc.Test(0.,5.,0.,1.,true);
}
+ {
+ PrintTest("Landau distribution");
+ CREATE_DIST(landau);
+ dist.SetParameters( 2);
+ // Landau is not very precise (put prec at 10-6)
+ // as indicated in Landau paper (
+ dist.ScaleTol1(10000);
+ dist.ScaleTol2(10000000000);
+ iret |= dist.Test(-1,10,-10.,1.E10);
+ CREATE_DIST_C(landau);
+ distc.SetParameters( 2);
+ distc.ScaleTol1(10000);
+ distc.ScaleTol2(10000000000);
+ iret |= distc.Test(-1,10,-10.,1.E10,true);
+ }
{
PrintTest("Uniform distribution");
diff --git a/test/stressMathMore.cxx b/test/stressMathMore.cxx
index 0b8903a095817aaa2d4cb6bb19a5a393f8d7be69..fcb276652a8ef0005d0449fd36e6604983521b53 100644
--- a/test/stressMathMore.cxx
+++ b/test/stressMathMore.cxx
@@ -270,9 +270,9 @@ int StatFunction::TestIntegral(IntegrationOneDim::Type algoType = IntegrationOne
double scale = std::max( fScaleIg * err / std::numeric_limits<double>::epsilon(), 1.);
r |= compare("test integral", q1, q2, scale );
if (r && debug) {
- std::cout << "Failed test for x = " << v1 << " p = ";
+ std::cout << "Failed test for x = " << v1 << " q1= " << q1 << " q2= " << q2 << " p = ";
for (int j = 0; j < NPAR; ++j) std::cout << fParams[j] << "\t";
- std::cout << "ig error is " << err << std::endl;
+ std::cout << "ig error is " << err << " status " << ig.Status() << std::endl;
}
iret |= r;