ROOT  6.06/09
Reference Guide
testSpecFunc.cxx
Go to the documentation of this file.
1 
2 /* MathCore/tests/test_SpecFunc.cpp
3  *
4  * Copyright (C) 2004 Andras Zsenei
5  *
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  */
21 
22 
23 /**
24 
25 Test file for the special functions implemented in MathMore. For
26 the moment nothing exceptional. Evaluates the functions and checks
27 if the value is right against values copied from the GSL tests.
28 
29 */
30 
31 
32 
33 
34 
35 #include <iostream>
36 #include <iomanip>
37 #include <string>
38 #include <limits>
39 
40 #include "Math/SpecFuncMathMore.h"
41 #ifndef NO_MATHCORE
42 #include "Math/SpecFuncMathCore.h"
43 #endif
44 
45 #ifdef CHECK_WITH_GSL
46 #include <gsl/gsl_sf.h>
47 #endif
48 
49 #ifndef PI
50 #define PI 3.14159265358979323846264338328 /* pi */
51 #endif
52 
53 
54 int compare( const std::string & name, double v1, double v2, double scale = 2.0) {
55  // ntest = ntest + 1;
56 
57  std::cout << std::setw(50) << std::left << name << ":\t";
58 
59  // numerical double limit for epsilon
60  double eps = scale* std::numeric_limits<double>::epsilon();
61  int iret = 0;
62  double delta = v2 - v1;
63  double d = 0;
64  if (delta < 0 ) delta = - delta;
65  if (v1 == 0 || v2 == 0) {
66  if (delta > eps ) {
67  iret = 1;
68  }
69  }
70  // skip case v1 or v2 is infinity
71  else {
72  d = v1;
73 
74  if ( v1 < 0) d = -d;
75  // add also case when delta is small by default
76  if ( delta/d > eps && delta > eps )
77  iret = 1;
78  }
79 
80  if (iret == 0)
81  std::cout <<" OK" << std::endl;
82  else {
83  std::cout <<" FAILED " << std::endl;
84  int pr = std::cout.precision (18);
85  std::cout << "\nDiscrepancy in " << name << "() :\n " << v1 << " != " << v2 << " discr = " << int(delta/d/eps)
86  << " (Allowed discrepancy is " << eps << ")\n\n";
87  std::cout.precision (pr);
88  //nfail = nfail + 1;
89  }
90  return iret;
91 }
92 
93 
94 
95 // void showDiff(std::string name, double calculatedValue, double expectedValue, double scale = 1.0) {
96 
97 // compare(calculatedValue, expectedValue, name, scale)
98 
99 // std::cout << name << calculatedValue << " expected value: " << expectedValue;
100 // int prec = std::cout.precision();
101 // std::cout.precision(5);
102 // std::cout << " diff: " << (calculatedValue-expectedValue) << " reldiff: " <<
103 // (calculatedValue-expectedValue)/expectedValue << std::endl;
104 // std::cout.precision(prec);
105 
106 // }
107 
108 
109 
111 
112  using namespace ROOT::Math;
113 
114  int iret = 0;
115 
116  std::cout.precision(20);
117 
118 #ifndef NO_MATHCORE
119 
120  // explicit put namespace to be sure to use right ones
121 
122  iret |= compare("tgamma(9.0) ", ROOT::Math::tgamma(9.0), 40320.0, 4);
123 
124  iret |= compare("lgamma(0.1) ", ROOT::Math::lgamma(0.1), 2.252712651734205, 4);
125 
126  iret |= compare("inc_gamma(1,0.001) ", ROOT::Math::inc_gamma(1.0,0.001), 0.0009995001666250083319, 1);
127 
128  // increase tolerance when using Cephes (test values are correctly checked with Mathematica
129  // GSL was more precise in this case
130  iret |= compare("inc_gamma(100,99) ", ROOT::Math::inc_gamma(100.,99.), 0.4733043303994607, 100);
131 
132  iret |= compare("inc_gamma_c(100,99) ", ROOT::Math::inc_gamma_c(100.,99.), 0.5266956696005394, 100);
133 
134  // need to increase here by a further factor of 5 for Windows
135  iret |= compare("inc_gamma_c(1000,1000.1) ", ROOT::Math::inc_gamma_c(1000.,1000.1), 0.4945333598559338247, 5000);
136 
137  iret |= compare("erf(0.5) ", ROOT::Math::erf(0.5), 0.5204998778130465377);
138 
139  iret |= compare("erfc(-1.0) ", ROOT::Math::erfc(-1.0), 1.8427007929497148693);
140 
141  iret |= compare("beta(1.0, 5.0) ", ROOT::Math::beta(1.0, 5.0), 0.2);
142 
143  iret |= compare("inc_beta(1,1,1) ", ROOT::Math::inc_beta(1.0, 1.0, 1.0), 1.0);
144 
145  iret |= compare("inc_beta(0.5,0.1,1.0) ", ROOT::Math::inc_beta(0.5, 0.1, 1.0), 0.9330329915368074160 );
146 
147 
148 #endif
149 
150  iret |= compare("assoc_laguerre(4, 2, 0.5) ", assoc_laguerre(4, 2, 0.5), 6.752604166666666667,8);
151 
152  iret |= compare("assoc_legendre(10, 1, -0.5) ", assoc_legendre(10, 1, -0.5), 2.0066877394361256516);
153 
154  iret |= compare("comp_ellint_1(0.50) ", comp_ellint_1(0.50), 1.6857503548125960429);
155 
156  iret |= compare("comp_ellint_2(0.50) ", comp_ellint_2(0.50), 1.4674622093394271555);
157 
158  iret |= compare("comp_ellint_3(0.5, 0.5) ", comp_ellint_3(0.5, 0.5), 2.41367150420119, 16);
159 
160  iret |= compare("conf_hyperg(1, 1.5, 1) ", conf_hyperg(1, 1.5, 1), 2.0300784692787049755);
161 
162  iret |= compare("cyl_bessel_i(1.0, 1.0) ", cyl_bessel_i(1.0, 1.0), 0.5651591039924850272);
163 
164  iret |= compare("cyl_bessel_j(0.75, 1.0) ", cyl_bessel_j(0.75, 1.0), 0.5586524932048917478, 16);
165 
166  iret |= compare("cyl_bessel_k(1.0, 1.0) ", cyl_bessel_k(1.0, 1.0), 0.6019072301972345747);
167 
168  iret |= compare("cyl_neumann(0.75, 1.0) ", cyl_neumann(0.75, 1.0), -0.6218694174429746383 );
169 
170  iret |= compare("ellint_1(0.50, PI/3.0) ", ellint_1(0.50, PI/3.0), 1.0895506700518854093);
171 
172  iret |= compare("ellint_2(0.50, PI/3.0) ", ellint_2(0.50, PI/3.0), 1.0075555551444720293);
173 
174  iret |= compare("ellint_3(-0.50, 0.5, PI/3.0) ", ellint_3(-0.50, 0.5, PI/3.0), 0.9570574331323584890);
175 
176  iret |= compare("expint(1.0) ", expint(1.0), 1.8951178163559367555);
177 
178  // std::cout << "Hermite polynomials: to do!" << std::endl;
179 
180  iret |= compare("hyperg(8, -8, 1, 0.5) ", hyperg(8, -8, 1, 0.5), 0.13671875);
181 
182  iret |= compare("laguerre(4, 1.) ", laguerre(4, 1.), -0.6250); // need to find more precise value
183 
184  iret |= compare("legendre(10, -0.5) ", legendre(10, -0.5), -0.1882286071777345);
185 
186  iret |= compare("riemann_zeta(-0.5) ", riemann_zeta(-0.5), -0.207886224977354566017307, 16);
187 
188  iret |= compare("sph_bessel(1, 10.0) ", sph_bessel(1, 10.0), 0.07846694179875154709000);
189 
190  iret |= compare("sph_legendre(3, 1, PI/2.) ", sph_legendre(3, 1, PI/2.), 0.323180184114150653007);
191 
192  iret |= compare("sph_neumann(0, 1.0) ", sph_neumann(0, 1.0), -0.54030230586813972);
193 
194  iret |= compare("airy_Ai(-0.5) ", airy_Ai(-0.5), 0.475728091610539583); // wolfram alpha: 0.47572809161053958880
195 
196  iret |= compare("airy_Bi(0.5) ", airy_Bi(0.5), 0.854277043103155553); // wolfram alpha: 0.85427704310315549330
197 
198  iret |= compare("airy_Ai_deriv(-2) ", airy_Ai_deriv(-2), 0.618259020741691145); // wolfram alpha: 0.61825902074169104141
199 
200  iret |= compare("airy_Bi_deriv(-3) ", airy_Bi_deriv(-3), -0.675611222685258639); // wolfram alpha: -0.67561122268525853767
201 
202  iret |= compare("airy_zero_Ai(2) ", airy_zero_Ai(2), -4.08794944413097028, 8); // mathworld: -4.08795
203 
204  iret |= compare("airy_zero_Bi(2) ", airy_zero_Bi(2), -3.27109330283635291, 8); // mathworld: -3.27109
205 
206  iret |= compare("airy_zero_Ai_deriv(2) ", airy_zero_Ai_deriv(2), -3.24819758217983656, 8);
207 
208  iret |= compare("airy_zero_Bi_deriv(2) ", airy_zero_Bi_deriv(2), -4.07315508907182799, 8);
209 
210  if (iret != 0) {
211  std::cout << "\n\nError: Special Functions Test FAILED !!!!!" << std::endl;
212  }
213  return iret;
214 
215 }
216 
217 void getGSLErrors() {
218 
219 #ifdef CHECK_WITH_GSL
220  gsl_sf_result r;
221  int iret;
222 
223  iret = gsl_sf_ellint_P_e(PI/2.0, 0.5, -0.5, GSL_PREC_DOUBLE, &r);
224  std::cout << "comp_ellint_3(0.50, 0.5) : " << r.val << " err: " << r.err << " iret: " << iret << std::endl;
225 
226  iret = gsl_sf_ellint_P_e(PI/3.0, 0.5, 0.5, GSL_PREC_DOUBLE, &r);
227  std::cout << "ellint_3(0.50, 0.5, PI/3.0) : " << r.val << " err: " << r.err << " iret: " << iret << std::endl;
228 
229  iret = gsl_sf_zeta_e(-0.5, &r);
230  std::cout << "riemann_zeta(-0.5) : " << r.val << " err: " << r.err << " iret: " << iret << std::endl;
231 #endif
232 
233 
234 }
235 
236 
237 int main() {
238 
239  int iret = 0;
240  iret |= testSpecFunc();
241 
242  getGSLErrors();
243  return iret;
244 }
245 
246 
double erf(double x)
Error function encountered in integrating the normal distribution.
double airy_zero_Bi(unsigned int s)
Calculates the zeroes of the Airy function Bi.
const Double_t * v1
Definition: TArcBall.cxx:33
double sph_legendre(unsigned l, unsigned m, double theta)
Computes the spherical (normalized) associated Legendre polynomials, or spherical harmonic without az...
double cyl_bessel_j(double nu, double x)
Calculates the (cylindrical) Bessel functions of the first kind (also called regular (cylindrical) Be...
double inc_beta(double x, double a, double b)
Calculates the normalized (regularized) incomplete beta function.
double airy_zero_Ai_deriv(unsigned int s)
Calculates the zeroes of the derivative of the Airy function Ai.
double ellint_2(double k, double phi)
Calculates the complete elliptic integral of the second kind.
double tgamma(double x)
The gamma function is defined to be the extension of the factorial to real numbers.
double legendre(unsigned l, double x)
Calculates the Legendre polynomials.
double beta(double x, double y)
Calculates the beta function.
double erfc(double x)
Complementary error function.
double inc_gamma_c(double a, double x)
Calculates the normalized (regularized) upper incomplete gamma function (upper integral) ...
double comp_ellint_1(double k)
Calculates the complete elliptic integral of the first kind.
double sph_bessel(unsigned n, double x)
Calculates the spherical Bessel functions of the first kind (also called regular spherical Bessel fun...
double conf_hyperg(double a, double b, double z)
Calculates the confluent hypergeometric functions of the first kind.
double hyperg(double a, double b, double c, double x)
Calculates Gauss' hypergeometric function.
double airy_zero_Bi_deriv(unsigned int s)
Calculates the zeroes of the derivative of the Airy function Bi.
double expint(double x)
Calculates the exponential integral.
double ellint_3(double n, double k, double phi)
Calculates the complete elliptic integral of the third kind.
double riemann_zeta(double x)
Calculates the Riemann zeta function.
double comp_ellint_2(double k)
Calculates the complete elliptic integral of the second kind.
ROOT::R::TRInterface & r
Definition: Object.C:4
double assoc_laguerre(unsigned n, double m, double x)
Computes the generalized Laguerre polynomials for .
double assoc_legendre(unsigned l, unsigned m, double x)
Computes the associated Legendre polynomials.
double airy_zero_Ai(unsigned int s)
Calculates the zeroes of the Airy function Ai.
double airy_Bi(double x)
Calculates the Airy function Bi.
double laguerre(unsigned n, double x)
Calculates the Laguerre polynomials.
REAL epsilon
Definition: triangle.c:617
double airy_Bi_deriv(double x)
Calculates the derivative of the Airy function Bi.
double comp_ellint_3(double n, double k)
Calculates the complete elliptic integral of the third kind.
int main()
double airy_Ai_deriv(double x)
Calculates the derivative of the Airy function Ai.
double ellint_1(double k, double phi)
Calculates the incomplete elliptic integral of the first kind.
int compare(const std::string &name, double v1, double v2, double scale=2.0)
double cyl_neumann(double nu, double x)
Calculates the (cylindrical) Bessel functions of the second kind (also called irregular (cylindrical)...
#define name(a, b)
Definition: linkTestLib0.cpp:5
#define PI
Test file for the special functions implemented in MathMore.
void getGSLErrors()
double sph_neumann(unsigned n, double x)
Calculates the spherical Bessel functions of the second kind (also called irregular spherical Bessel ...
double cyl_bessel_i(double nu, double x)
Calculates the modified Bessel function of the first kind (also called regular modified (cylindrical)...
double lgamma(double x)
Calculates the logarithm of the gamma function.
double airy_Ai(double x)
Calculates the Airy function Ai.
double inc_gamma(double a, double x)
Calculates the normalized (regularized) lower incomplete gamma function (lower integral) ...
int testSpecFunc()
double cyl_bessel_k(double nu, double x)
Calculates the modified Bessel functions of the second kind (also called irregular modified (cylindri...