Logo ROOT   6.14/05
Reference Guide
TMath.h
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: Rene Brun, Anna Kreshuk, Eddy Offermann, Fons Rademakers 29/07/95
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2004, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #ifndef ROOT_TMath
13 #define ROOT_TMath
14 
15 #include "Rtypes.h"
16 #include "TMathBase.h"
17 
18 #include "TError.h"
19 #include <algorithm>
20 #include <limits>
21 #include <cmath>
22 
23 ////////////////////////////////////////////////////////////////////////////////
24 ///
25 /// TMath
26 ///
27 /// Encapsulate most frequently used Math functions.
28 /// NB. The basic functions Min, Max, Abs and Sign are defined
29 /// in TMathBase.
30 
31 namespace TMath {
32 
33 ////////////////////////////////////////////////////////////////////////////////
34 // Fundamental constants
35 
36 ////////////////////////////////////////////////////////////////////////////////
37 /// \f[ \pi\f]
38 constexpr Double_t Pi()
39 {
40  return 3.14159265358979323846;
41 }
42 
43 ////////////////////////////////////////////////////////////////////////////////
44 /// \f[ 2\pi\f]
45 constexpr Double_t TwoPi()
46 {
47  return 2.0 * Pi();
48 }
49 
50 ////////////////////////////////////////////////////////////////////////////////
51 /// \f[ \frac{\pi}{2} \f]
52 constexpr Double_t PiOver2()
53 {
54  return Pi() / 2.0;
55 }
56 
57 ////////////////////////////////////////////////////////////////////////////////
58 /// \f[ \frac{\pi}{4} \f]
59 constexpr Double_t PiOver4()
60 {
61  return Pi() / 4.0;
62 }
63 
64 ////////////////////////////////////////////////////////////////////////////////
65 /// \f$ \frac{1.}{\pi}\f$
66 constexpr Double_t InvPi()
67 {
68  return 1.0 / Pi();
69 }
70 
71 ////////////////////////////////////////////////////////////////////////////////
72 /// Conversion from radian to degree:
73 /// \f[ \frac{180}{\pi} \f]
74 constexpr Double_t RadToDeg()
75 {
76  return 180.0 / Pi();
77 }
78 
79 ////////////////////////////////////////////////////////////////////////////////
80 /// Conversion from degree to radian:
81 /// \f[ \frac{\pi}{180} \f]
82 constexpr Double_t DegToRad()
83 {
84  return Pi() / 180.0;
85 }
86 
87 ////////////////////////////////////////////////////////////////////////////////
88 /// \f[ \sqrt{2} \f]
89 constexpr Double_t Sqrt2()
90 {
91  return 1.4142135623730950488016887242097;
92 }
93 
94 ////////////////////////////////////////////////////////////////////////////////
95 /// Base of natural log:
96 /// \f[ e \f]
97 constexpr Double_t E()
98 {
99  return 2.71828182845904523536;
100 }
101 
102 ////////////////////////////////////////////////////////////////////////////////
103 /// Natural log of 10 (to convert log to ln)
104 constexpr Double_t Ln10()
105 {
106  return 2.30258509299404568402;
107 }
108 
109 ////////////////////////////////////////////////////////////////////////////////
110 /// Base-10 log of e (to convert ln to log)
111 constexpr Double_t LogE()
112 {
113  return 0.43429448190325182765;
114 }
115 
116 ////////////////////////////////////////////////////////////////////////////////
117 /// Velocity of light in \f$ m s^{-1} \f$
118 constexpr Double_t C()
119 {
120  return 2.99792458e8;
121 }
122 
123 ////////////////////////////////////////////////////////////////////////////////
124 /// \f$ cm s^{-1} \f$
125 constexpr Double_t Ccgs()
126 {
127  return 100.0 * C();
128 }
129 
130 ////////////////////////////////////////////////////////////////////////////////
131 /// Speed of light uncertainty.
133 {
134  return 0.0;
135 }
136 
137 ////////////////////////////////////////////////////////////////////////////////
138 /// Gravitational constant in: \f$ m^{3} kg^{-1} s^{-2} \f$
139 constexpr Double_t G()
140 {
141  return 6.673e-11;
142 }
143 
144 ////////////////////////////////////////////////////////////////////////////////
145 /// \f$ cm^{3} g^{-1} s^{-2} \f$
146 constexpr Double_t Gcgs()
147 {
148  return G() / 1000.0;
149 }
150 
151 ////////////////////////////////////////////////////////////////////////////////
152 /// Gravitational constant uncertainty.
154 {
155  return 0.010e-11;
156 }
157 
158 ////////////////////////////////////////////////////////////////////////////////
159 /// \f$ \frac{G}{\hbar C} \f$ in \f$ (GeV/c^{2})^{-2} \f$
160 constexpr Double_t GhbarC()
161 {
162  return 6.707e-39;
163 }
164 
165 ////////////////////////////////////////////////////////////////////////////////
166 /// \f$ \frac{G}{\hbar C} \f$ uncertainty.
168 {
169  return 0.010e-39;
170 }
171 
172 ////////////////////////////////////////////////////////////////////////////////
173 /// Standard acceleration of gravity in \f$ m s^{-2} \f$
174 constexpr Double_t Gn()
175 {
176  return 9.80665;
177 }
178 
179 ////////////////////////////////////////////////////////////////////////////////
180 /// Standard acceleration of gravity uncertainty.
182 {
183  return 0.0;
184 }
185 
186 ////////////////////////////////////////////////////////////////////////////////
187 /// Planck's constant in \f$ J s \f$
188 /// \f[ h \f]
189 constexpr Double_t H()
190 {
191  return 6.62606876e-34;
192 }
193 
194 ////////////////////////////////////////////////////////////////////////////////
195 /// \f$ erg s \f$
196 constexpr Double_t Hcgs()
197 {
198  return 1.0e7 * H();
199 }
200 
201 ////////////////////////////////////////////////////////////////////////////////
202 /// Planck's constant uncertainty.
204 {
205  return 0.00000052e-34;
206 }
207 
208 ////////////////////////////////////////////////////////////////////////////////
209 /// \f$ \hbar \f$ in \f$ J s \f$
210 /// \f[ \hbar = \frac{h}{2\pi} \f]
211 constexpr Double_t Hbar()
212 {
213  return 1.054571596e-34;
214 }
215 
216 ////////////////////////////////////////////////////////////////////////////////
217 /// \f$ erg s \f$
218 constexpr Double_t Hbarcgs()
219 {
220  return 1.0e7 * Hbar();
221 }
222 
223 ////////////////////////////////////////////////////////////////////////////////
224 /// \f$ \hbar \f$ uncertainty.
226 {
227  return 0.000000082e-34;
228 }
229 
230 ////////////////////////////////////////////////////////////////////////////////
231 /// \f$ hc \f$ in \f$ J m \f$
232 constexpr Double_t HC()
233 {
234  return H() * C();
235 }
236 
237 ////////////////////////////////////////////////////////////////////////////////
238 /// \f$ erg cm \f$
239 constexpr Double_t HCcgs()
240 {
241  return Hcgs() * Ccgs();
242 }
243 
244 ////////////////////////////////////////////////////////////////////////////////
245 /// Boltzmann's constant in \f$ J K^{-1} \f$
246 /// \f[ k \f]
247 constexpr Double_t K()
248 {
249  return 1.3806503e-23;
250 }
251 
252 ////////////////////////////////////////////////////////////////////////////////
253 /// \f$ erg K^{-1} \f$
254 constexpr Double_t Kcgs()
255 {
256  return 1.0e7 * K();
257 }
258 
259 ////////////////////////////////////////////////////////////////////////////////
260 /// Boltzmann's constant uncertainty.
262 {
263  return 0.0000024e-23;
264 }
265 
266 ////////////////////////////////////////////////////////////////////////////////
267 /// Stefan-Boltzmann constant in \f$ W m^{-2} K^{-4}\f$
268 /// \f[ \sigma \f]
269 constexpr Double_t Sigma()
270 {
271  return 5.6704e-8;
272 }
273 
274 ////////////////////////////////////////////////////////////////////////////////
275 /// Stefan-Boltzmann constant uncertainty.
277 {
278  return 0.000040e-8;
279 }
280 
281 ////////////////////////////////////////////////////////////////////////////////
282 /// Avogadro constant (Avogadro's Number) in \f$ mol^{-1} \f$
283 constexpr Double_t Na()
284 {
285  return 6.02214199e+23;
286 }
287 
288 ////////////////////////////////////////////////////////////////////////////////
289 /// Avogadro constant (Avogadro's Number) uncertainty.
291 {
292  return 0.00000047e+23;
293 }
294 
295 ////////////////////////////////////////////////////////////////////////////////
296 /// [Universal gas constant](http://scienceworld.wolfram.com/physics/UniversalGasConstant.html)
297 /// (\f$ Na K \f$) in \f$ J K^{-1} mol^{-1} \f$
298 //
299 constexpr Double_t R()
300 {
301  return K() * Na();
302 }
303 
304 ////////////////////////////////////////////////////////////////////////////////
305 /// Universal gas constant uncertainty.
307 {
308  return R() * ((KUncertainty() / K()) + (NaUncertainty() / Na()));
309 }
310 
311 ////////////////////////////////////////////////////////////////////////////////
312 /// [Molecular weight of dry air 1976 US Standard Atmosphere](http://atmos.nmsu.edu/jsdap/encyclopediawork.html)
313 /// in \f$ kg kmol^{-1} \f$ or \f$ gm mol^{-1} \f$
314 constexpr Double_t MWair()
315 {
316  return 28.9644;
317 }
318 
319 ////////////////////////////////////////////////////////////////////////////////
320 /// [Dry Air Gas Constant (R / MWair)](http://atmos.nmsu.edu/education_and_outreach/encyclopedia/gas_constant.htm)
321 /// in \f$ J kg^{-1} K^{-1} \f$
322 constexpr Double_t Rgair()
323 {
324  return (1000.0 * R()) / MWair();
325 }
326 
327 ////////////////////////////////////////////////////////////////////////////////
328 /// Euler-Mascheroni Constant.
329 constexpr Double_t EulerGamma()
330 {
331  return 0.577215664901532860606512090082402431042;
332 }
333 
334 ////////////////////////////////////////////////////////////////////////////////
335 /// Elementary charge in \f$ C \f$ .
336 constexpr Double_t Qe()
337 {
338  return 1.602176462e-19;
339 }
340 
341 ////////////////////////////////////////////////////////////////////////////////
342 /// Elementary charge uncertainty.
344 {
345  return 0.000000063e-19;
346 }
347 
348 ////////////////////////////////////////////////////////////////////////////////
349 // Mathematical Functions
350 
351 ////////////////////////////////////////////////////////////////////////////////
352 // Trigonometrical Functions
353 
354 inline Double_t Sin(Double_t);
355 inline Double_t Cos(Double_t);
356 inline Double_t Tan(Double_t);
357 inline Double_t SinH(Double_t);
358 inline Double_t CosH(Double_t);
359 inline Double_t TanH(Double_t);
360 inline Double_t ASin(Double_t);
361 inline Double_t ACos(Double_t);
362 inline Double_t ATan(Double_t);
368 
369 ////////////////////////////////////////////////////////////////////////////////
370 // Elementary Functions
371 
372 inline Double_t Ceil(Double_t x);
373 inline Int_t CeilNint(Double_t x);
374 inline Double_t Floor(Double_t x);
375 inline Int_t FloorNint(Double_t x);
376 template <typename T>
377 inline Int_t Nint(T x);
378 
379 inline Double_t Sq(Double_t x);
380 inline Double_t Sqrt(Double_t x);
381 inline Double_t Exp(Double_t x);
382 inline Double_t Ldexp(Double_t x, Int_t exp);
386 inline LongDouble_t Power(Long64_t x, Long64_t y);
387 inline Double_t Power(Double_t x, Double_t y);
388 inline Double_t Power(Double_t x, Int_t y);
389 inline Double_t Log(Double_t x);
391 inline Double_t Log10(Double_t x);
392 inline Int_t Finite(Double_t x);
393 inline Int_t Finite(Float_t x);
394 inline Bool_t IsNaN(Double_t x);
395 inline Bool_t IsNaN(Float_t x);
396 
397 inline Double_t QuietNaN();
398 inline Double_t SignalingNaN();
399 inline Double_t Infinity();
400 
401 template <typename T>
402 struct Limits {
403  inline static T Min();
404  inline static T Max();
405  inline static T Epsilon();
406  };
407 
408  // Some integer math
409  Long_t Hypot(Long_t x, Long_t y); // sqrt(px*px + py*py)
410 
411  // Comparing floating points
413  //return kTRUE if absolute difference between af and bf is less than epsilon
414  return TMath::Abs(af-bf) < epsilon;
415  }
416  inline Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec) {
417  //return kTRUE if relative difference between af and bf is less than relPrec
418  return TMath::Abs(af - bf) <= 0.5 * relPrec * (TMath::Abs(af) + TMath::Abs(bf)) ||
419  TMath::Abs(af - bf) < Limits<Double_t>::Min(); // handle denormals
420  }
421 
422  /////////////////////////////////////////////////////////////////////////////
423  // Array Algorithms
424 
425  // Min, Max of an array
426  template <typename T> T MinElement(Long64_t n, const T *a);
427  template <typename T> T MaxElement(Long64_t n, const T *a);
428 
429  // Locate Min, Max element number in an array
430  template <typename T> Long64_t LocMin(Long64_t n, const T *a);
431  template <typename Iterator> Iterator LocMin(Iterator first, Iterator last);
432  template <typename T> Long64_t LocMax(Long64_t n, const T *a);
433  template <typename Iterator> Iterator LocMax(Iterator first, Iterator last);
434 
435  // Binary search
436  template <typename T> Long64_t BinarySearch(Long64_t n, const T *array, T value);
437  template <typename T> Long64_t BinarySearch(Long64_t n, const T **array, T value);
438  template <typename Iterator, typename Element> Iterator BinarySearch(Iterator first, Iterator last, Element value);
439 
440  // Hashing
441  ULong_t Hash(const void *txt, Int_t ntxt);
442  ULong_t Hash(const char *str);
443 
444  // Sorting
445  template <typename Element, typename Index>
446  void Sort(Index n, const Element* a, Index* index, Bool_t down=kTRUE);
447  template <typename Iterator, typename IndexIterator>
448  void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE);
449 
450  void BubbleHigh(Int_t Narr, Double_t *arr1, Int_t *arr2);
451  void BubbleLow (Int_t Narr, Double_t *arr1, Int_t *arr2);
452 
453  Bool_t Permute(Int_t n, Int_t *a); // Find permutations
454 
455  /////////////////////////////////////////////////////////////////////////////
456  // Geometrical Functions
457 
458  //Sample quantiles
459  void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob,
460  Bool_t isSorted=kTRUE, Int_t *index = 0, Int_t type=7);
461 
462  // IsInside
463  template <typename T> Bool_t IsInside(T xp, T yp, Int_t np, T *x, T *y);
464 
465  // Calculate the Cross Product of two vectors
466  template <typename T> T *Cross(const T v1[3],const T v2[3], T out[3]);
467 
468  Float_t Normalize(Float_t v[3]); // Normalize a vector
469  Double_t Normalize(Double_t v[3]); // Normalize a vector
470 
471  //Calculate the Normalized Cross Product of two vectors
472  template <typename T> inline T NormCross(const T v1[3],const T v2[3],T out[3]);
473 
474  // Calculate a normal vector of a plane
475  template <typename T> T *Normal2Plane(const T v1[3],const T v2[3],const T v3[3], T normal[3]);
476 
477  /////////////////////////////////////////////////////////////////////////////
478  // Polynomial Functions
479 
480  Bool_t RootsCubic(const Double_t coef[4],Double_t &a, Double_t &b, Double_t &c);
481 
482  /////////////////////////////////////////////////////////////////////////////
483  // Statistic Functions
484 
485  Double_t Binomial(Int_t n,Int_t k); // Calculate the binomial coefficient n over k
494  Double_t KolmogorovTest(Int_t na, const Double_t *a, Int_t nb, const Double_t *b, Option_t *option);
503  Double_t Prob(Double_t chi2,Int_t ndf);
507  Double_t Vavilov(Double_t x, Double_t kappa, Double_t beta2);
508  Double_t VavilovI(Double_t x, Double_t kappa, Double_t beta2);
509  Double_t Voigt(Double_t x, Double_t sigma, Double_t lg, Int_t r = 4);
510 
511  /////////////////////////////////////////////////////////////////////////////
512  // Statistics over arrays
513 
514  //Mean, Geometric Mean, Median, RMS(sigma)
515 
516  template <typename T> Double_t Mean(Long64_t n, const T *a, const Double_t *w=0);
517  template <typename Iterator> Double_t Mean(Iterator first, Iterator last);
518  template <typename Iterator, typename WeightIterator> Double_t Mean(Iterator first, Iterator last, WeightIterator wfirst);
519 
520  template <typename T> Double_t GeomMean(Long64_t n, const T *a);
521  template <typename Iterator> Double_t GeomMean(Iterator first, Iterator last);
522 
523  template <typename T> Double_t RMS(Long64_t n, const T *a, const Double_t *w=0);
524  template <typename Iterator> Double_t RMS(Iterator first, Iterator last);
525  template <typename Iterator, typename WeightIterator> Double_t RMS(Iterator first, Iterator last, WeightIterator wfirst);
526 
527  template <typename T> Double_t StdDev(Long64_t n, const T *a, const Double_t * w = 0) { return RMS<T>(n,a,w); }
528  template <typename Iterator> Double_t StdDev(Iterator first, Iterator last) { return RMS<Iterator>(first,last); }
529  template <typename Iterator, typename WeightIterator> Double_t StdDev(Iterator first, Iterator last, WeightIterator wfirst) { return RMS<Iterator,WeightIterator>(first,last,wfirst); }
530 
531  template <typename T> Double_t Median(Long64_t n, const T *a, const Double_t *w=0, Long64_t *work=0);
532 
533  //k-th order statistic
534  template <class Element, typename Size> Element KOrdStat(Size n, const Element *a, Size k, Size *work = 0);
535 
536  /////////////////////////////////////////////////////////////////////////////
537  // Special Functions
538 
544 
545  // Bessel functions
546  Double_t BesselI(Int_t n,Double_t x); /// integer order modified Bessel function I_n(x)
547  Double_t BesselK(Int_t n,Double_t x); /// integer order modified Bessel function K_n(x)
548  Double_t BesselI0(Double_t x); /// modified Bessel function I_0(x)
549  Double_t BesselK0(Double_t x); /// modified Bessel function K_0(x)
550  Double_t BesselI1(Double_t x); /// modified Bessel function I_1(x)
551  Double_t BesselK1(Double_t x); /// modified Bessel function K_1(x)
552  Double_t BesselJ0(Double_t x); /// Bessel function J0(x) for any real x
553  Double_t BesselJ1(Double_t x); /// Bessel function J1(x) for any real x
554  Double_t BesselY0(Double_t x); /// Bessel function Y0(x) for positive x
555  Double_t BesselY1(Double_t x); /// Bessel function Y1(x) for positive x
556  Double_t StruveH0(Double_t x); /// Struve functions of order 0
557  Double_t StruveH1(Double_t x); /// Struve functions of order 1
558  Double_t StruveL0(Double_t x); /// Modified Struve functions of order 0
559  Double_t StruveL1(Double_t x); /// Modified Struve functions of order 1
560 
562  Double_t Erf(Double_t x);
571 }
572 
573 ////////////////////////////////////////////////////////////////////////////////
574 // Trig and other functions
575 
576 #include <float.h>
577 
578 #if defined(R__WIN32) && !defined(__CINT__)
579 # ifndef finite
580 # define finite _finite
581 # endif
582 #endif
583 #if defined(R__AIX) || defined(R__SOLARIS_CC50) || \
584  defined(R__HPUX11) || defined(R__GLIBC) || \
585  (defined(R__MACOSX) )
586 // math functions are defined inline so we have to include them here
587 # include <math.h>
588 # ifdef R__SOLARIS_CC50
589  extern "C" { int finite(double); }
590 # endif
591 // # if defined(R__GLIBC) && defined(__STRICT_ANSI__)
592 // # ifndef finite
593 // # define finite __finite
594 // # endif
595 // # ifndef isnan
596 // # define isnan __isnan
597 // # endif
598 // # endif
599 #else
600 // don't want to include complete <math.h>
601 extern "C" {
602  extern double sin(double);
603  extern double cos(double);
604  extern double tan(double);
605  extern double sinh(double);
606  extern double cosh(double);
607  extern double tanh(double);
608  extern double asin(double);
609  extern double acos(double);
610  extern double atan(double);
611  extern double atan2(double, double);
612  extern double sqrt(double);
613  extern double exp(double);
614  extern double pow(double, double);
615  extern double log(double);
616  extern double log10(double);
617 #ifndef R__WIN32
618 # if !defined(finite)
619  extern int finite(double);
620 # endif
621 # if !defined(isnan)
622  extern int isnan(double);
623 # endif
624  extern double ldexp(double, int);
625  extern double ceil(double);
626  extern double floor(double);
627 #else
628  _CRTIMP double ldexp(double, int);
629  _CRTIMP double ceil(double);
630  _CRTIMP double floor(double);
631 #endif
632 }
633 #endif
634 
635 ////////////////////////////////////////////////////////////////////////////////
637  { return sin(x); }
638 
639 ////////////////////////////////////////////////////////////////////////////////
641  { return cos(x); }
642 
643 ////////////////////////////////////////////////////////////////////////////////
645  { return tan(x); }
646 
647 ////////////////////////////////////////////////////////////////////////////////
649  { return sinh(x); }
650 
651 ////////////////////////////////////////////////////////////////////////////////
653  { return cosh(x); }
654 
655 ////////////////////////////////////////////////////////////////////////////////
657  { return tanh(x); }
658 
659 ////////////////////////////////////////////////////////////////////////////////
661  { if (x < -1.) return -TMath::Pi()/2;
662  if (x > 1.) return TMath::Pi()/2;
663  return asin(x);
664  }
665 
666 ////////////////////////////////////////////////////////////////////////////////
668  { if (x < -1.) return TMath::Pi();
669  if (x > 1.) return 0;
670  return acos(x);
671  }
672 
673 ////////////////////////////////////////////////////////////////////////////////
675  { return atan(x); }
676 
677 ////////////////////////////////////////////////////////////////////////////////
679  { if (x != 0) return atan2(y, x);
680  if (y == 0) return 0;
681  if (y > 0) return Pi()/2;
682  else return -Pi()/2;
683  }
684 
685 ////////////////////////////////////////////////////////////////////////////////
687  { return x*x; }
688 
689 ////////////////////////////////////////////////////////////////////////////////
691  { return sqrt(x); }
692 
693 ////////////////////////////////////////////////////////////////////////////////
695  { return ceil(x); }
696 
697 ////////////////////////////////////////////////////////////////////////////////
699  { return TMath::Nint(ceil(x)); }
700 
701 ////////////////////////////////////////////////////////////////////////////////
703  { return floor(x); }
704 
705 ////////////////////////////////////////////////////////////////////////////////
707  { return TMath::Nint(floor(x)); }
708 
709 ////////////////////////////////////////////////////////////////////////////////
710 /// Round to nearest integer. Rounds half integers to the nearest even integer.
711 template<typename T>
713 {
714  int i;
715  if (x >= 0) {
716  i = int(x + 0.5);
717  if ( i & 1 && x + 0.5 == T(i) ) i--;
718  } else {
719  i = int(x - 0.5);
720  if ( i & 1 && x - 0.5 == T(i) ) i++;
721  }
722  return i;
723 }
724 
725 ////////////////////////////////////////////////////////////////////////////////
727  { return exp(x); }
728 
729 ////////////////////////////////////////////////////////////////////////////////
731  { return ldexp(x, exp); }
732 
733 ////////////////////////////////////////////////////////////////////////////////
735  { return std::pow(x,y); }
736 
737 ////////////////////////////////////////////////////////////////////////////////
739  { return std::pow(x,(LongDouble_t)y); }
740 
741 ////////////////////////////////////////////////////////////////////////////////
743  { return std::pow(x,y); }
744 
745 ////////////////////////////////////////////////////////////////////////////////
747  { return pow(x, y); }
748 
749 ////////////////////////////////////////////////////////////////////////////////
751 #ifdef R__ANSISTREAM
752  return std::pow(x, y);
753 #else
754  return pow(x, (Double_t) y);
755 #endif
756 }
757 
758 ////////////////////////////////////////////////////////////////////////////////
760  { return log(x); }
761 
762 ////////////////////////////////////////////////////////////////////////////////
764  { return log10(x); }
765 
766 ////////////////////////////////////////////////////////////////////////////////
767 /// Check if it is finite with a mask in order to be consistent in presence of
768 /// fast math.
769 /// Inspired from the CMSSW FWCore/Utilities package
771 #if defined(R__FAST_MATH)
772 
773 {
774  const unsigned long long mask = 0x7FF0000000000000LL;
775  union { unsigned long long l; double d;} v;
776  v.d =x;
777  return (v.l&mask)!=mask;
778 }
779 #else
780 # if defined(R__HPUX11)
781  { return isfinite(x); }
782 # elif defined(R__MACOSX)
783 # ifdef isfinite
784  // from math.h
785  { return isfinite(x); }
786 # else
787  // from cmath
788  { return std::isfinite(x); }
789 # endif
790 # else
791  { return finite(x); }
792 # endif
793 #endif
794 
795 ////////////////////////////////////////////////////////////////////////////////
796 /// Check if it is finite with a mask in order to be consistent in presence of
797 /// fast math.
798 /// Inspired from the CMSSW FWCore/Utilities package
800 #if defined(R__FAST_MATH)
801 
802 {
803  const unsigned int mask = 0x7f800000;
804  union { unsigned int l; float d;} v;
805  v.d =x;
806  return (v.l&mask)!=mask;
807 }
808 #else
809 { return std::isfinite(x); }
810 #endif
811 
812 // This namespace provides all the routines necessary for checking if a number
813 // is a NaN also in presence of optimisations affecting the behaviour of the
814 // floating point calculations.
815 // Inspired from the CMSSW FWCore/Utilities package
816 
817 #if defined (R__FAST_MATH)
818 namespace ROOT {
819 namespace Internal {
820 namespace Math {
821 // abridged from GNU libc 2.6.1 - in detail from
822 // math/math_private.h
823 // sysdeps/ieee754/ldbl-96/math_ldbl.h
824 
825 // part of ths file:
826  /*
827  * ====================================================
828  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
829  *
830  * Developed at SunPro, a Sun Microsystems, Inc. business.
831  * Permission to use, copy, modify, and distribute this
832  * software is freely granted, provided that this notice
833  * is preserved.
834  * ====================================================
835  */
836 
837  // A union which permits us to convert between a double and two 32 bit ints.
838  typedef union {
839  Double_t value;
840  struct {
841  UInt_t lsw;
842  UInt_t msw;
843  } parts;
844  } ieee_double_shape_type;
845 
846 #define EXTRACT_WORDS(ix0,ix1,d) \
847  do { \
848  ieee_double_shape_type ew_u; \
849  ew_u.value = (d); \
850  (ix0) = ew_u.parts.msw; \
851  (ix1) = ew_u.parts.lsw; \
852  } while (0)
853 
854  inline Bool_t IsNaN(Double_t x)
855  {
856  UInt_t hx, lx;
857 
858  EXTRACT_WORDS(hx, lx, x);
859 
860  lx |= hx & 0xfffff;
861  hx &= 0x7ff00000;
862  return (hx == 0x7ff00000) && (lx != 0);
863  }
864 
865  typedef union {
866  Float_t value;
867  UInt_t word;
868  } ieee_float_shape_type;
869 
870 #define GET_FLOAT_WORD(i,d) \
871  do { \
872  ieee_float_shape_type gf_u; \
873  gf_u.value = (d); \
874  (i) = gf_u.word; \
875  } while (0)
876 
877  inline Bool_t IsNaN(Float_t x)
878  {
879  UInt_t wx;
880  GET_FLOAT_WORD (wx, x);
881  wx &= 0x7fffffff;
882  return (Bool_t)(wx > 0x7f800000);
883  }
884 } } } // end NS ROOT::Internal::Math
885 #endif // End R__FAST_MATH
886 
887 #if defined(R__FAST_MATH)
889  inline Bool_t TMath::IsNaN(Float_t x) { return ROOT::Internal::Math::IsNaN(x); }
890 #else
891  inline Bool_t TMath::IsNaN(Double_t x) { return std::isnan(x); }
892  inline Bool_t TMath::IsNaN(Float_t x) { return std::isnan(x); }
893 #endif
894 
895 ////////////////////////////////////////////////////////////////////////////////
896 // Wrapper to numeric_limits
897 
898 ////////////////////////////////////////////////////////////////////////////////
899 /// Returns a quiet NaN as [defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Quiet_NaN)
901 
902  return std::numeric_limits<Double_t>::quiet_NaN();
903 }
904 
905 ////////////////////////////////////////////////////////////////////////////////
906 /// Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN)
908  return std::numeric_limits<Double_t>::signaling_NaN();
909 }
910 
911 ////////////////////////////////////////////////////////////////////////////////
912 /// Returns an infinity as defined by the IEEE standard
914  return std::numeric_limits<Double_t>::infinity();
915 }
916 
917 ////////////////////////////////////////////////////////////////////////////////
918 /// Returns maximum representation for type T
919 template<typename T>
921  return (std::numeric_limits<T>::min)(); //N.B. use this signature to avoid class with macro min() on Windows
922 }
923 
924 ////////////////////////////////////////////////////////////////////////////////
925 /// Returns minimum double representation
926 template<typename T>
928  return (std::numeric_limits<T>::max)(); //N.B. use this signature to avoid class with macro max() on Windows
929 }
930 
931 ////////////////////////////////////////////////////////////////////////////////
932 /// Returns minimum double representation
933 template<typename T>
936 }
937 
938 ////////////////////////////////////////////////////////////////////////////////
939 // Advanced.
940 
941 ////////////////////////////////////////////////////////////////////////////////
942 /// Calculate the Normalized Cross Product of two vectors
943 template <typename T> inline T TMath::NormCross(const T v1[3],const T v2[3],T out[3])
944 {
945  return Normalize(Cross(v1,v2,out));
946 }
947 
948 ////////////////////////////////////////////////////////////////////////////////
949 /// Return minimum of array a of length n.
950 template <typename T>
952  return *std::min_element(a,a+n);
953 }
954 
955 ////////////////////////////////////////////////////////////////////////////////
956 /// Return maximum of array a of length n.
957 template <typename T>
959  return *std::max_element(a,a+n);
960 }
961 
962 ////////////////////////////////////////////////////////////////////////////////
963 /// Return index of array with the minimum element.
964 /// If more than one element is minimum returns first found.
965 ///
966 /// Implement here since this one is found to be faster (mainly on 64 bit machines)
967 /// than stl generic implementation.
968 /// When performing the comparison, the STL implementation needs to de-reference both the array iterator
969 /// and the iterator pointing to the resulting minimum location
970 template <typename T>
972  if (n <= 0 || !a) return -1;
973  T xmin = a[0];
974  Long64_t loc = 0;
975  for (Long64_t i = 1; i < n; i++) {
976  if (xmin > a[i]) {
977  xmin = a[i];
978  loc = i;
979  }
980  }
981  return loc;
982 }
983 
984 ////////////////////////////////////////////////////////////////////////////////
985 /// Return index of array with the minimum element.
986 /// If more than one element is minimum returns first found.
987 template <typename Iterator>
988 Iterator TMath::LocMin(Iterator first, Iterator last) {
989 
990  return std::min_element(first, last);
991 }
992 
993 ////////////////////////////////////////////////////////////////////////////////
994 /// Return index of array with the maximum element.
995 /// If more than one element is maximum returns first found.
996 ///
997 /// Implement here since it is faster (see comment in LocMin function)
998 template <typename T>
1000  if (n <= 0 || !a) return -1;
1001  T xmax = a[0];
1002  Long64_t loc = 0;
1003  for (Long64_t i = 1; i < n; i++) {
1004  if (xmax < a[i]) {
1005  xmax = a[i];
1006  loc = i;
1007  }
1008  }
1009  return loc;
1010 }
1011 
1012 ////////////////////////////////////////////////////////////////////////////////
1013 /// Return index of array with the maximum element.
1014 /// If more than one element is maximum returns first found.
1015 template <typename Iterator>
1016 Iterator TMath::LocMax(Iterator first, Iterator last)
1017 {
1018 
1019  return std::max_element(first, last);
1020 }
1021 
1022 template<typename T>
1023 struct CompareDesc {
1024 
1025  CompareDesc(T d) : fData(d) {}
1026 
1027  template<typename Index>
1028  bool operator()(Index i1, Index i2) {
1029  return *(fData + i1) > *(fData + i2);
1030  }
1031 
1033 };
1034 
1035 template<typename T>
1036 struct CompareAsc {
1037 
1038  CompareAsc(T d) : fData(d) {}
1039 
1040  template<typename Index>
1041  bool operator()(Index i1, Index i2) {
1042  return *(fData + i1) < *(fData + i2);
1043  }
1044 
1046 };
1047 
1048 ////////////////////////////////////////////////////////////////////////////////
1049 /// Return the weighted mean of an array defined by the iterators.
1050 template <typename Iterator>
1051 Double_t TMath::Mean(Iterator first, Iterator last)
1052 {
1053  Double_t sum = 0;
1054  Double_t sumw = 0;
1055  while ( first != last )
1056  {
1057  sum += *first;
1058  sumw += 1;
1059  first++;
1060  }
1061 
1062  return sum/sumw;
1063 }
1064 
1065 ////////////////////////////////////////////////////////////////////////////////
1066 /// Return the weighted mean of an array defined by the first and
1067 /// last iterators. The w iterator should point to the first element
1068 /// of a vector of weights of the same size as the main array.
1069 template <typename Iterator, typename WeightIterator>
1070 Double_t TMath::Mean(Iterator first, Iterator last, WeightIterator w)
1071 {
1072 
1073  Double_t sum = 0;
1074  Double_t sumw = 0;
1075  int i = 0;
1076  while ( first != last ) {
1077  if ( *w < 0) {
1078  ::Error("TMath::Mean","w[%d] = %.4e < 0 ?!",i,*w);
1079  return 0;
1080  }
1081  sum += (*w) * (*first);
1082  sumw += (*w) ;
1083  ++w;
1084  ++first;
1085  ++i;
1086  }
1087  if (sumw <= 0) {
1088  ::Error("TMath::Mean","sum of weights == 0 ?!");
1089  return 0;
1090  }
1091 
1092  return sum/sumw;
1093 }
1094 
1095 ////////////////////////////////////////////////////////////////////////////////
1096 /// Return the weighted mean of an array a with length n.
1097 template <typename T>
1099 {
1100  if (w) {
1101  return TMath::Mean(a, a+n, w);
1102  } else {
1103  return TMath::Mean(a, a+n);
1104  }
1105 }
1106 
1107 ////////////////////////////////////////////////////////////////////////////////
1108 /// Return the geometric mean of an array defined by the iterators.
1109 /// \f[ GeomMean = (\prod_{i=0}^{n-1} |a[i]|)^{1/n} \f]
1110 template <typename Iterator>
1111 Double_t TMath::GeomMean(Iterator first, Iterator last)
1112 {
1113  Double_t logsum = 0.;
1114  Long64_t n = 0;
1115  while ( first != last ) {
1116  if (*first == 0) return 0.;
1117  Double_t absa = (Double_t) TMath::Abs(*first);
1118  logsum += TMath::Log(absa);
1119  ++first;
1120  ++n;
1121  }
1122 
1123  return TMath::Exp(logsum/n);
1124 }
1125 
1126 ////////////////////////////////////////////////////////////////////////////////
1127 /// Return the geometric mean of an array a of size n.
1128 /// \f[ GeomMean = (\prod_{i=0}^{n-1} |a[i]|)^{1/n} \f]
1129 template <typename T>
1131 {
1132  return TMath::GeomMean(a, a+n);
1133 }
1134 
1135 ////////////////////////////////////////////////////////////////////////////////
1136 /// Return the Standard Deviation of an array defined by the iterators.
1137 /// Note that this function returns the sigma(standard deviation) and
1138 /// not the root mean square of the array.
1139 ///
1140 /// Use the two pass algorithm, which is slower (! a factor of 2) but much more
1141 /// precise. Since we have a vector the 2 pass algorithm is still faster than the
1142 /// Welford algorithm. (See also ROOT-5545)
1143 template <typename Iterator>
1144 Double_t TMath::RMS(Iterator first, Iterator last)
1145 {
1146 
1147  Double_t n = 0;
1148 
1149  Double_t tot = 0;
1150  Double_t mean = TMath::Mean(first,last);
1151  while ( first != last ) {
1152  Double_t x = Double_t(*first);
1153  tot += (x - mean)*(x - mean);
1154  ++first;
1155  ++n;
1156  }
1157  Double_t rms = (n > 1) ? TMath::Sqrt(tot/(n-1)) : 0.0;
1158  return rms;
1159 }
1160 
1161 ////////////////////////////////////////////////////////////////////////////////
1162 /// Return the weighted Standard Deviation of an array defined by the iterators.
1163 /// Note that this function returns the sigma(standard deviation) and
1164 /// not the root mean square of the array.
1165 ///
1166 /// As in the unweighted case use the two pass algorithm
1167 template <typename Iterator, typename WeightIterator>
1168 Double_t TMath::RMS(Iterator first, Iterator last, WeightIterator w)
1169 {
1170  Double_t tot = 0;
1171  Double_t sumw = 0;
1172  Double_t sumw2 = 0;
1173  Double_t mean = TMath::Mean(first,last,w);
1174  while ( first != last ) {
1175  Double_t x = Double_t(*first);
1176  sumw += *w;
1177  sumw2 += (*w) * (*w);
1178  tot += (*w) * (x - mean)*(x - mean);
1179  ++first;
1180  ++w;
1181  }
1182  // use the correction neff/(neff -1) for the unbiased formula
1183  Double_t rms = TMath::Sqrt(tot * sumw/ (sumw*sumw - sumw2) );
1184  return rms;
1185 }
1186 
1187 ////////////////////////////////////////////////////////////////////////////////
1188 /// Return the Standard Deviation of an array a with length n.
1189 /// Note that this function returns the sigma(standard deviation) and
1190 /// not the root mean square of the array.
1191 template <typename T>
1192 Double_t TMath::RMS(Long64_t n, const T *a, const Double_t * w)
1193 {
1194  return (w) ? TMath::RMS(a, a+n, w) : TMath::RMS(a, a+n);
1195 }
1196 
1197 ////////////////////////////////////////////////////////////////////////////////
1198 /// Binary search in an array defined by its iterators.
1199 ///
1200 /// The values in the iterators range are supposed to be sorted
1201 /// prior to this call. If match is found, function returns
1202 /// position of element. If no match found, function gives nearest
1203 /// element smaller than value.
1204 template <typename Iterator, typename Element>
1205 Iterator TMath::BinarySearch(Iterator first, Iterator last, Element value)
1206 {
1207  Iterator pind;
1208  pind = std::lower_bound(first, last, value);
1209  if ( (pind != last) && (*pind == value) )
1210  return pind;
1211  else
1212  return ( pind - 1);
1213 }
1214 
1215 ////////////////////////////////////////////////////////////////////////////////
1216 /// Binary search in an array of n values to locate value.
1217 ///
1218 /// Array is supposed to be sorted prior to this call.
1219 /// If match is found, function returns position of element.
1220 /// If no match found, function gives nearest element smaller than value.
1221 template <typename T> Long64_t TMath::BinarySearch(Long64_t n, const T *array, T value)
1222 {
1223  const T* pind;
1224  pind = std::lower_bound(array, array + n, value);
1225  if ( (pind != array + n) && (*pind == value) )
1226  return (pind - array);
1227  else
1228  return ( pind - array - 1);
1229 }
1230 
1231 ////////////////////////////////////////////////////////////////////////////////
1232 /// Binary search in an array of n values to locate value.
1233 ///
1234 /// Array is supposed to be sorted prior to this call.
1235 /// If match is found, function returns position of element.
1236 /// If no match found, function gives nearest element smaller than value.
1237 template <typename T> Long64_t TMath::BinarySearch(Long64_t n, const T **array, T value)
1238 {
1239  const T* pind;
1240  pind = std::lower_bound(*array, *array + n, value);
1241  if ( (pind != *array + n) && (*pind == value) )
1242  return (pind - *array);
1243  else
1244  return ( pind - *array - 1);
1245 }
1246 
1247 ////////////////////////////////////////////////////////////////////////////////
1248 /// Sort the n1 elements of the Short_t array defined by its
1249 /// iterators. In output the array index contains the indices of
1250 /// the sorted array. If down is false sort in increasing order
1251 /// (default is decreasing order).
1252 ///
1253 /// NOTE that the array index must be created with a length bigger
1254 /// or equal than the main array before calling this function.
1255 template <typename Iterator, typename IndexIterator>
1256 void TMath::SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down)
1257 {
1258  int i = 0;
1259 
1260  IndexIterator cindex = index;
1261  for ( Iterator cfirst = first; cfirst != last; ++cfirst )
1262  {
1263  *cindex = i++;
1264  ++cindex;
1265  }
1266 
1267  if ( down )
1268  std::sort(index, cindex, CompareDesc<Iterator>(first) );
1269  else
1270  std::sort(index, cindex, CompareAsc<Iterator>(first) );
1271 }
1272 
1273 ////////////////////////////////////////////////////////////////////////////////
1274 /// Sort the n elements of the array a of generic templated type Element.
1275 /// In output the array index of type Index contains the indices of the sorted array.
1276 /// If down is false sort in increasing order (default is decreasing order).
1277 ///
1278 /// NOTE that the array index must be created with a length >= n
1279 /// before calling this function.
1280 /// NOTE also that the size type for n must be the same type used for the index array
1281 /// (templated type Index)
1282 template <typename Element, typename Index> void TMath::Sort(Index n, const Element* a, Index* index, Bool_t down)
1283 {
1284 
1285  for(Index i = 0; i < n; i++) { index[i] = i; }
1286  if ( down )
1287  std::sort(index, index + n, CompareDesc<const Element*>(a) );
1288  else
1289  std::sort(index, index + n, CompareAsc<const Element*>(a) );
1290 }
1291 
1292 ////////////////////////////////////////////////////////////////////////////////
1293 /// Calculate the Cross Product of two vectors:
1294 /// out = [v1 x v2]
1295 template <typename T> T *TMath::Cross(const T v1[3],const T v2[3], T out[3])
1296 {
1297  out[0] = v1[1] * v2[2] - v1[2] * v2[1];
1298  out[1] = v1[2] * v2[0] - v1[0] * v2[2];
1299  out[2] = v1[0] * v2[1] - v1[1] * v2[0];
1300 
1301  return out;
1302 }
1303 
1304 ////////////////////////////////////////////////////////////////////////////////
1305 /// Calculate a normal vector of a plane.
1306 ///
1307 /// \param[in] p1, p2,p3 3 3D points belonged the plane to define it.
1308 /// \param[out] normal Pointer to 3D normal vector (normalized)
1309 template <typename T> T * TMath::Normal2Plane(const T p1[3],const T p2[3],const T p3[3], T normal[3])
1310 {
1311  T v1[3], v2[3];
1312 
1313  v1[0] = p2[0] - p1[0];
1314  v1[1] = p2[1] - p1[1];
1315  v1[2] = p2[2] - p1[2];
1316 
1317  v2[0] = p3[0] - p1[0];
1318  v2[1] = p3[1] - p1[1];
1319  v2[2] = p3[2] - p1[2];
1320 
1321  NormCross(v1,v2,normal);
1322  return normal;
1323 }
1324 
1325 ////////////////////////////////////////////////////////////////////////////////
1326 /// Function which returns kTRUE if point xp,yp lies inside the
1327 /// polygon defined by the np points in arrays x and y, kFALSE otherwise.
1328 /// Note that the polygon may be open or closed.
1329 template <typename T> Bool_t TMath::IsInside(T xp, T yp, Int_t np, T *x, T *y)
1330 {
1331  Int_t i, j = np-1 ;
1332  Bool_t oddNodes = kFALSE;
1333 
1334  for (i=0; i<np; i++) {
1335  if ((y[i]<yp && y[j]>=yp) || (y[j]<yp && y[i]>=yp)) {
1336  if (x[i]+(yp-y[i])/(y[j]-y[i])*(x[j]-x[i])<xp) {
1337  oddNodes = !oddNodes;
1338  }
1339  }
1340  j=i;
1341  }
1342 
1343  return oddNodes;
1344 }
1345 
1346 ////////////////////////////////////////////////////////////////////////////////
1347 /// Return the median of the array a where each entry i has weight w[i] .
1348 /// Both arrays have a length of at least n . The median is a number obtained
1349 /// from the sorted array a through
1350 ///
1351 /// median = (a[jl]+a[jh])/2. where (using also the sorted index on the array w)
1352 ///
1353 /// sum_i=0,jl w[i] <= sumTot/2
1354 /// sum_i=0,jh w[i] >= sumTot/2
1355 /// sumTot = sum_i=0,n w[i]
1356 ///
1357 /// If w=0, the algorithm defaults to the median definition where it is
1358 /// a number that divides the sorted sequence into 2 halves.
1359 /// When n is odd or n > 1000, the median is kth element k = (n + 1) / 2.
1360 /// when n is even and n < 1000the median is a mean of the elements k = n/2 and k = n/2 + 1.
1361 ///
1362 /// If the weights are supplied (w not 0) all weights must be >= 0
1363 ///
1364 /// If work is supplied, it is used to store the sorting index and assumed to be
1365 /// >= n . If work=0, local storage is used, either on the stack if n < kWorkMax
1366 /// or on the heap for n >= kWorkMax .
1367 template <typename T> Double_t TMath::Median(Long64_t n, const T *a, const Double_t *w, Long64_t *work)
1368 {
1369 
1370  const Int_t kWorkMax = 100;
1371 
1372  if (n <= 0 || !a) return 0;
1373  Bool_t isAllocated = kFALSE;
1374  Double_t median;
1375  Long64_t *ind;
1376  Long64_t workLocal[kWorkMax];
1377 
1378  if (work) {
1379  ind = work;
1380  } else {
1381  ind = workLocal;
1382  if (n > kWorkMax) {
1383  isAllocated = kTRUE;
1384  ind = new Long64_t[n];
1385  }
1386  }
1387 
1388  if (w) {
1389  Double_t sumTot2 = 0;
1390  for (Int_t j = 0; j < n; j++) {
1391  if (w[j] < 0) {
1392  ::Error("TMath::Median","w[%d] = %.4e < 0 ?!",j,w[j]);
1393  if (isAllocated) delete [] ind;
1394  return 0;
1395  }
1396  sumTot2 += w[j];
1397  }
1398 
1399  sumTot2 /= 2.;
1400 
1401  Sort(n, a, ind, kFALSE);
1402 
1403  Double_t sum = 0.;
1404  Int_t jl;
1405  for (jl = 0; jl < n; jl++) {
1406  sum += w[ind[jl]];
1407  if (sum >= sumTot2) break;
1408  }
1409 
1410  Int_t jh;
1411  sum = 2.*sumTot2;
1412  for (jh = n-1; jh >= 0; jh--) {
1413  sum -= w[ind[jh]];
1414  if (sum <= sumTot2) break;
1415  }
1416 
1417  median = 0.5*(a[ind[jl]]+a[ind[jh]]);
1418 
1419  } else {
1420 
1421  if (n%2 == 1)
1422  median = KOrdStat(n, a,n/2, ind);
1423  else {
1424  median = 0.5*(KOrdStat(n, a, n/2 -1, ind)+KOrdStat(n, a, n/2, ind));
1425  }
1426  }
1427 
1428  if (isAllocated)
1429  delete [] ind;
1430  return median;
1431 }
1432 
1433 ////////////////////////////////////////////////////////////////////////////////
1434 /// Returns k_th order statistic of the array a of size n
1435 /// (k_th smallest element out of n elements).
1436 ///
1437 /// C-convention is used for array indexing, so if you want
1438 /// the second smallest element, call KOrdStat(n, a, 1).
1439 ///
1440 /// If work is supplied, it is used to store the sorting index and
1441 /// assumed to be >= n. If work=0, local storage is used, either on
1442 /// the stack if n < kWorkMax or on the heap for n >= kWorkMax.
1443 /// Note that the work index array will not contain the sorted indices but
1444 /// all indices of the smaller element in arbitrary order in work[0,...,k-1] and
1445 /// all indices of the larger element in arbitrary order in work[k+1,..,n-1]
1446 /// work[k] will contain instead the index of the returned element.
1447 ///
1448 /// Taken from "Numerical Recipes in C++" without the index array
1449 /// implemented by Anna Khreshuk.
1450 ///
1451 /// See also the declarations at the top of this file
1452 template <class Element, typename Size>
1453 Element TMath::KOrdStat(Size n, const Element *a, Size k, Size *work)
1454 {
1455 
1456  const Int_t kWorkMax = 100;
1457 
1458  typedef Size Index;
1459 
1460  Bool_t isAllocated = kFALSE;
1461  Size i, ir, j, l, mid;
1462  Index arr;
1463  Index *ind;
1464  Index workLocal[kWorkMax];
1465  Index temp;
1466 
1467  if (work) {
1468  ind = work;
1469  } else {
1470  ind = workLocal;
1471  if (n > kWorkMax) {
1472  isAllocated = kTRUE;
1473  ind = new Index[n];
1474  }
1475  }
1476 
1477  for (Size ii=0; ii<n; ii++) {
1478  ind[ii]=ii;
1479  }
1480  Size rk = k;
1481  l=0;
1482  ir = n-1;
1483  for(;;) {
1484  if (ir<=l+1) { //active partition contains 1 or 2 elements
1485  if (ir == l+1 && a[ind[ir]]<a[ind[l]])
1486  {temp = ind[l]; ind[l]=ind[ir]; ind[ir]=temp;}
1487  Element tmp = a[ind[rk]];
1488  if (isAllocated)
1489  delete [] ind;
1490  return tmp;
1491  } else {
1492  mid = (l+ir) >> 1; //choose median of left, center and right
1493  {temp = ind[mid]; ind[mid]=ind[l+1]; ind[l+1]=temp;}//elements as partitioning element arr.
1494  if (a[ind[l]]>a[ind[ir]]) //also rearrange so that a[l]<=a[l+1]
1495  {temp = ind[l]; ind[l]=ind[ir]; ind[ir]=temp;}
1496 
1497  if (a[ind[l+1]]>a[ind[ir]])
1498  {temp=ind[l+1]; ind[l+1]=ind[ir]; ind[ir]=temp;}
1499 
1500  if (a[ind[l]]>a[ind[l+1]])
1501  {temp = ind[l]; ind[l]=ind[l+1]; ind[l+1]=temp;}
1502 
1503  i=l+1; //initialize pointers for partitioning
1504  j=ir;
1505  arr = ind[l+1];
1506  for (;;){
1507  do i++; while (a[ind[i]]<a[arr]);
1508  do j--; while (a[ind[j]]>a[arr]);
1509  if (j<i) break; //pointers crossed, partitioning complete
1510  {temp=ind[i]; ind[i]=ind[j]; ind[j]=temp;}
1511  }
1512  ind[l+1]=ind[j];
1513  ind[j]=arr;
1514  if (j>=rk) ir = j-1; //keep active the partition that
1515  if (j<=rk) l=i; //contains the k_th element
1516  }
1517  }
1518 }
1519 
1520 #endif
constexpr Double_t H()
Planck&#39;s constant in .
Definition: TMath.h:189
Double_t ACosH(Double_t)
Definition: TMath.cxx:77
Double_t BesselI(Int_t n, Double_t x)
constexpr Double_t QeUncertainty()
Elementary charge uncertainty.
Definition: TMath.h:343
constexpr Double_t HC()
in
Definition: TMath.h:232
Double_t FDist(Double_t F, Double_t N, Double_t M)
Double_t Landau(Double_t x, Double_t mpv=0, Double_t sigma=1, Bool_t norm=kFALSE)
The LANDAU function.
Definition: TMath.cxx:469
Double_t ErfInverse(Double_t x)
returns the inverse error function x must be <-1<x<1
Definition: TMath.cxx:203
static long int sum(long int i)
Definition: Factory.cxx:2258
Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1)
Calculate a Breit Wigner function with mean and gamma.
Definition: TMath.cxx:437
float xmin
Definition: THbookFile.cxx:93
constexpr Double_t K()
Boltzmann&#39;s constant in .
Definition: TMath.h:247
Double_t FDistI(Double_t F, Double_t N, Double_t M)
Double_t PoissonI(Double_t x, Double_t par)
Double_t Floor(Double_t x)
Definition: TMath.h:702
double tanh(double)
Double_t TanH(Double_t)
Definition: TMath.h:656
Double_t Sq(Double_t x)
Definition: TMath.h:686
long long Long64_t
Definition: RtypesCore.h:69
Float_t Normalize(Float_t v[3])
Normalize a vector v in place.
Definition: TMath.cxx:495
auto * m
Definition: textangle.C:8
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t mid
Definition: TRolke.cxx:630
Long64_t LocMax(Long64_t n, const T *a)
Return index of array with the maximum element.
Definition: TMath.h:999
bool operator()(Index i1, Index i2)
Definition: TMath.h:1028
static double p3(double t, double a, double b, double c, double d)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
Double_t Log(Double_t x)
Definition: TMath.h:759
float Float_t
Definition: RtypesCore.h:53
constexpr Double_t Sqrt2()
Definition: TMath.h:89
Bool_t IsNaN(Float_t x)
Definition: TMath.h:892
const char Option_t
Definition: RtypesCore.h:62
Double_t KolmogorovProb(Double_t z)
constexpr Double_t Hcgs()
Definition: TMath.h:196
constexpr Double_t TwoPi()
Definition: TMath.h:45
const char * Size
Definition: TXMLSetup.cxx:55
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition: TMath.h:900
double T(double x)
Definition: ChebyshevPol.h:34
Double_t LaplaceDistI(Double_t x, Double_t alpha=0, Double_t beta=1)
Double_t StudentQuantile(Double_t p, Double_t ndf, Bool_t lower_tail=kTRUE)
Double_t NormQuantile(Double_t p)
Double_t BetaDist(Double_t x, Double_t p, Double_t q)
Double_t LaplaceDist(Double_t x, Double_t alpha=0, Double_t beta=1)
constexpr Double_t Rgair()
Dry Air Gas Constant (R / MWair) in
Definition: TMath.h:322
Double_t Log2(Double_t x)
Definition: TMath.cxx:101
#define N
Double_t StruveH1(Double_t x)
Struve functions of order 0.
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
CompareAsc(T d)
Definition: TMath.h:1038
constexpr Double_t SigmaUncertainty()
Stefan-Boltzmann constant uncertainty.
Definition: TMath.h:276
Bool_t IsNaN(Double_t x)
Definition: TMath.h:891
void BubbleHigh(Int_t Narr, Double_t *arr1, Int_t *arr2)
Int_t FloorNint(Double_t x)
Definition: TMath.h:706
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition: TMath.cxx:348
constexpr Double_t InvPi()
Definition: TMath.h:66
Double_t GeomMean(Long64_t n, const T *a)
Return the geometric mean of an array a of size n.
Definition: TMath.h:1130
Double_t BesselJ1(Double_t x)
Bessel function J0(x) for any real x.
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
double cos(double)
Double_t Prob(Double_t chi2, Int_t ndf)
Double_t StudentI(Double_t T, Double_t ndf)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:734
double beta(double x, double y)
Calculates the beta function.
Double_t Ldexp(Double_t x, Int_t exp)
Definition: TMath.h:730
T * Normal2Plane(const T v1[3], const T v2[3], const T v3[3], T normal[3])
Calculate a normal vector of a plane.
Definition: TMath.h:1309
constexpr Double_t Gn()
Standard acceleration of gravity in .
Definition: TMath.h:174
constexpr Double_t Ln10()
Natural log of 10 (to convert log to ln)
Definition: TMath.h:104
constexpr Double_t Hbarcgs()
Definition: TMath.h:218
Double_t RMS(Long64_t n, const T *a, const Double_t *w=0)
Return the Standard Deviation of an array a with length n.
Definition: TMath.h:1192
T NormCross(const T v1[3], const T v2[3], T out[3])
Calculate the Normalized Cross Product of two vectors.
Definition: TMath.h:943
double sqrt(double)
Double_t BesselY1(Double_t x)
Bessel function Y0(x) for positive x.
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
Double_t x[n]
Definition: legend1.C:17
double acos(double)
double sinh(double)
void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted=kTRUE, Int_t *index=0, Int_t type=7)
double log10(double)
Double_t ASinH(Double_t)
Definition: TMath.cxx:64
Double_t Log10(Double_t x)
Definition: TMath.h:763
Double_t BesselI1(Double_t x)
modified Bessel function K_0(x)
static double p2(double t, double a, double b, double c)
Double_t Freq(Double_t x)
Computation of the normal frequency function freq(x).
Definition: TMath.cxx:265
constexpr Double_t GhbarCUncertainty()
uncertainty.
Definition: TMath.h:167
constexpr Double_t HCcgs()
Definition: TMath.h:239
double pow(double, double)
bool operator()(Index i1, Index i2)
Definition: TMath.h:1041
constexpr Double_t DegToRad()
Conversion from degree to radian: .
Definition: TMath.h:82
constexpr Double_t MWair()
Molecular weight of dry air 1976 US Standard Atmosphere in or
Definition: TMath.h:314
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Sort the n elements of the array a of generic templated type Element.
Definition: TMath.h:1282
constexpr Double_t Ccgs()
Definition: TMath.h:125
Int_t Finite(Double_t x)
Check if it is finite with a mask in order to be consistent in presence of fast math.
Definition: TMath.h:770
const Double_t sigma
constexpr Double_t PiOver4()
Definition: TMath.h:59
Double_t ATan2(Double_t, Double_t)
Definition: TMath.h:678
Bool_t IsInside(T xp, T yp, Int_t np, T *x, T *y)
Function which returns kTRUE if point xp,yp lies inside the polygon defined by the np points in array...
Definition: TMath.h:1329
constexpr Double_t Pi()
Definition: TMath.h:38
double sin(double)
Double_t Infinity()
Returns an infinity as defined by the IEEE standard.
Definition: TMath.h:913
void BubbleLow(Int_t Narr, Double_t *arr1, Int_t *arr2)
void Error(const char *location, const char *msgfmt,...)
Double_t Erfc(Double_t x)
Compute the complementary error function erfc(x).
Definition: TMath.cxx:194
constexpr Double_t GnUncertainty()
Standard acceleration of gravity uncertainty.
Definition: TMath.h:181
Double_t DiLog(Double_t x)
Modified Struve functions of order 1.
Definition: TMath.cxx:110
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Definition: TMath.h:412
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Definition: TMath.h:416
#define F(x, y, z)
int isnan(double)
double gamma(double x)
constexpr Double_t Qe()
Elementary charge in .
Definition: TMath.h:336
constexpr Double_t G()
Gravitational constant in: .
Definition: TMath.h:139
ROOT::R::TRInterface & r
Definition: Object.C:4
Double_t VavilovI(Double_t x, Double_t kappa, Double_t beta2)
SVector< double, 2 > v
Definition: Dict.h:5
Double_t BesselY0(Double_t x)
Bessel function J1(x) for any real x.
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:184
auto * a
Definition: textangle.C:12
Double_t BetaDistI(Double_t x, Double_t p, Double_t q)
constexpr Double_t NaUncertainty()
Avogadro constant (Avogadro&#39;s Number) uncertainty.
Definition: TMath.h:290
constexpr Double_t RUncertainty()
Universal gas constant uncertainty.
Definition: TMath.h:306
Double_t StruveL1(Double_t x)
Modified Struve functions of order 0.
double cosh(double)
Double_t Voigt(Double_t x, Double_t sigma, Double_t lg, Int_t r=4)
Bool_t Permute(Int_t n, Int_t *a)
Double_t StruveL0(Double_t x)
Struve functions of order 1.
Double_t Mean(Long64_t n, const T *a, const Double_t *w=0)
Return the weighted mean of an array a with length n.
Definition: TMath.h:1098
unsigned int UInt_t
Definition: RtypesCore.h:42
Double_t KolmogorovTest(Int_t na, const Double_t *a, Int_t nb, const Double_t *b, Option_t *option)
Double_t Binomial(Int_t n, Int_t k)
double floor(double)
Double_t LandauI(Double_t x)
long double LongDouble_t
Definition: RtypesCore.h:57
Double_t ACos(Double_t)
Definition: TMath.h:667
Double_t BesselJ0(Double_t x)
modified Bessel function K_1(x)
static double p1(double t, double a, double b)
float xmax
Definition: THbookFile.cxx:93
Double_t ErfcInverse(Double_t x)
returns the inverse of the complementary error function x must be 0<x<2 implement using the quantile ...
Definition: TMath.cxx:237
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN) ...
Definition: TMath.h:907
T * Cross(const T v1[3], const T v2[3], T out[3])
Calculate the Cross Product of two vectors: out = [v1 x v2].
Definition: TMath.h:1295
double asin(double)
ULong_t Hash(const void *txt, Int_t ntxt)
constexpr Double_t Na()
Avogadro constant (Avogadro&#39;s Number) in .
Definition: TMath.h:283
constexpr Double_t HUncertainty()
Planck&#39;s constant uncertainty.
Definition: TMath.h:203
constexpr Double_t GUncertainty()
Gravitational constant uncertainty.
Definition: TMath.h:153
REAL epsilon
Definition: triangle.c:617
static T Min()
Returns maximum representation for type T.
Definition: TMath.h:920
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:97
Double_t Poisson(Double_t x, Double_t par)
Double_t Student(Double_t T, Double_t ndf)
Double_t Beta(Double_t p, Double_t q)
Double_t Cos(Double_t)
Definition: TMath.h:640
Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE)
Calculate a gaussian function with mean and sigma.
Definition: TMath.cxx:448
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t BesselI0(Double_t x)
integer order modified Bessel function K_n(x)
Double_t BesselK(Int_t n, Double_t x)
integer order modified Bessel function I_n(x)
constexpr Double_t HbarUncertainty()
uncertainty.
Definition: TMath.h:225
Double_t StdDev(Long64_t n, const T *a, const Double_t *w=0)
Definition: TMath.h:527
long Long_t
Definition: RtypesCore.h:50
RooCmdArg Index(RooCategory &icat)
#define d(i)
Definition: RSha256.hxx:102
Double_t Exp(Double_t x)
Definition: TMath.h:726
double Double_t
Definition: RtypesCore.h:55
Double_t StruveH0(Double_t x)
Bessel function Y1(x) for positive x.
double atan2(double, double)
int type
Definition: TGX11.cxx:120
Double_t Median(Long64_t n, const T *a, const Double_t *w=0, Long64_t *work=0)
Return the median of the array a where each entry i has weight w[i] .
Definition: TMath.h:1367
unsigned long ULong_t
Definition: RtypesCore.h:51
T MaxElement(Long64_t n, const T *a)
Return maximum of array a of length n.
Definition: TMath.h:958
Double_t y[n]
Definition: legend1.C:17
double atan(double)
static constexpr double s
Double_t Hypot(Double_t x, Double_t y)
Definition: TMath.cxx:57
constexpr Double_t C()
Velocity of light in .
Definition: TMath.h:118
constexpr Double_t Hbar()
in
Definition: TMath.h:211
Namespace for new Math classes and functions.
Double_t BesselK0(Double_t x)
modified Bessel function I_0(x)
Double_t BinomialI(Double_t p, Int_t n, Int_t k)
constexpr Double_t EulerGamma()
Euler-Mascheroni Constant.
Definition: TMath.h:329
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
constexpr Double_t Kcgs()
Definition: TMath.h:254
constexpr Double_t RadToDeg()
Conversion from radian to degree: .
Definition: TMath.h:74
constexpr Double_t CUncertainty()
Speed of light uncertainty.
Definition: TMath.h:132
Double_t BetaCf(Double_t x, Double_t a, Double_t b)
Element KOrdStat(Size n, const Element *a, Size k, Size *work=0)
Returns k_th order statistic of the array a of size n (k_th smallest element out of n elements)...
Definition: TMath.h:1453
auto * l
Definition: textangle.C:4
CompareDesc(T d)
Definition: TMath.h:1025
Double_t GammaDist(Double_t x, Double_t gamma, Double_t mu=0, Double_t beta=1)
double ceil(double)
int finite(double)
double tan(double)
Double_t Factorial(Int_t i)
Compute factorial(n).
Definition: TMath.cxx:247
Double_t Sin(Double_t)
Definition: TMath.h:636
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Double_t Ceil(Double_t x)
Definition: TMath.h:694
constexpr Double_t LogE()
Base-10 log of e (to convert ln to log)
Definition: TMath.h:111
Double_t ASin(Double_t)
Definition: TMath.h:660
Double_t BesselK1(Double_t x)
modified Bessel function I_1(x)
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition: TMath.cxx:486
constexpr Double_t Gcgs()
Definition: TMath.h:146
#define c(i)
Definition: RSha256.hxx:101
Double_t LogNormal(Double_t x, Double_t sigma, Double_t theta=0, Double_t m=1)
constexpr Double_t Sigma()
Stefan-Boltzmann constant in .
Definition: TMath.h:269
Double_t SinH(Double_t)
Definition: TMath.h:648
constexpr Double_t GhbarC()
in
Definition: TMath.h:160
Definition: first.py:1
Bool_t RootsCubic(const Double_t coef[4], Double_t &a, Double_t &b, Double_t &c)
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
double exp(double)
Double_t BetaIncomplete(Double_t x, Double_t a, Double_t b)
Long64_t LocMin(Long64_t n, const T *a)
Return index of array with the minimum element.
Definition: TMath.h:971
static T Max()
Returns minimum double representation.
Definition: TMath.h:927
float * q
Definition: THbookFile.cxx:87
TMath.
Definition: TMathBase.h:34
Double_t CauchyDist(Double_t x, Double_t t=0, Double_t s=1)
constexpr Double_t R()
Universal gas constant ( ) in
Definition: TMath.h:299
static T Epsilon()
Returns minimum double representation.
Definition: TMath.h:934
const Bool_t kTRUE
Definition: RtypesCore.h:87
Int_t Nint(T x)
Round to nearest integer. Rounds half integers to the nearest even integer.
Definition: TMath.h:712
const Int_t n
Definition: legend1.C:16
Double_t ATanH(Double_t)
Definition: TMath.cxx:90
Double_t CosH(Double_t)
Definition: TMath.h:652
Int_t CeilNint(Double_t x)
Definition: TMath.h:698
Double_t Tan(Double_t)
Definition: TMath.h:644
constexpr Double_t PiOver2()
Definition: TMath.h:52
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Binary search in an array of n values to locate value.
Definition: TMath.h:1221
double log(double)
double ldexp(double, int)
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE)
Sort the n1 elements of the Short_t array defined by its iterators.
Definition: TMath.h:1256
Double_t Vavilov(Double_t x, Double_t kappa, Double_t beta2)
constexpr Double_t KUncertainty()
Boltzmann&#39;s constant uncertainty.
Definition: TMath.h:261
T MinElement(Long64_t n, const T *a)
Return minimum of array a of length n.
Definition: TMath.h:951
Double_t ATan(Double_t)
Definition: TMath.h:674