Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
VavilovAccurate.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: B. List 29.4.2010
3 
4 
5  /**********************************************************************
6  * *
7  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
8  * *
9  * This library is free software; you can redistribute it and/or *
10  * modify it under the terms of the GNU General Public License *
11  * as published by the Free Software Foundation; either version 2 *
12  * of the License, or (at your option) any later version. *
13  * *
14  * This library is distributed in the hope that it will be useful, *
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
17  * General Public License for more details. *
18  * *
19  * You should have received a copy of the GNU General Public License *
20  * along with this library (see file COPYING); if not, write *
21  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
22  * 330, Boston, MA 02111-1307 USA, or contact the author. *
23  * *
24  **********************************************************************/
25 
26 // Implementation file for class VavilovAccurate
27 //
28 // Created by: blist at Thu Apr 29 11:19:00 2010
29 //
30 // Last update: Thu Apr 29 11:19:00 2010
31 //
32 
33 
34 #include "Math/VavilovAccurate.h"
35 #include "Math/SpecFuncMathCore.h"
36 #include "Math/SpecFuncMathMore.h"
37 #include "Math/QuantFuncMathCore.h"
38 
39 #include <cassert>
40 #include <iostream>
41 #include <iomanip>
42 #include <cmath>
43 #include <limits>
44 
45 
46 namespace ROOT {
47 namespace Math {
48 
49 VavilovAccurate *VavilovAccurate::fgInstance = 0;
50 
51 
52 VavilovAccurate::VavilovAccurate(double kappa, double beta2, double epsilonPM, double epsilon)
53 {
54  Set (kappa, beta2, epsilonPM, epsilon);
55 }
56 
57 
58 VavilovAccurate::~VavilovAccurate()
59 {
60  // desctructor (clean up resources)
61 }
62 
63 void VavilovAccurate::SetKappaBeta2 (double kappa, double beta2) {
64  Set (kappa, beta2);
65 }
66 
67 void VavilovAccurate::Set(double kappa, double beta2, double epsilonPM, double epsilon) {
68  // Method described in
69  // B. Schorr, Programs for the Landau and the Vavilov distributions and the corresponding random numbers,
70  // <A HREF="http://dx.doi.org/10.1016/0010-4655(74)90091-5">Computer Phys. Comm. 7 (1974) 215-224</A>.
71  fQuantileInit = false;
72 
73  fKappa = kappa;
74  fBeta2 = beta2;
75  fEpsilonPM = epsilonPM; // epsilon_+ = epsilon_-: determines support (T0, T1)
76  fEpsilon = epsilon;
77 
78  static const double eu = 0.577215664901532860606; // Euler's constant
79  static const double pi2 = 6.28318530717958647693, // 2pi
80  rpi = 0.318309886183790671538, // 1/pi
81  pih = 1.57079632679489661923; // pi/2
82  double h1 = -std::log(fEpsilon)-1.59631259113885503887; // -ln(fEpsilon) + ln(2/pi**2)
83  double deltaEpsilon = 0.001;
84  static const double logdeltaEpsilon = -std::log(deltaEpsilon); // 3 ln 10 = -ln(.001);
85  double logEpsilonPM = std::log(fEpsilonPM);
86  static const double eps = 1e-5; // accuracy of root finding for x0
87 
88  double xp[9] = {0,
89  9.29, 2.47, 0.89, 0.36, 0.15, 0.07, 0.03, 0.02};
90  double xq[7] = {0,
91  0.012, 0.03, 0.08, 0.26, 0.87, 3.83};
92 
93  if (kappa < 0.001) {
94  std::cerr << "VavilovAccurate::Set: kappa = " << kappa << " - out of range" << std::endl;
95  if (kappa < 0.001) kappa = 0.001;
96  }
97  if (beta2 < 0 || beta2 > 1) {
98  std::cerr << "VavilovAccurate::Set: beta2 = " << beta2 << " - out of range" << std::endl;
99  if (beta2 < 0) beta2 = -beta2;
100  if (beta2 > 1) beta2 = 1;
101  }
102 
103  // approximation of x_-
104  fH[5] = 1-beta2*(1-eu)-logEpsilonPM/kappa; // eq. 3.9
105  fH[6] = beta2;
106  fH[7] = 1-beta2;
107  double h4 = logEpsilonPM/kappa-(1+beta2*eu);
108  double logKappa = std::log(kappa);
109  double kappaInv = 1/kappa;
110  // Calculate T0 from Eq. (3.6), using x_- = fH[5]
111 // double e1h5 = (fH[5] > 40 ) ? 0 : -ROOT::Math::expint (-fH[5]);
112 // fT0 = (h4-fH[5]*logKappa-(fH[5]+beta2)*(std::log(fH[5])+e1h5)+std::exp(-fH[5]))/fH[5];
113  fT0 = (h4-fH[5]*logKappa-(fH[5]+beta2)*E1plLog(fH[5])+std::exp(-fH[5]))/fH[5];
114  int lp = 1;
115  while (lp < 9 && kappa < xp[lp]) ++lp;
116  int lq = 1;
117  while (lq < 7 && kappa >= xq[lq]) ++lq;
118  // Solve eq. 3.7 to get fH[0] = x_+
119  double delta = 0;
120  int ifail = 0;
121  do {
122  ifail = Rzero(-lp-0.5-delta,lq-7.5+delta,fH[0],eps,1000,&ROOT::Math::VavilovAccurate::G116f2);
123  delta += 0.5;
124  } while (ifail == 2);
125 
126  double q = 1/fH[0];
127  // Calculate T1 from Eq. (3.6)
128 // double e1h0 = (fH[0] > 40 ) ? 0 : -ROOT::Math::expint (-fH[0]);
129 // fT1 = h4*q-logKappa-(1+beta2*q)*(std::log(std::fabs(fH[0]))+e1h0)+std::exp(-fH[0])*q;
130  fT1 = h4*q-logKappa-(1+beta2*q)*E1plLog(fH[0])+std::exp(-fH[0])*q;
131 
132  fT = fT1-fT0; // Eq. (2.5)
133  fOmega = pi2/fT; // Eq. (2.5)
134  fH[1] = kappa*(2+beta2*eu)+h1;
135  if(kappa >= 0.07) fH[1] += logdeltaEpsilon; // reduce fEpsilon by a factor .001 for large kappa
136  fH[2] = beta2*kappa;
137  fH[3] = kappaInv*fOmega;
138  fH[4] = pih*fOmega;
139 
140  // Solve log(eq. (4.10)) to get fX0 = N
141  ifail = Rzero(5.,MAXTERMS,fX0,eps,1000,&ROOT::Math::VavilovAccurate::G116f1);
142 // if (ifail) {
143 // std::cerr << "Rzero failed for x0: ifail=" << ifail << ", kappa=" << kappa << ", beta2=" << beta2 << std::endl;
144 // std::cerr << "G116f1(" << 5. << ")=" << G116f1(5.) << ", G116f1(" << MAXTERMS << ")=" << G116f1(MAXTERMS) << std::endl;
145 // std::cerr << "fH[0]=" << fH[0] << ", fH[1]=" << fH[1] << ", fH[2]=" << fH[2] << ", fH[3]=" << fH[3] << ", fH[4]=" << fH[4] << std::endl;
146 // std::cerr << "fH[5]=" << fH[5] << ", fH[6]=" << fH[6] << ", fH[7]=" << fH[7] << std::endl;
147 // std::cerr << "fT0=" << fT0 << ", fT1=" << fT1 << std::endl;
148 // std::cerr << "G116f2(" << fH[0] << ")=" << G116f2(fH[0]) << std::endl;
149 // }
150  if (ifail == 2) {
151  fX0 = (G116f1(5) > G116f1(MAXTERMS)) ? MAXTERMS : 5;
152  }
153  if (fX0 < 5) fX0 = 5;
154  else if (fX0 > MAXTERMS) fX0 = MAXTERMS;
155  int n = int(fX0+1);
156  // logKappa=log(kappa)
157  double d = rpi*std::exp(kappa*(1+beta2*(eu-logKappa)));
158  fA_pdf[n] = rpi*fOmega;
159  fA_cdf[n] = 0;
160  q = -1;
161  double q2 = 2;
162  for (int k = 1; k < n; ++k) {
163  int l = n-k;
164  double x = fOmega*k;
165  double x1 = kappaInv*x;
166  double c1 = std::log(x)-ROOT::Math::cosint(x1);
167  double c2 = ROOT::Math::sinint(x1);
168  double c3 = std::sin(x1);
169  double c4 = std::cos(x1);
170  double xf1 = kappa*(beta2*c1-c4)-x*c2;
171  double xf2 = x*(c1 + fT0) + kappa*(c3+beta2*c2);
172  double d1 = q*d*fOmega*std::exp(xf1);
173  double s = std::sin(xf2);
174  double c = std::cos(xf2);
175  fA_pdf[l] = d1*c;
176  fB_pdf[l] = -d1*s;
177  d1 = q*d*std::exp(xf1)/k;
178  fA_cdf[l] = d1*s;
179  fB_cdf[l] = d1*c;
180  fA_cdf[n] += q2*fA_cdf[l];
181  q = -q;
182  q2 = -q2;
183  }
184 }
185 
186 void VavilovAccurate::InitQuantile() const {
187  fQuantileInit = true;
188 
189  fNQuant = 16;
190  // for kappa<0.02: use landau_quantile as first approximation
191  if (fKappa < 0.02) return;
192  else if (fKappa < 0.05) fNQuant = 32;
193 
194  // crude approximation for the median:
195 
196  double estmedian = -4.22784335098467134e-01-std::log(fKappa)-fBeta2;
197  if (estmedian>1.3) estmedian = 1.3;
198 
199  // distribute test values evenly below and above the median
200  for (int i = 1; i < fNQuant/2; ++i) {
201  double x = fT0 + i*(estmedian-fT0)/(fNQuant/2);
202  fQuant[i] = Cdf(x);
203  fLambda[i] = x;
204  }
205  for (int i = fNQuant/2; i < fNQuant-1; ++i) {
206  double x = estmedian + (i-fNQuant/2)*(fT1-estmedian)/(fNQuant/2-1);
207  fQuant[i] = Cdf(x);
208  fLambda[i] = x;
209  }
210 
211  fQuant[0] = 0;
212  fLambda[0] = fT0;
213  fQuant[fNQuant-1] = 1;
214  fLambda[fNQuant-1] = fT1;
215 
216 }
217 
218 double VavilovAccurate::Pdf (double x) const {
219 
220  static const double pi = 3.14159265358979323846; // pi
221 
222  int n = int(fX0);
223  double f;
224  if (x < fT0) {
225  f = 0;
226  } else if (x <= fT1) {
227  double y = x-fT0;
228  double u = fOmega*y-pi;
229  double cof = 2*cos(u);
230  double a1 = 0;
231  double a0 = fA_pdf[1];
232  double a2 = 0;
233  for (int k = 2; k <= n+1; ++k) {
234  a2 = a1;
235  a1 = a0;
236  a0 = fA_pdf[k]+cof*a1-a2;
237  }
238  double b1 = 0;
239  double b0 = fB_pdf[1];
240  for (int k = 2; k <= n; ++k) {
241  double b2 = b1;
242  b1 = b0;
243  b0 = fB_pdf[k]+cof*b1-b2;
244  }
245  f = 0.5*(a0-a2)+b0*sin(u);
246  } else {
247  f = 0;
248  }
249  return f;
250 }
251 
252 double VavilovAccurate::Pdf (double x, double kappa, double beta2) {
253  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
254  return Pdf (x);
255 }
256 
257 double VavilovAccurate::Cdf (double x) const {
258 
259  static const double pi = 3.14159265358979323846; // pi
260 
261  int n = int(fX0);
262  double f;
263  if (x < fT0) {
264  f = 0;
265  } else if (x <= fT1) {
266  double y = x-fT0;
267  double u = fOmega*y-pi;
268  double cof = 2*cos(u);
269  double a1 = 0;
270  double a0 = fA_cdf[1];
271  double a2 = 0;
272  for (int k = 2; k <= n+1; ++k) {
273  a2 = a1;
274  a1 = a0;
275  a0 = fA_cdf[k]+cof*a1-a2;
276  }
277  double b1 = 0;
278  double b0 = fB_cdf[1];
279  for (int k = 2; k <= n; ++k) {
280  double b2 = b1;
281  b1 = b0;
282  b0 = fB_cdf[k]+cof*b1-b2;
283  }
284  f = 0.5*(a0-a2)+b0*sin(u);
285  f += y/fT;
286  } else {
287  f = 1;
288  }
289  return f;
290 }
291 
292 double VavilovAccurate::Cdf (double x, double kappa, double beta2) {
293  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
294  return Cdf (x);
295 }
296 
297 double VavilovAccurate::Cdf_c (double x) const {
298 
299  static const double pi = 3.14159265358979323846; // pi
300 
301  int n = int(fX0);
302  double f;
303  if (x < fT0) {
304  f = 1;
305  } else if (x <= fT1) {
306  double y = fT1-x;
307  double u = fOmega*y-pi;
308  double cof = 2*cos(u);
309  double a1 = 0;
310  double a0 = fA_cdf[1];
311  double a2 = 0;
312  for (int k = 2; k <= n+1; ++k) {
313  a2 = a1;
314  a1 = a0;
315  a0 = fA_cdf[k]+cof*a1-a2;
316  }
317  double b1 = 0;
318  double b0 = fB_cdf[1];
319  for (int k = 2; k <= n; ++k) {
320  double b2 = b1;
321  b1 = b0;
322  b0 = fB_cdf[k]+cof*b1-b2;
323  }
324  f = -0.5*(a0-a2)+b0*sin(u);
325  f += y/fT;
326  } else {
327  f = 0;
328  }
329  return f;
330 }
331 
332 double VavilovAccurate::Cdf_c (double x, double kappa, double beta2) {
333  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
334  return Cdf_c (x);
335 }
336 
337 double VavilovAccurate::Quantile (double z) const {
338  if (z < 0 || z > 1) return std::numeric_limits<double>::signaling_NaN();
339 
340  if (!fQuantileInit) InitQuantile();
341 
342  double x;
343  if (fKappa < 0.02) {
344  x = ROOT::Math::landau_quantile (z*(1-2*fEpsilonPM) + fEpsilonPM);
345  if (x < fT0+5*fEpsilon) x = fT0+5*fEpsilon;
346  else if (x > fT1-10*fEpsilon) x = fT1-10*fEpsilon;
347  }
348  else {
349  // yes, I know what a binary search is, but linear search is faster for small n!
350  int i = 1;
351  while (z > fQuant[i]) ++i;
352  assert (i < fNQuant);
353 
354  assert (i >= 1);
355  assert (i < fNQuant);
356 
357  // initial solution
358  double f = (z-fQuant[i-1])/(fQuant[i]-fQuant[i-1]);
359  assert (f >= 0);
360  assert (f <= 1);
361  assert (fQuant[i] > fQuant[i-1]);
362 
363  x = f*fLambda[i] + (1-f)*fLambda[i-1];
364  }
365  if (fabs(x-fT0) < fEpsilon || fabs(x-fT1) < fEpsilon) return x;
366 
367  assert (x > fT0 && x < fT1);
368  double dx;
369  int n = 0;
370  do {
371  ++n;
372  double y = Cdf(x)-z;
373  double y1 = Pdf (x);
374  dx = - y/y1;
375  x = x + dx;
376  // protect against shooting beyond the support
377  if (x < fT0) x = 0.5*(fT0+x-dx);
378  else if (x > fT1) x = 0.5*(fT1+x-dx);
379  assert (x > fT0 && x < fT1);
380  } while (fabs(dx) > fEpsilon && n < 100);
381  return x;
382 }
383 
384 double VavilovAccurate::Quantile (double z, double kappa, double beta2) {
385  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
386  return Quantile (z);
387 }
388 
389 double VavilovAccurate::Quantile_c (double z) const {
390  if (z < 0 || z > 1) return std::numeric_limits<double>::signaling_NaN();
391 
392  if (!fQuantileInit) InitQuantile();
393 
394  double z1 = 1-z;
395 
396  double x;
397  if (fKappa < 0.02) {
398  x = ROOT::Math::landau_quantile (z1*(1-2*fEpsilonPM) + fEpsilonPM);
399  if (x < fT0+5*fEpsilon) x = fT0+5*fEpsilon;
400  else if (x > fT1-10*fEpsilon) x = fT1-10*fEpsilon;
401  }
402  else {
403  // yes, I know what a binary search is, but linear search is faster for small n!
404  int i = 1;
405  while (z1 > fQuant[i]) ++i;
406  assert (i < fNQuant);
407 
408 // int i0=0, i1=fNQuant, i;
409 // for (int it = 0; it < LOG2fNQuant; ++it) {
410 // i = (i0+i1)/2;
411 // if (z > fQuant[i]) i0 = i;
412 // else i1 = i;
413 // }
414 // assert (i1-i0 == 1);
415 
416  assert (i >= 1);
417  assert (i < fNQuant);
418 
419  // initial solution
420  double f = (z1-fQuant[i-1])/(fQuant[i]-fQuant[i-1]);
421  assert (f >= 0);
422  assert (f <= 1);
423  assert (fQuant[i] > fQuant[i-1]);
424 
425  x = f*fLambda[i] + (1-f)*fLambda[i-1];
426  }
427  if (fabs(x-fT0) < fEpsilon || fabs(x-fT1) < fEpsilon) return x;
428 
429  assert (x > fT0 && x < fT1);
430  double dx;
431  int n = 0;
432  do {
433  ++n;
434  double y = Cdf_c(x)-z;
435  double y1 = -Pdf (x);
436  dx = - y/y1;
437  x = x + dx;
438  // protect against shooting beyond the support
439  if (x < fT0) x = 0.5*(fT0+x-dx);
440  else if (x > fT1) x = 0.5*(fT1+x-dx);
441  assert (x > fT0 && x < fT1);
442  } while (fabs(dx) > fEpsilon && n < 100);
443  return x;
444 }
445 
446 double VavilovAccurate::Quantile_c (double z, double kappa, double beta2) {
447  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
448  return Quantile_c (z);
449 }
450 
451 VavilovAccurate *VavilovAccurate::GetInstance() {
452  if (!fgInstance) fgInstance = new VavilovAccurate (1, 1);
453  return fgInstance;
454 }
455 
456 VavilovAccurate *VavilovAccurate::GetInstance(double kappa, double beta2) {
457  if (!fgInstance) fgInstance = new VavilovAccurate (kappa, beta2);
458  else if (kappa != fgInstance->fKappa || beta2 != fgInstance->fBeta2) fgInstance->Set (kappa, beta2);
459  return fgInstance;
460 }
461 
462 double vavilov_accurate_pdf (double x, double kappa, double beta2) {
463  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
464  return vavilov->Pdf (x);
465 }
466 
467 double vavilov_accurate_cdf_c (double x, double kappa, double beta2) {
468  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
469  return vavilov->Cdf_c (x);
470 }
471 
472 double vavilov_accurate_cdf (double x, double kappa, double beta2) {
473  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
474  return vavilov->Cdf (x);
475 }
476 
477 double vavilov_accurate_quantile (double z, double kappa, double beta2) {
478  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
479  return vavilov->Quantile (z);
480 }
481 
482 double vavilov_accurate_quantile_c (double z, double kappa, double beta2) {
483  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
484  return vavilov->Quantile_c (z);
485 }
486 
487 double VavilovAccurate::G116f1 (double x) const {
488  // fH[1] = kappa*(2+beta2*eu) -ln(fEpsilon) + ln(2/pi**2)
489  // fH[2] = beta2*kappa
490  // fH[3] = omwga/kappa
491  // fH[4] = pi/2 *fOmega
492  // log of Eq. (4.10)
493  return fH[1]+fH[2]*std::log(fH[3]*x)-fH[4]*x;
494 }
495 
496 double VavilovAccurate::G116f2 (double x) const {
497  // fH[5] = 1-beta2*(1-eu)-logEpsilonPM/kappa; // eq. 3.9
498  // fH[6] = beta2;
499  // fH[7] = 1-beta2;
500  // Eq. 3.7 of Schorr
501 // return fH[5]-x+fH[6]*(std::log(std::fabs(x))-ROOT::Math::expint (-x))-fH[7]*std::exp(-x);
502  return fH[5]-x+fH[6]*E1plLog(x)-fH[7]*std::exp(-x);
503 }
504 
505 int VavilovAccurate::Rzero (double a, double b, double& x0,
506  double eps, int mxf, double (VavilovAccurate::*f)(double)const) const {
507 
508  double xa, xb, fa, fb, r;
509 
510  if (a <= b) {
511  xa = a;
512  xb = b;
513  } else {
514  xa = b;
515  xb = a;
516  }
517  fa = (this->*f)(xa);
518  fb = (this->*f)(xb);
519 
520  if(fa*fb > 0) {
521  r = -2*(xb-xa);
522  x0 = 0;
523 // std::cerr << "VavilovAccurate::Rzero: fail=2, f(" << a << ")=" << (this->*f) (a)
524 // << ", f(" << b << ")=" << (this->*f) (b) << std::endl;
525  return 2;
526  }
527  int mc = 0;
528 
529  bool recalcF12 = true;
530  bool recalcFab = true;
531  bool fail = false;
532 
533 
534  double x1=0, x2=0, f1=0, f2=0, fx=0, ee=0;
535  do {
536  if (recalcF12) {
537  x0 = 0.5*(xa+xb);
538  r = x0-xa;
539  ee = eps*(std::abs(x0)+1);
540  if(r <= ee) break;
541  f1 = fa;
542  x1 = xa;
543  f2 = fb;
544  x2 = xb;
545  }
546  if (recalcFab) {
547  fx = (this->*f)(x0);
548  ++mc;
549  if(mc > mxf) {
550  fail = true;
551  break;
552  }
553  if(fx*fa > 0) {
554  xa = x0;
555  fa = fx;
556  } else {
557  xb = x0;
558  fb = fx;
559  }
560  }
561  recalcF12 = true;
562  recalcFab = true;
563 
564  double u1 = f1-f2;
565  double u2 = x1-x2;
566  double u3 = f2-fx;
567  double u4 = x2-x0;
568  if(u2 == 0 || u4 == 0) continue;
569  double f3 = fx;
570  double x3 = x0;
571  u1 = u1/u2;
572  u2 = u3/u4;
573  double ca = u1-u2;
574  double cb = (x1+x2)*u2-(x2+x0)*u1;
575  double cc = (x1-x0)*f1-x1*(ca*x1+cb);
576  if(ca == 0) {
577  if(cb == 0) continue;
578  x0 = -cc/cb;
579  } else {
580  u3 = cb/(2*ca);
581  u4 = u3*u3-cc/ca;
582  if(u4 < 0) continue;
583  x0 = -u3 + (x0+u3 >= 0 ? +1 : -1)*std::sqrt(u4);
584  }
585  if(x0 < xa || x0 > xb) continue;
586 
587  recalcF12 = false;
588 
589  r = std::abs(x0-x3) < std::abs(x0-x2) ? std::abs(x0-x3) : std::abs(x0-x2);
590  ee = eps*(std::abs(x0)+1);
591  if (r > ee) {
592  f1 = f2;
593  x1 = x2;
594  f2 = f3;
595  x2 = x3;
596  continue;
597  }
598 
599  recalcFab = false;
600 
601  fx = (this->*f) (x0);
602  if (fx == 0) break;
603  double xx, ff;
604  if (fx*fa < 0) {
605  xx = x0-ee;
606  if (xx <= xa) break;
607  ff = (this->*f)(xx);
608  fb = ff;
609  xb = xx;
610  } else {
611  xx = x0+ee;
612  if(xx >= xb) break;
613  ff = (this->*f)(xx);
614  fa = ff;
615  xa = xx;
616  }
617  if (fx*ff <= 0) break;
618  mc += 2;
619  if (mc > mxf) {
620  fail = true;
621  break;
622  }
623  f1 = f3;
624  x1 = x3;
625  f2 = fx;
626  x2 = x0;
627  x0 = xx;
628  fx = ff;
629  }
630  while (true);
631 
632  if (fail) {
633  r = -0.5*std::abs(xb-xa);
634  x0 = 0;
635  std::cerr << "VavilovAccurate::Rzero: fail=" << fail << ", f(" << a << ")=" << (this->*f) (a)
636  << ", f(" << b << ")=" << (this->*f) (b) << std::endl;
637  return 1;
638  }
639 
640  r = ee;
641  return 0;
642 }
643 
644 // Calculates log(|x|)+E_1(x)
645 double VavilovAccurate::E1plLog (double x) {
646  static const double eu = 0.577215664901532860606; // Euler's constant
647  double absx = std::fabs(x);
648  if (absx < 1E-4) {
649  return (x-0.25*x)*x-eu;
650  }
651  else if (x > 35) {
652  return log (x);
653  }
654  else if (x < -50) {
655  return -ROOT::Math::expint (-x);
656  }
657  return log(absx) -ROOT::Math::expint (-x);
658 }
659 
660 double VavilovAccurate::GetLambdaMin() const {
661  return fT0;
662 }
663 
664 double VavilovAccurate::GetLambdaMax() const {
665  return fT1;
666 }
667 
668 double VavilovAccurate::GetKappa() const {
669  return fKappa;
670 }
671 
672 double VavilovAccurate::GetBeta2() const {
673  return fBeta2;
674 }
675 
676 double VavilovAccurate::Mode() const {
677  double x = -4.22784335098467134e-01-std::log(fKappa)-fBeta2;
678  if (x>-0.223172) x = -0.223172;
679  double eps = 0.01;
680  double dx;
681 
682  do {
683  double p0 = Pdf (x - eps);
684  double p1 = Pdf (x);
685  double p2 = Pdf (x + eps);
686  double y1 = 0.5*(p2-p0)/eps;
687  double y2 = (p2-2*p1+p0)/(eps*eps);
688  dx = - y1/y2;
689  x += dx;
690  if (fabs(dx) < eps) eps = 0.1*fabs(dx);
691  } while (fabs(dx) > 1E-5);
692  return x;
693 }
694 
695 double VavilovAccurate::Mode(double kappa, double beta2) {
696  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
697  return Mode();
698 }
699 
700 double VavilovAccurate::GetEpsilonPM() const {
701  return fEpsilonPM;
702 }
703 
704 double VavilovAccurate::GetEpsilon() const {
705  return fEpsilon;
706 }
707 
708 double VavilovAccurate::GetNTerms() const {
709  return fX0;
710 }
711 
712 
713 
714 } // namespace Math
715 } // namespace ROOT