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
16 //////////////////////////////////////////////////////////////////////////
17 // //
18 // TMath //
19 // //
20 // Encapsulate most frequently used Math functions. //
21 // NB. The basic functions Min, Max, Abs and Sign are defined //
22 // in TMathBase. //
23 // //
24 //////////////////////////////////////////////////////////////////////////
25
26 #ifndef ROOT_Rtypes
27 #include "Rtypes.h"
28 #endif
29 #ifndef ROOT_TMathBase
30 #include "TMathBase.h"
31 #endif
32
33 #include "TError.h"
34 #include <algorithm>
35 #include <limits>
36 #include <cmath>
37
38 namespace TMath {
39
40 /* ************************* */
41 /* * Fundamental constants * */
42 /* ************************* */
43
44 inline Double_t Pi() { return 3.14159265358979323846; }
45 inline Double_t TwoPi() { return 2.0 * Pi(); }
46 inline Double_t PiOver2() { return Pi() / 2.0; }
47 inline Double_t PiOver4() { return Pi() / 4.0; }
48 inline Double_t InvPi() { return 1.0 / Pi(); }
49 inline Double_t RadToDeg() { return 180.0 / Pi(); }
50 inline Double_t DegToRad() { return Pi() / 180.0; }
51 inline Double_t Sqrt2() { return 1.4142135623730950488016887242097; }
52
53 // e (base of natural log)
54 inline Double_t E() { return 2.71828182845904523536; }
55
56 // natural log of 10 (to convert log to ln)
57 inline Double_t Ln10() { return 2.30258509299404568402; }
58
59 // base-10 log of e (to convert ln to log)
60 inline Double_t LogE() { return 0.43429448190325182765; }
61
62 // velocity of light
63 inline Double_t C() { return 2.99792458e8; } // m s^-1
64 inline Double_t Ccgs() { return 100.0 * C(); } // cm s^-1
65 inline Double_t CUncertainty() { return 0.0; } // exact
66
67 // gravitational constant
68 inline Double_t G() { return 6.673e-11; } // m^3 kg^-1 s^-2
69 inline Double_t Gcgs() { return G() / 1000.0; } // cm^3 g^-1 s^-2
70 inline Double_t GUncertainty() { return 0.010e-11; }
71
72 // G over h-bar C
73 inline Double_t GhbarC() { return 6.707e-39; } // (GeV/c^2)^-2
74 inline Double_t GhbarCUncertainty() { return 0.010e-39; }
75
76 // standard acceleration of gravity
77 inline Double_t Gn() { return 9.80665; } // m s^-2
78 inline Double_t GnUncertainty() { return 0.0; } // exact
79
80 // Planck's constant
81 inline Double_t H() { return 6.62606876e-34; } // J s
82 inline Double_t Hcgs() { return 1.0e7 * H(); } // erg s
83 inline Double_t HUncertainty() { return 0.00000052e-34; }
84
85 // h-bar (h over 2 pi)
86 inline Double_t Hbar() { return 1.054571596e-34; } // J s
87 inline Double_t Hbarcgs() { return 1.0e7 * Hbar(); } // erg s
88 inline Double_t HbarUncertainty() { return 0.000000082e-34; }
89
90 // hc (h * c)
91 inline Double_t HC() { return H() * C(); } // J m
92 inline Double_t HCcgs() { return Hcgs() * Ccgs(); } // erg cm
93
94 // Boltzmann's constant
95 inline Double_t K() { return 1.3806503e-23; } // J K^-1
96 inline Double_t Kcgs() { return 1.0e7 * K(); } // erg K^-1
97 inline Double_t KUncertainty() { return 0.0000024e-23; }
98
99 // Stefan-Boltzmann constant
100 inline Double_t Sigma() { return 5.6704e-8; } // W m^-2 K^-4
101 inline Double_t SigmaUncertainty() { return 0.000040e-8; }
102
103 // Avogadro constant (Avogadro's Number)
104 inline Double_t Na() { return 6.02214199e+23; } // mol^-1
105 inline Double_t NaUncertainty() { return 0.00000047e+23; }
106
107 // universal gas constant (Na * K)
108 // http://scienceworld.wolfram.com/physics/UniversalGasConstant.html
109 inline Double_t R() { return K() * Na(); } // J K^-1 mol^-1
110 inline Double_t RUncertainty() { return R()*((KUncertainty()/K()) + (NaUncertainty()/Na())); }
111
112 // Molecular weight of dry air
113 // 1976 US Standard Atmosphere,
114 // also see http://atmos.nmsu.edu/jsdap/encyclopediawork.html
115 inline Double_t MWair() { return 28.9644; } // kg kmol^-1 (or gm mol^-1)
116
117 // Dry Air Gas Constant (R / MWair)
118 // http://atmos.nmsu.edu/education_and_outreach/encyclopedia/gas_constant.htm
119 inline Double_t Rgair() { return (1000.0 * R()) / MWair(); } // J kg^-1 K^-1
120
121 // Euler-Mascheroni Constant
122 inline Double_t EulerGamma() { return 0.577215664901532860606512090082402431042; }
123
124 // Elementary charge
125 inline Double_t Qe() { return 1.602176462e-19; } // C
126 inline Double_t QeUncertainty() { return 0.000000063e-19; }
127
128 /* ************************** */
129 /* * Mathematical Functions * */
130 /* ************************** */
131
132 /* ***************************** */
133 /* * Trigonometrical Functions * */
134 /* ***************************** */
135 inline Double_t Sin(Double_t);
136 inline Double_t Cos(Double_t);
137 inline Double_t Tan(Double_t);
138 inline Double_t SinH(Double_t);
139 inline Double_t CosH(Double_t);
140 inline Double_t TanH(Double_t);
141 inline Double_t ASin(Double_t);
142 inline Double_t ACos(Double_t);
143 inline Double_t ATan(Double_t);
144 inline Double_t ATan2(Double_t, Double_t);
145 Double_t ASinH(Double_t);
146 Double_t ACosH(Double_t);
147 Double_t ATanH(Double_t);
148 Double_t Hypot(Double_t x, Double_t y);
149
150
151 /* ************************ */
152 /* * Elementary Functions * */
153 /* ************************ */
154 inline Double_t Ceil(Double_t x);
155 inline Int_t CeilNint(Double_t x);
156 inline Double_t Floor(Double_t x);
157 inline Int_t FloorNint(Double_t x);
158 template<typename T>
159 inline Int_t Nint(T x);
160
161 inline Double_t Sqrt(Double_t x);
162 inline Double_t Exp(Double_t x);
163 inline Double_t Ldexp(Double_t x, Int_t exp);
164 Double_t Factorial(Int_t i);
165 inline LongDouble_t Power(LongDouble_t x, LongDouble_t y);
166 inline LongDouble_t Power(LongDouble_t x, Long64_t y);
167 inline LongDouble_t Power(Long64_t x, Long64_t y);
168 inline Double_t Power(Double_t x, Double_t y);
169 inline Double_t Power(Double_t x, Int_t y);
170 inline Double_t Log(Double_t x);
171 Double_t Log2(Double_t x);
172 inline Double_t Log10(Double_t x);
173 inline Int_t Finite(Double_t x);
174 inline Int_t IsNaN(Double_t x);
175
176 inline Double_t QuietNaN();
177 inline Double_t SignalingNaN();
178 inline Double_t Infinity();
179
180 template <typename T>
181 struct Limits {
182 inline static T Min();
183 inline static T Max();
184 inline static T Epsilon();
185 };
186
187 // Some integer math
188 Long_t Hypot(Long_t x, Long_t y); // sqrt(px*px + py*py)
189
190 // Comparing floating points
191 inline Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon) {
192 //return kTRUE if absolute difference between af and bf is less than epsilon
193 return TMath::Abs(af-bf) < epsilon;
194 }
195 inline Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec) {
196 //return kTRUE if relative difference between af and bf is less than relPrec
197 return TMath::Abs(af-bf) <= 0.5*relPrec*(TMath::Abs(af)+TMath::Abs(bf));
198 }
199
200 /* ******************** */
201 /* * Array Algorithms * */
202 /* ******************** */
203
204 // Min, Max of an array
205 template <typename T> T MinElement(Long64_t n, const T *a);
206 template <typename T> T MaxElement(Long64_t n, const T *a);
207
208 // Locate Min, Max element number in an array
209 template <typename T> Long64_t LocMin(Long64_t n, const T *a);
210 template <typename Iterator> Iterator LocMin(Iterator first, Iterator last);
211 template <typename T> Long64_t LocMax(Long64_t n, const T *a);
212 template <typename Iterator> Iterator LocMax(Iterator first, Iterator last);
213
214 // Binary search
215 template <typename T> Long64_t BinarySearch(Long64_t n, const T *array, T value);
216 template <typename T> Long64_t BinarySearch(Long64_t n, const T **array, T value);
217 template <typename Iterator, typename Element> Iterator BinarySearch(Iterator first, Iterator last, Element value);
218
219 // Hashing
220 ULong_t Hash(const void *txt, Int_t ntxt);
221 ULong_t Hash(const char *str);
222
223 // Sorting
224 template <typename Element, typename Index>
225 void Sort(Index n, const Element* a, Index* index, Bool_t down=kTRUE);
226 template <typename Iterator, typename IndexIterator>
227 void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE);
228
229 void BubbleHigh(Int_t Narr, Double_t *arr1, Int_t *arr2);
230 void BubbleLow (Int_t Narr, Double_t *arr1, Int_t *arr2);
231
232 Bool_t Permute(Int_t n, Int_t *a); // Find permutations
233
234 /* ************************* */
235 /* * Geometrical Functions * */
236 /* ************************* */
237
238 //Sample quantiles
239 void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob,
240 Bool_t isSorted=kTRUE, Int_t *index = 0, Int_t type=7);
241
242 // IsInside
243 template <typename T> Bool_t IsInside(T xp, T yp, Int_t np, T *x, T *y);
244
245 // Calculate the Cross Product of two vectors
246 template <typename T> T *Cross(const T v1[3],const T v2[3], T out[3]);
247
248 Float_t Normalize(Float_t v[3]); // Normalize a vector
249 Double_t Normalize(Double_t v[3]); // Normalize a vector
250
251 //Calculate the Normalized Cross Product of two vectors
252 template <typename T> inline T NormCross(const T v1[3],const T v2[3],T out[3]);
253
254 // Calculate a normal vector of a plane
255 template <typename T> T *Normal2Plane(const T v1[3],const T v2[3],const T v3[3], T normal[3]);
256
257 /* ************************ */
258 /* * Polynomial Functions * */
259 /* ************************ */
260
261 Bool_t RootsCubic(const Double_t coef[4],Double_t &a, Double_t &b, Double_t &c);
262
263 /* *********************** */
264 /* * Statistic Functions * */
265 /* *********************** */
266
267 Double_t Binomial(Int_t n,Int_t k); // Calculate the binomial coefficient n over k
268 Double_t BinomialI(Double_t p, Int_t n, Int_t k);
269 Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1);
270 Double_t CauchyDist(Double_t x, Double_t t=0, Double_t s=1);
271 Double_t ChisquareQuantile(Double_t p, Double_t ndf);
272 Double_t FDist(Double_t F, Double_t N, Double_t M);
273 Double_t FDistI(Double_t F, Double_t N, Double_t M);
274 Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE);
275 Double_t KolmogorovProb(Double_t z);
276 Double_t KolmogorovTest(Int_t na, const Double_t *a, Int_t nb, const Double_t *b, Option_t *option);
277 Double_t Landau(Double_t x, Double_t mpv=0, Double_t sigma=1, Bool_t norm=kFALSE);
278 Double_t LandauI(Double_t x);
279 Double_t LaplaceDist(Double_t x, Double_t alpha=0, Double_t beta=1);
280 Double_t LaplaceDistI(Double_t x, Double_t alpha=0, Double_t beta=1);
281 Double_t LogNormal(Double_t x, Double_t sigma, Double_t theta=0, Double_t m=1);
282 Double_t NormQuantile(Double_t p);
283 Double_t Poisson(Double_t x, Double_t par);
284 Double_t PoissonI(Double_t x, Double_t par);
285 Double_t Prob(Double_t chi2,Int_t ndf);
286 Double_t Student(Double_t T, Double_t ndf);
287 Double_t StudentI(Double_t T, Double_t ndf);
288 Double_t StudentQuantile(Double_t p, Double_t ndf, Bool_t lower_tail=kTRUE);
289 Double_t Vavilov(Double_t x, Double_t kappa, Double_t beta2);
290 Double_t VavilovI(Double_t x, Double_t kappa, Double_t beta2);
291 Double_t Voigt(Double_t x, Double_t sigma, Double_t lg, Int_t r = 4);
292
293 /* ************************** */
294 /* * Statistics over arrays * */
295 /* ************************** */
296
297 //Mean, Geometric Mean, Median, RMS(sigma)
298
299 template <typename T> Double_t Mean(Long64_t n, const T *a, const Double_t *w=0);
300 template <typename Iterator> Double_t Mean(Iterator first, Iterator last);
301 template <typename Iterator, typename WeightIterator> Double_t Mean(Iterator first, Iterator last, WeightIterator wfirst);
302
303 template <typename T> Double_t GeomMean(Long64_t n, const T *a);
304 template <typename Iterator> Double_t GeomMean(Iterator first, Iterator last);
305
306 template <typename T> Double_t RMS(Long64_t n, const T *a, const Double_t *w=0);
307 template <typename Iterator> Double_t RMS(Iterator first, Iterator last);
308 template <typename Iterator, typename WeightIterator> Double_t RMS(Iterator first, Iterator last, WeightIterator wfirst);
309
310 template <typename T> Double_t StdDev(Long64_t n, const T *a, const Double_t * w = 0) { return RMS<T>(n,a,w); }
311 template <typename Iterator> Double_t StdDev(Iterator first, Iterator last) { return RMS<Iterator>(first,last); }
312 template <typename Iterator, typename WeightIterator> Double_t StdDev(Iterator first, Iterator last, WeightIterator wfirst) { return RMS<Iterator,WeightIterator>(first,last,wfirst); }
313
314 template <typename T> Double_t Median(Long64_t n, const T *a, const Double_t *w=0, Long64_t *work=0);
315
316 //k-th order statistic
317 template <class Element, typename Size> Element KOrdStat(Size n, const Element *a, Size k, Size *work = 0);
318
319 /* ******************* */
320 /* * Special Functions */
321 /* ******************* */
322
323 Double_t Beta(Double_t p, Double_t q);
324 Double_t BetaCf(Double_t x, Double_t a, Double_t b);
325 Double_t BetaDist(Double_t x, Double_t p, Double_t q);
326 Double_t BetaDistI(Double_t x, Double_t p, Double_t q);
327 Double_t BetaIncomplete(Double_t x, Double_t a, Double_t b);
328
329 // Bessel functions
330 Double_t BesselI(Int_t n,Double_t x); // integer order modified Bessel function I_n(x)
331 Double_t BesselK(Int_t n,Double_t x); // integer order modified Bessel function K_n(x)
332 Double_t BesselI0(Double_t x); // modified Bessel function I_0(x)
333 Double_t BesselK0(Double_t x); // modified Bessel function K_0(x)
334 Double_t BesselI1(Double_t x); // modified Bessel function I_1(x)
335 Double_t BesselK1(Double_t x); // modified Bessel function K_1(x)
336 Double_t BesselJ0(Double_t x); // Bessel function J0(x) for any real x
337 Double_t BesselJ1(Double_t x); // Bessel function J1(x) for any real x
338 Double_t BesselY0(Double_t x); // Bessel function Y0(x) for positive x
339 Double_t BesselY1(Double_t x); // Bessel function Y1(x) for positive x
340 Double_t StruveH0(Double_t x); // Struve functions of order 0
341 Double_t StruveH1(Double_t x); // Struve functions of order 1
342 Double_t StruveL0(Double_t x); // Modified Struve functions of order 0
343 Double_t StruveL1(Double_t x); // Modified Struve functions of order 1
344
345 Double_t DiLog(Double_t x);
346 Double_t Erf(Double_t x);
347 Double_t ErfInverse(Double_t x);
348 Double_t Erfc(Double_t x);
349 Double_t ErfcInverse(Double_t x);
350 Double_t Freq(Double_t x);
351 Double_t Gamma(Double_t z);
352 Double_t Gamma(Double_t a,Double_t x);
353 Double_t GammaDist(Double_t x, Double_t gamma, Double_t mu=0, Double_t beta=1);
354 Double_t LnGamma(Double_t z);
355 }
356
357
358 //---- Trig and other functions ------------------------------------------------
359
360 #include <float.h>
361
362 #if defined(R__WIN32) && !defined(__CINT__)
363 # ifndef finite
364 # define finite _finite
365 # define isnan _isnan
366 # endif
367 #endif
368 #if defined(R__AIX) || defined(R__SOLARIS_CC50) || \
369 defined(R__HPUX11) || defined(R__GLIBC) || \
370 (defined(R__MACOSX) && (defined(__INTEL_COMPILER) || defined(__arm__)))
371 // math functions are defined inline so we have to include them here
372 # include <math.h>
373 # ifdef R__SOLARIS_CC50
374 extern "C" { int finite(double); }
375 # endif
376 // # if defined(R__GLIBC) && defined(__STRICT_ANSI__)
377 // # ifndef finite
378 // # define finite __finite
379 // # endif
380 // # ifndef isnan
381 // # define isnan __isnan
382 // # endif
383 // # endif
384 #else
385 // don't want to include complete <math.h>
386 extern "C" {
387 extern double sin(double);
388 extern double cos(double);
389 extern double tan(double);
390 extern double sinh(double);
391 extern double cosh(double);
392 extern double tanh(double);
393 extern double asin(double);
394 extern double acos(double);
395 extern double atan(double);
396 extern double atan2(double, double);
397 extern double sqrt(double);
398 extern double exp(double);
399 extern double pow(double, double);
400 extern double log(double);
401 extern double log10(double);
402 #ifndef R__WIN32
403 # if !defined(finite)
404 extern int finite(double);
405 # endif
406 # if !defined(isnan)
407 extern int isnan(double);
408 # endif
409 extern double ldexp(double, int);
410 extern double ceil(double);
411 extern double floor(double);
412 #else
413 _CRTIMP double ldexp(double, int);
414 _CRTIMP double ceil(double);
415 _CRTIMP double floor(double);
416 #endif
417 }
418 #endif
419
420 inline Double_t TMath::Sin(Double_t x)
421 { return sin(x); }
422
423 inline Double_t TMath::Cos(Double_t x)
424 { return cos(x); }
425
426 inline Double_t TMath::Tan(Double_t x)
427 { return tan(x); }
428
429 inline Double_t TMath::SinH(Double_t x)
430 { return sinh(x); }
431
432 inline Double_t TMath::CosH(Double_t x)
433 { return cosh(x); }
434
435 inline Double_t TMath::TanH(Double_t x)
436 { return tanh(x); }
437
438 inline Double_t TMath::ASin(Double_t x)
439 { if (x < -1.) return -TMath::Pi()/2;
440 if (x > 1.) return TMath::Pi()/2;
441 return asin(x);
442 }
443
444 inline Double_t TMath::ACos(Double_t x)
445 { if (x < -1.) return TMath::Pi();
446 if (x > 1.) return 0;
447 return acos(x);
448 }
449
450 inline Double_t TMath::ATan(Double_t x)
451 { return atan(x); }
452
453 inline Double_t TMath::ATan2(Double_t y, Double_t x)
454 { if (x != 0) return atan2(y, x);
455 if (y == 0) return 0;
456 if (y > 0) return Pi()/2;
457 else return -Pi()/2;
458 }
459
460 inline Double_t TMath::Sqrt(Double_t x)
461 { return sqrt(x); }
462
463 inline Double_t TMath::Ceil(Double_t x)
464 { return ceil(x); }
465
466 inline Int_t TMath::CeilNint(Double_t x)
467 { return TMath::Nint(ceil(x)); }
468
469 inline Double_t TMath::Floor(Double_t x)
470 { return floor(x); }
471
472 inline Int_t TMath::FloorNint(Double_t x)
473 { return TMath::Nint(floor(x)); }
474
475 template<typename T>
476 inline Int_t TMath::Nint(T x)
477 {
478 // Round to nearest integer. Rounds half integers to the nearest
479 // even integer.
480 int i;
481 if (x >= 0) {
482 i = int(x + 0.5);
483 if ( i & 1 && x + 0.5 == T(i) ) i--;
484 } else {
485 i = int(x - 0.5);
486 if ( i & 1 && x - 0.5 == T(i) ) i++;
487 }
488 return i;
489 }
490
491 inline Double_t TMath::Exp(Double_t x)
492 { return exp(x); }
493
494 inline Double_t TMath::Ldexp(Double_t x, Int_t exp)
495 { return ldexp(x, exp); }
496
497 inline LongDouble_t TMath::Power(LongDouble_t x, LongDouble_t y)
498 { return std::pow(x,y); }
499
500 inline LongDouble_t TMath::Power(LongDouble_t x, Long64_t y)
501 { return std::pow(x,(LongDouble_t)y); }
502
503 inline LongDouble_t TMath::Power(Long64_t x, Long64_t y)
504 #if __cplusplus >= 201103 /*C++11*/
505 { return std::pow(x,y); }
506 #else
507 { return std::pow((LongDouble_t)x,(LongDouble_t)y); }
508 #endif
509
510 inline Double_t TMath::Power(Double_t x, Double_t y)
511 { return pow(x, y); }
512
513 inline Double_t TMath::Power(Double_t x, Int_t y) {
514 #ifdef R__ANSISTREAM
515 return std::pow(x, y);
516 #else
517 return pow(x, (Double_t) y);
518 #endif
519 }
520
521
522 inline Double_t TMath::Log(Double_t x)
523 { return log(x); }
524
525 inline Double_t TMath::Log10(Double_t x)
526 { return log10(x); }
527
528 inline Int_t TMath::Finite(Double_t x)
529 #if defined(R__HPUX11)
530 { return isfinite(x); }
531 #elif defined(R__MACOSX)
532 #ifdef isfinite
533 // from math.h
534 { return isfinite(x); }
535 #else
536 // from cmath
537 { return std::isfinite(x); }
538 #endif
539 #else
540 { return finite(x); }
541 #endif
542
543 inline Int_t TMath::IsNaN(Double_t x)
544 #if (defined(R__ANSISTREAM) || (defined(R__MACOSX) && defined(__arm__))) && !defined(_AIX) && !defined(__CUDACC__)
545 #if defined(isnan) || defined(R__SOLARIS_CC50) || defined(__INTEL_COMPILER)
546 // from math.h
547 { return ::isnan(x); }
548 #else
549 // from cmath
550 { return std::isnan(x); }
551 #endif
552 #else
553 { return isnan(x); }
554 #endif
555
556 //--------wrapper to numeric_limits
557 //____________________________________________________________________________
558 inline Double_t TMath::QuietNaN() {
559 // returns a quiet NaN as defined by IEEE 754
560 // see http://en.wikipedia.org/wiki/NaN#Quiet_NaN
561 return std::numeric_limits<Double_t>::quiet_NaN();
562 }
563
564 //____________________________________________________________________________
565 inline Double_t TMath::SignalingNaN() {
566 // returns a signaling NaN as defined by IEEE 754
567 // see http://en.wikipedia.org/wiki/NaN#Signaling_NaN
568 return std::numeric_limits<Double_t>::signaling_NaN();
569 }
570
571 inline Double_t TMath::Infinity() {
572 // returns an infinity as defined by the IEEE standard
573 return std::numeric_limits<Double_t>::infinity();
574 }
575
576 template<typename T>
577 inline T TMath::Limits<T>::Min() {
578 // returns maximum representation for type T
579 return (std::numeric_limits<T>::min)(); //N.B. use this signature to avoid class with macro min() on Windows
580 }
581
582 template<typename T>
583 inline T TMath::Limits<T>::Max() {
584 // returns minimum double representation
585 return (std::numeric_limits<T>::max)(); //N.B. use this signature to avoid class with macro max() on Windows
586 }
587
588 template<typename T>
589 inline T TMath::Limits<T>::Epsilon() {
590 // returns minimum double representation
591 return std::numeric_limits<T>::epsilon();
592 }
593
594
595 //-------- Advanced -------------
596
597 template <typename T> inline T TMath::NormCross(const T v1[3],const T v2[3],T out[3])
598 {
599 // Calculate the Normalized Cross Product of two vectors
600 return Normalize(Cross(v1,v2,out));
601 }
602
603 template <typename T>
604 T TMath::MinElement(Long64_t n, const T *a) {
605 // Return minimum of array a of length n.
606
607 return *std::min_element(a,a+n);
608 }
609
610 template <typename T>
611 T TMath::MaxElement(Long64_t n, const T *a) {
612 // Return maximum of array a of length n.
613
614 return *std::max_element(a,a+n);
615 }
616
617 template <typename T>
618 Long64_t TMath::LocMin(Long64_t n, const T *a) {
619 // Return index of array with the minimum element.
620 // If more than one element is minimum returns first found.
621
622 // Implement here since this one is found to be faster (mainly on 64 bit machines)
623 // than stl generic implementation.
624 // When performing the comparison, the STL implementation needs to de-reference both the array iterator
625 // and the iterator pointing to the resulting minimum location
626
627 if (n <= 0 || !a) return -1;
628 T xmin = a[0];
629 Long64_t loc = 0;
630 for (Long64_t i = 1; i < n; i++) {
631 if (xmin > a[i]) {
632 xmin = a[i];
633 loc = i;
634 }
635 }
636 return loc;
637 }
638
639 template <typename Iterator>
640 Iterator TMath::LocMin(Iterator first, Iterator last) {
641 // Return index of array with the minimum element.
642 // If more than one element is minimum returns first found.
643 return std::min_element(first, last);
644 }
645
646 template <typename T>
647 Long64_t TMath::LocMax(Long64_t n, const T *a) {
648 // Return index of array with the maximum element.
649 // If more than one element is maximum returns first found.
650
651 // Implement here since it is faster (see comment in LocMin function)
652
653 if (n <= 0 || !a) return -1;
654 T xmax = a[0];
655 Long64_t loc = 0;
656 for (Long64_t i = 1; i < n; i++) {
657 if (xmax < a[i]) {
658 xmax = a[i];
659 loc = i;
660 }
661 }
662 return loc;
663 }
664
665 template <typename Iterator>
666 Iterator TMath::LocMax(Iterator first, Iterator last)
667 {
668 // Return index of array with the maximum element.
669 // If more than one element is maximum returns first found.
670
671 return std::max_element(first, last);
672 }
673
674 template<typename T>
675 struct CompareDesc {
676
677 CompareDesc(T d) : fData(d) {}
678
679 template<typename Index>
680 bool operator()(Index i1, Index i2) {
681 return *(fData + i1) > *(fData + i2);
682 }
683
684 T fData;
685 };
686
687 template<typename T>
688 struct CompareAsc {
689
690 CompareAsc(T d) : fData(d) {}
691
692 template<typename Index>
693 bool operator()(Index i1, Index i2) {
694 return *(fData + i1) < *(fData + i2);
695 }
696
697 T fData;
698 };
699
700 template <typename Iterator>
701 Double_t TMath::Mean(Iterator first, Iterator last)
702 {
703 // Return the weighted mean of an array defined by the iterators.
704
705 Double_t sum = 0;
706 Double_t sumw = 0;
707 while ( first != last )
708 {
709 sum += *first;
710 sumw += 1;
711 first++;
712 }
713
714 return sum/sumw;
715 }
716
717 template <typename Iterator, typename WeightIterator>
718 Double_t TMath::Mean(Iterator first, Iterator last, WeightIterator w)
719 {
720 // Return the weighted mean of an array defined by the first and
721 // last iterators. The w iterator should point to the first element
722 // of a vector of weights of the same size as the main array.
723
724 Double_t sum = 0;
725 Double_t sumw = 0;
726 int i = 0;
727 while ( first != last ) {
728 if ( *w < 0) {
729 ::Error("TMath::Mean","w[%d] = %.4e < 0 ?!",i,*w);
730 return 0;
731 }
732 sum += (*w) * (*first);
733 sumw += (*w) ;
734 ++w;
735 ++first;
736 ++i;
737 }
738 if (sumw <= 0) {
739 ::Error("TMath::Mean","sum of weights == 0 ?!");
740 return 0;
741 }
742
743 return sum/sumw;
744 }
745
746 template <typename T>
747 Double_t TMath::Mean(Long64_t n, const T *a, const Double_t *w)
748 {
749 // Return the weighted mean of an array a with length n.
750
751 if (w) {
752 return TMath::Mean(a, a+n, w);
753 } else {
754 return TMath::Mean(a, a+n);
755 }
756 }
757
758 template <typename Iterator>
759 Double_t TMath::GeomMean(Iterator first, Iterator last)
760 {
761 // Return the geometric mean of an array defined by the iterators.
762 // geometric_mean = (Prod_i=0,n-1 |a[i]|)^1/n
763
764 Double_t logsum = 0.;
765 Long64_t n = 0;
766 while ( first != last ) {
767 if (*first == 0) return 0.;
768 Double_t absa = (Double_t) TMath::Abs(*first);
769 logsum += TMath::Log(absa);
770 ++first;
771 ++n;
772 }
773
774 return TMath::Exp(logsum/n);
775 }
776
777 template <typename T>
778 Double_t TMath::GeomMean(Long64_t n, const T *a)
779 {
780 // Return the geometric mean of an array a of size n.
781 // geometric_mean = (Prod_i=0,n-1 |a[i]|)^1/n
782
783 return TMath::GeomMean(a, a+n);
784 }
785
786 template <typename Iterator>
787 Double_t TMath::RMS(Iterator first, Iterator last)
788 {
789 // Return the Standard Deviation of an array defined by the iterators.
790 // Note that this function returns the sigma(standard deviation) and
791 // not the root mean square of the array.
792
793 // Use the two pass algorithm, which is slower (! a factor of 2) but much more
794 // precise. Since we have a vector the 2 pass algorithm is still faster than the
795 // Welford algorithm. (See also ROOT-5545)
796
797 Double_t n = 0;
798
799 Double_t tot = 0;
800 Double_t mean = TMath::Mean(first,last);
801 while ( first != last ) {
802 Double_t x = Double_t(*first);
803 tot += (x - mean)*(x - mean);
804 ++first;
805 ++n;
806 }
807 Double_t rms = (n > 1) ? TMath::Sqrt(tot/(n-1)) : 0.0;
808 return rms;
809 }
810
811 template <typename Iterator, typename WeightIterator>
812 Double_t TMath::RMS(Iterator first, Iterator last, WeightIterator w)
813 {
814 // Return the weighted Standard Deviation of an array defined by the iterators.
815 // Note that this function returns the sigma(standard deviation) and
816 // not the root mean square of the array.
817
818 // As in the unweighted case use the two pass algorithm
819
820 Double_t tot = 0;
821 Double_t sumw = 0;
822 Double_t sumw2 = 0;
823 Double_t mean = TMath::Mean(first,last,w);
824 while ( first != last ) {
825 Double_t x = Double_t(*first);
826 sumw += *w;
827 sumw2 += (*w) * (*w);
828 tot += (*w) * (x - mean)*(x - mean);
829 ++first;
830 ++w;
831 }
832 // use the correction neff/(neff -1) for the unbiased formula
833 Double_t rms = TMath::Sqrt(tot * sumw/ (sumw*sumw - sumw2) );
834 return rms;
835 }
836
837
838 template <typename T>
839 Double_t TMath::RMS(Long64_t n, const T *a, const Double_t * w)
840 {
841 // Return the Standard Deviation of an array a with length n.
842 // Note that this function returns the sigma(standard deviation) and
843 // not the root mean square of the array.
844
845 return (w) ? TMath::RMS(a, a+n, w) : TMath::RMS(a, a+n);
846 }
847
848 template <typename Iterator, typename Element>
849 Iterator TMath::BinarySearch(Iterator first, Iterator last, Element value)
850 {
851 // Binary search in an array defined by its iterators.
852 //
853 // The values in the iterators range are supposed to be sorted
854 // prior to this call. If match is found, function returns
855 // position of element. If no match found, function gives nearest
856 // element smaller than value.
857
858 Iterator pind;
859 pind = std::lower_bound(first, last, value);
860 if ( (pind != last) && (*pind == value) )
861 return pind;
862 else
863 return ( pind - 1);
864 }
865
866
867 template <typename T> Long64_t TMath::BinarySearch(Long64_t n, const T *array, T value)
868 {
869 // Binary search in an array of n values to locate value.
870 //
871 // Array is supposed to be sorted prior to this call.
872 // If match is found, function returns position of element.
873 // If no match found, function gives nearest element smaller than value.
874
875 const T* pind;
876 pind = std::lower_bound(array, array + n, value);
877 if ( (pind != array + n) && (*pind == value) )
878 return (pind - array);
879 else
880 return ( pind - array - 1);
881 }
882
883 template <typename T> Long64_t TMath::BinarySearch(Long64_t n, const T **array, T value)
884 {
885 // Binary search in an array of n values to locate value.
886 //
887 // Array is supposed to be sorted prior to this call.
888 // If match is found, function returns position of element.
889 // If no match found, function gives nearest element smaller than value.
890
891 const T* pind;
892 pind = std::lower_bound(*array, *array + n, value);
893 if ( (pind != *array + n) && (*pind == value) )
894 return (pind - *array);
895 else
896 return ( pind - *array - 1);
897 }
898
899 template <typename Iterator, typename IndexIterator>
900 void TMath::SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down)
901 {
902 // Sort the n1 elements of the Short_t array defined by its
903 // iterators. In output the array index contains the indices of
904 // the sorted array. If down is false sort in increasing order
905 // (default is decreasing order).
906
907 // NOTE that the array index must be created with a length bigger
908 // or equal than the main array before calling this function.
909
910 int i = 0;
911
912 IndexIterator cindex = index;
913 for ( Iterator cfirst = first; cfirst != last; ++cfirst )
914 {
915 *cindex = i++;
916 ++cindex;
917 }
918
919 if ( down )
920 std::sort(index, cindex, CompareDesc<Iterator>(first) );
921 else
922 std::sort(index, cindex, CompareAsc<Iterator>(first) );
923 }
924
925 template <typename Element, typename Index> void TMath::Sort(Index n, const Element* a, Index* index, Bool_t down)
926 {
927 // Sort the n elements of the array a of generic templated type Element.
928 // In output the array index of type Index contains the indices of the sorted array.
929 // If down is false sort in increasing order (default is decreasing order).
930
931 // NOTE that the array index must be created with a length >= n
932 // before calling this function.
933 // NOTE also that the size type for n must be the same type used for the index array
934 // (templated type Index)
935
936 for(Index i = 0; i < n; i++) { index[i] = i; }
937 if ( down )
938 std::sort(index, index + n, CompareDesc<const Element*>(a) );
939 else
940 std::sort(index, index + n, CompareAsc<const Element*>(a) );
941 }
942
943 template <typename T> T *TMath::Cross(const T v1[3],const T v2[3], T out[3])
944 {
945 // Calculate the Cross Product of two vectors:
946 // out = [v1 x v2]
947
948 out[0] = v1[1] * v2[2] - v1[2] * v2[1];
949 out[1] = v1[2] * v2[0] - v1[0] * v2[2];
950 out[2] = v1[0] * v2[1] - v1[1] * v2[0];
951
952 return out;
953 }
954
955 template <typename T> T * TMath::Normal2Plane(const T p1[3],const T p2[3],const T p3[3], T normal[3])
956 {
957 // Calculate a normal vector of a plane.
958 //
959 // Input:
960 // Float_t *p1,*p2,*p3 - 3 3D points belonged the plane to define it.
961 //
962 // Return:
963 // Pointer to 3D normal vector (normalized)
964
965 T v1[3], v2[3];
966
967 v1[0] = p2[0] - p1[0];
968 v1[1] = p2[1] - p1[1];
969 v1[2] = p2[2] - p1[2];
970
971 v2[0] = p3[0] - p1[0];
972 v2[1] = p3[1] - p1[1];
973 v2[2] = p3[2] - p1[2];
974
975 NormCross(v1,v2,normal);
976 return normal;
977 }
978
979 template <typename T> Bool_t TMath::IsInside(T xp, T yp, Int_t np, T *x, T *y)
980 {
981 // Function which returns kTRUE if point xp,yp lies inside the
982 // polygon defined by the np points in arrays x and y, kFALSE otherwise.
983 // Note that the polygon may be open or closed.
984
985 Int_t i, j = np-1 ;
986 Bool_t oddNodes = kFALSE;
987
988 for (i=0; i<np; i++) {
989 if ((y[i]<yp && y[j]>=yp) || (y[j]<yp && y[i]>=yp)) {
990 if (x[i]+(yp-y[i])/(y[j]-y[i])*(x[j]-x[i])<xp) {
991 oddNodes = !oddNodes;
992 }
993 }
994 j=i;
995 }
996
997 return oddNodes;
998 }
999
1000 template <typename T> Double_t TMath::Median(Long64_t n, const T *a, const Double_t *w, Long64_t *work)
1001 {
1002 // Return the median of the array a where each entry i has weight w[i] .
1003 // Both arrays have a length of at least n . The median is a number obtained
1004 // from the sorted array a through
1005 //
1006 // median = (a[jl]+a[jh])/2. where (using also the sorted index on the array w)
1007 //
1008 // sum_i=0,jl w[i] <= sumTot/2
1009 // sum_i=0,jh w[i] >= sumTot/2
1010 // sumTot = sum_i=0,n w[i]
1011 //
1012 // If w=0, the algorithm defaults to the median definition where it is
1013 // a number that divides the sorted sequence into 2 halves.
1014 // When n is odd or n > 1000, the median is kth element k = (n + 1) / 2.
1015 // when n is even and n < 1000the median is a mean of the elements k = n/2 and k = n/2 + 1.
1016 //
1017 // If the weights are supplied (w not 0) all weights must be >= 0
1018 //
1019 // If work is supplied, it is used to store the sorting index and assumed to be
1020 // >= n . If work=0, local storage is used, either on the stack if n < kWorkMax
1021 // or on the heap for n >= kWorkMax .
1022
1023 const Int_t kWorkMax = 100;
1024
1025 if (n <= 0 || !a) return 0;
1026 Bool_t isAllocated = kFALSE;
1027 Double_t median;
1028 Long64_t *ind;
1029 Long64_t workLocal[kWorkMax];
1030
1031 if (work) {
1032 ind = work;
1033 } else {
1034 ind = workLocal;
1035 if (n > kWorkMax) {
1036 isAllocated = kTRUE;
1037 ind = new Long64_t[n];
1038 }
1039 }
1040
1041 if (w) {
1042 Double_t sumTot2 = 0;
1043 for (Int_t j = 0; j < n; j++) {
1044 if (w[j] < 0) {
1045 ::Error("TMath::Median","w[%d] = %.4e < 0 ?!",j,w[j]);
1046 if (isAllocated) delete [] ind;
1047 return 0;
1048 }
1049 sumTot2 += w[j];
1050 }
1051
1052 sumTot2 /= 2.;
1053
1054 Sort(n, a, ind, kFALSE);
1055
1056 Double_t sum = 0.;
1057 Int_t jl;
1058 for (jl = 0; jl < n; jl++) {
1059 sum += w[ind[jl]];
1060 if (sum >= sumTot2) break;
1061 }
1062
1063 Int_t jh;
1064 sum = 2.*sumTot2;
1065 for (jh = n-1; jh >= 0; jh--) {
1066 sum -= w[ind[jh]];
1067 if (sum <= sumTot2) break;
1068 }
1069
1070 median = 0.5*(a[ind[jl]]+a[ind[jh]]);
1071
1072 } else {
1073
1074 if (n%2 == 1)
1075 median = KOrdStat(n, a,n/2, ind);
1076 else {
1077 median = 0.5*(KOrdStat(n, a, n/2 -1, ind)+KOrdStat(n, a, n/2, ind));
1078 }
1079 }
1080
1081 if (isAllocated)
1082 delete [] ind;
1083 return median;
1084 }
1085
1086
1087
1088
1089 template <class Element, typename Size>
1090 Element TMath::KOrdStat(Size n, const Element *a, Size k, Size *work)
1091 {
1092 // Returns k_th order statistic of the array a of size n
1093 // (k_th smallest element out of n elements).
1094 //
1095 // C-convention is used for array indexing, so if you want
1096 // the second smallest element, call KOrdStat(n, a, 1).
1097 //
1098 // If work is supplied, it is used to store the sorting index and
1099 // assumed to be >= n. If work=0, local storage is used, either on
1100 // the stack if n < kWorkMax or on the heap for n >= kWorkMax.
1101 // Note that the work index array will not contain the sorted indices but
1102 // all indeces of the smaller element in arbitrary order in work[0,...,k-1] and
1103 // all indeces of the larger element in arbitrary order in work[k+1,..,n-1]
1104 // work[k] will contain instead the index of the returned element.
1105 //
1106 // Taken from "Numerical Recipes in C++" without the index array
1107 // implemented by Anna Khreshuk.
1108 //
1109 // See also the declarations at the top of this file
1110
1111 const Int_t kWorkMax = 100;
1112
1113 typedef Size Index;
1114
1115 Bool_t isAllocated = kFALSE;
1116 Size i, ir, j, l, mid;
1117 Index arr;
1118 Index *ind;
1119 Index workLocal[kWorkMax];
1120 Index temp;
1121
1122 if (work) {
1123 ind = work;
1124 } else {
1125 ind = workLocal;
1126 if (n > kWorkMax) {
1127 isAllocated = kTRUE;
1128 ind = new Index[n];
1129 }
1130 }
1131
1132 for (Size ii=0; ii<n; ii++) {
1133 ind[ii]=ii;
1134 }
1135 Size rk = k;
1136 l=0;
1137 ir = n-1;
1138 for(;;) {
1139 if (ir<=l+1) { //active partition contains 1 or 2 elements
1140 if (ir == l+1 && a[ind[ir]]<a[ind[l]])
1141 {temp = ind[l]; ind[l]=ind[ir]; ind[ir]=temp;}
1142 Element tmp = a[ind[rk]];
1143 if (isAllocated)
1144 delete [] ind;
1145 return tmp;
1146 } else {
1147 mid = (l+ir) >> 1; //choose median of left, center and right
1148 {temp = ind[mid]; ind[mid]=ind[l+1]; ind[l+1]=temp;}//elements as partitioning element arr.
1149 if (a[ind[l]]>a[ind[ir]]) //also rearrange so that a[l]<=a[l+1]
1150 {temp = ind[l]; ind[l]=ind[ir]; ind[ir]=temp;}
1151
1152 if (a[ind[l+1]]>a[ind[ir]])
1153 {temp=ind[l+1]; ind[l+1]=ind[ir]; ind[ir]=temp;}
1154
1155 if (a[ind[l]]>a[ind[l+1]])
1156 {temp = ind[l]; ind[l]=ind[l+1]; ind[l+1]=temp;}
1157
1158 i=l+1; //initialize pointers for partitioning
1159 j=ir;
1160 arr = ind[l+1];
1161 for (;;){
1162 do i++; while (a[ind[i]]<a[arr]);
1163 do j--; while (a[ind[j]]>a[arr]);
1164 if (j<i) break; //pointers crossed, partitioning complete
1165 {temp=ind[i]; ind[i]=ind[j]; ind[j]=temp;}
1166 }
1167 ind[l+1]=ind[j];
1168 ind[j]=arr;
1169 if (j>=rk) ir = j-1; //keep active the partition that
1170 if (j<=rk) l=i; //contains the k_th element
1171 }
1172 }
1173 }
1174
1175 #endif
1176