Logo ROOT  
Reference Guide
VectorizedTMath.cxx
Go to the documentation of this file.
1 #include "VectorizedTMath.h"
2 
3 #if defined(R__HAS_VECCORE) && defined(R__HAS_VC)
4 
5 namespace TMath {
6 ////////////////////////////////////////////////////////////////////////////////
8 {
9  return vecCore::math::Log2(x);
10 }
11 
12 ////////////////////////////////////////////////////////////////////////////////
13 /// Calculate a Breit Wigner function with mean and gamma.
15 {
16  return 0.5 * M_1_PI * (gamma / (0.25 * gamma * gamma + (x - mean) * (x - mean)));
17 }
18 
19 ////////////////////////////////////////////////////////////////////////////////
20 /// Calculate a gaussian function with mean and sigma.
21 /// If norm=kTRUE (default is kFALSE) the result is divided
22 /// by sqrt(2*Pi)*sigma.
24 {
25  if (sigma == 0)
27 
28  ::ROOT::Double_v inv_sigma = 1.0 / ::ROOT::Double_v(sigma);
29  ::ROOT::Double_v arg = (x - ::ROOT::Double_v(mean)) * inv_sigma;
30 
31  // For those entries of |arg| > 39 result is zero in double precision
32  ::ROOT::Double_v out =
33  vecCore::Blend<::ROOT::Double_v>(vecCore::math::Abs(arg) < ::ROOT::Double_v(39.0),
34  vecCore::math::Exp(::ROOT::Double_v(-0.5) * arg * arg), ::ROOT::Double_v(0.0));
35  if (norm)
36  out *= 0.3989422804014327 * inv_sigma; // 1/sqrt(2*Pi)=0.3989422804014327
37  return out;
38 }
39 
40 ////////////////////////////////////////////////////////////////////////////////
41 /// Computes the probability density function of Laplace distribution
42 /// at point x, with location parameter alpha and shape parameter beta.
43 /// By default, alpha=0, beta=1
44 /// This distribution is known under different names, most common is
45 /// double exponential distribution, but it also appears as
46 /// the two-tailed exponential or the bilateral exponential distribution
48 {
50  ::ROOT::Double_v out = vecCore::math::Exp(-vecCore::math::Abs((x - ::ROOT::Double_v(alpha)) * beta_v_inv));
51  out *= ::ROOT::Double_v(0.5) * beta_v_inv;
52  return out;
53 }
54 
55 ////////////////////////////////////////////////////////////////////////////////
56 /// Computes the distribution function of Laplace distribution
57 /// at point x, with location parameter alpha and shape parameter beta.
58 /// By default, alpha=0, beta=1
59 /// This distribution is known under different names, most common is
60 /// double exponential distribution, but it also appears as
61 /// the two-tailed exponential or the bilateral exponential distribution
63 {
64  ::ROOT::Double_v alpha_v = ::ROOT::Double_v(alpha);
66  return vecCore::Blend<::ROOT::Double_v>(
67  x <= alpha_v, 0.5 * vecCore::math::Exp(-vecCore::math::Abs((x - alpha_v) * beta_v_inv)),
68  1 - 0.5 * vecCore::math::Exp(-vecCore::math::Abs((x - alpha_v) * beta_v_inv)));
69 }
70 
71 ////////////////////////////////////////////////////////////////////////////////
72 /// Computation of the normal frequency function freq(x).
73 /// Freq(x) = (1/sqrt(2pi)) Integral(exp(-t^2/2))dt between -infinity and x.
74 ///
75 /// Translated from CERNLIB C300 by Rene Brun.
77 {
78  Double_t c1 = 0.56418958354775629;
79  Double_t w2 = 1.41421356237309505;
80 
81  Double_t p10 = 2.4266795523053175e+2, q10 = 2.1505887586986120e+2, p11 = 2.1979261618294152e+1,
82  q11 = 9.1164905404514901e+1, p12 = 6.9963834886191355e+0, q12 = 1.5082797630407787e+1,
83  p13 = -3.5609843701815385e-2;
84 
85  Double_t p20 = 3.00459261020161601e+2, q20 = 3.00459260956983293e+2, p21 = 4.51918953711872942e+2,
86  q21 = 7.90950925327898027e+2, p22 = 3.39320816734343687e+2, q22 = 9.31354094850609621e+2,
87  p23 = 1.52989285046940404e+2, q23 = 6.38980264465631167e+2, p24 = 4.31622272220567353e+1,
88  q24 = 2.77585444743987643e+2, p25 = 7.21175825088309366e+0, q25 = 7.70001529352294730e+1,
89  p26 = 5.64195517478973971e-1, q26 = 1.27827273196294235e+1, p27 = -1.36864857382716707e-7;
90 
91  Double_t p30 = -2.99610707703542174e-3, q30 = 1.06209230528467918e-2, p31 = -4.94730910623250734e-2,
92  q31 = 1.91308926107829841e-1, p32 = -2.26956593539686930e-1, q32 = 1.05167510706793207e+0,
93  p33 = -2.78661308609647788e-1, q33 = 1.98733201817135256e+0, p34 = -2.23192459734184686e-2, q34 = 1;
94 
96 
97  ::ROOT::Double_v ap, aq, h, hc, y, result;
98 
99  vecCore::Mask<::ROOT::Double_v> mask1 = v < ::ROOT::Double_v(0.5);
100  vecCore::Mask<::ROOT::Double_v> mask2 = !mask1 && v < ::ROOT::Double_v(4.0);
101  vecCore::Mask<::ROOT::Double_v> mask3 = !(mask1 || mask2);
102 
103  ::ROOT::Double_v v2 = v * v;
104  ::ROOT::Double_v v3 = v2 * v;
105  ::ROOT::Double_v v4 = v3 * v;
106  ::ROOT::Double_v v5 = v4 * v;
107  ::ROOT::Double_v v6 = v5 * v;
108  ::ROOT::Double_v v7 = v6 * v;
109  ::ROOT::Double_v v8 = v7 * v;
110 
111  vecCore::MaskedAssign<::ROOT::Double_v>(
112  result, mask1, v * (p10 + p11 * v2 + p12 * v4 + p13 * v6) / (q10 + q11 * v2 + q12 * v4 + v6));
113  vecCore::MaskedAssign<::ROOT::Double_v>(
114  result, mask2,
115  ::ROOT::Double_v(1.0) -
116  (p20 + p21 * v + p22 * v2 + p23 * v3 + p24 * v4 + p25 * v5 + p26 * v6 + p27 * v7) /
117  (vecCore::math::Exp(v2) * (q20 + q21 * v + q22 * v2 + q23 * v3 + q24 * v4 + q25 * v5 + q26 * v6 + v7)));
118  vecCore::MaskedAssign<::ROOT::Double_v>(result, mask3,
119  ::ROOT::Double_v(1.0) -
120  (c1 + (p30 * v8 + p31 * v6 + p32 * v4 + p33 * v2 + p34) /
121  ((q30 * v8 + q31 * v6 + q32 * v4 + q33 * v2 + q34) * v2)) /
122  (v * vecCore::math::Exp(v2)));
123 
124  return vecCore::Blend<::ROOT::Double_v>(x > 0, ::ROOT::Double_v(0.5) + ::ROOT::Double_v(0.5) * result,
125  ::ROOT::Double_v(0.5) * (::ROOT::Double_v(1) - result));
126 }
127 
128 ////////////////////////////////////////////////////////////////////////////////
129 /// Vectorized implementation of Bessel function I_0(x) for a vector x.
130 ::ROOT::Double_v BesselI0_Split_More(::ROOT::Double_v &ax)
131 {
132  ::ROOT::Double_v y = 3.75 / ax;
133  return (vecCore::math::Exp(ax) / vecCore::math::Sqrt(ax)) *
134  (0.39894228 +
135  y * (1.328592e-2 +
136  y * (2.25319e-3 +
137  y * (-1.57565e-3 +
138  y * (9.16281e-3 +
139  y * (-2.057706e-2 + y * (2.635537e-2 + y * (-1.647633e-2 + y * 3.92377e-3))))))));
140 }
141 
142 ::ROOT::Double_v BesselI0_Split_Less(::ROOT::Double_v &x)
143 {
144  ::ROOT::Double_v y = x * x * 0.071111111;
145 
146  return 1.0 +
147  y * (3.5156229 + y * (3.0899424 + y * (1.2067492 + y * (0.2659732 + y * (3.60768e-2 + y * 4.5813e-3)))));
148 }
149 
151 {
153 
154  return vecCore::Blend<::ROOT::Double_v>(ax <= 3.75, BesselI0_Split_Less(x), BesselI0_Split_More(ax));
155 }
156 
157 ////////////////////////////////////////////////////////////////////////////////
158 /// Vectorized implementation of modified Bessel function I_1(x) for a vector x.
159 ::ROOT::Double_v BesselI1_Split_More(::ROOT::Double_v &ax, ::ROOT::Double_v &x)
160 {
161  ::ROOT::Double_v y = 3.75 / ax;
162  ::ROOT::Double_v result =
164  (0.39894228 + y * (-3.988024e-2 +
165  y * (-3.62018e-3 +
166  y * (1.63801e-3 + y * (-1.031555e-2 +
167  y * (2.282967e-2 + y * (-2.895312e-2 +
168  y * (1.787654e-2 + y * -4.20059e-3))))))));
169  return vecCore::Blend<::ROOT::Double_v>(x < 0, ::ROOT::Double_v(-1.0) * result, result);
170 }
171 
172 ::ROOT::Double_v BesselI1_Split_Less(::ROOT::Double_v &x)
173 {
174  ::ROOT::Double_v y = x * x * 0.071111111;
175 
176  return x * (0.5 + y * (0.87890594 +
177  y * (0.51498869 + y * (0.15084934 + y * (2.658733e-2 + y * (3.01532e-3 + y * 3.2411e-4))))));
178 }
179 
181 {
183 
184  return vecCore::Blend<::ROOT::Double_v>(ax <= 3.75, BesselI1_Split_Less(x), BesselI1_Split_More(ax, x));
185 }
186 
187 ////////////////////////////////////////////////////////////////////////////////
188 /// Vectorized implementation of Bessel function J0(x) for a vector x.
189 ::ROOT::Double_v BesselJ0_Split1_More(::ROOT::Double_v &ax)
190 {
191  ::ROOT::Double_v z = 8 / ax;
192  ::ROOT::Double_v y = z * z;
193  ::ROOT::Double_v xx = ax - 0.785398164;
194  ::ROOT::Double_v result1 =
195  1 + y * (-0.1098628627e-2 + y * (0.2734510407e-4 + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
196  ::ROOT::Double_v result2 =
197  -0.1562499995e-1 + y * (0.1430488765e-3 + y * (-0.6911147651e-5 + y * (0.7621095161e-6 - y * 0.934935152e-7)));
198  return vecCore::math::Sqrt(0.636619772 / ax) *
199  (vecCore::math::Cos(xx) * result1 - z * vecCore::math::Sin(xx) * result2);
200 }
201 
202 ::ROOT::Double_v BesselJ0_Split1_Less(::ROOT::Double_v &x)
203 {
204  ::ROOT::Double_v y = x * x;
205  return (57568490574.0 +
206  y * (-13362590354.0 + y * (651619640.7 + y * (-11214424.18 + y * (77392.33017 + y * -184.9052456))))) /
207  (57568490411.0 + y * (1029532985.0 + y * (9494680.718 + y * (59272.64853 + y * (267.8532712 + y)))));
208 }
209 
211 {
213  return vecCore::Blend<::ROOT::Double_v>(ax < 8, BesselJ0_Split1_Less(x), BesselJ0_Split1_More(ax));
214 }
215 
216 ////////////////////////////////////////////////////////////////////////////////
217 /// Vectorized implementation of Bessel function J1(x) for a vector x.
218 ::ROOT::Double_v BesselJ1_Split1_More(::ROOT::Double_v &ax, ::ROOT::Double_v &x)
219 {
220  ::ROOT::Double_v z = 8 / ax;
221  ::ROOT::Double_v y = z * z;
222  ::ROOT::Double_v xx = ax - 2.356194491;
223  ::ROOT::Double_v result1 =
224  1 + y * (0.183105e-2 + y * (-0.3516396496e-4 + y * (0.2457520174e-5 + y * -0.240337019e-6)));
225  ::ROOT::Double_v result2 =
226  0.04687499995 + y * (-0.2002690873e-3 + y * (0.8449199096e-5 + y * (-0.88228987e-6 - y * 0.105787412e-6)));
227  ::ROOT::Double_v result =
228  vecCore::math::Sqrt(0.636619772 / ax) * (vecCore::math::Cos(xx) * result1 - z * vecCore::math::Sin(xx) * result2);
229  vecCore::MaskedAssign<::ROOT::Double_v>(result, x < 0, -result);
230  return result;
231 }
232 
233 ::ROOT::Double_v BesselJ1_Split1_Less(::ROOT::Double_v &x)
234 {
235  ::ROOT::Double_v y = x * x;
236  return x *
237  (72362614232.0 +
238  y * (-7895059235.0 + y * (242396853.1 + y * (-2972611.439 + y * (15704.48260 + y * -30.16036606))))) /
239  (144725228442.0 + y * (2300535178.0 + y * (18583304.74 + y * (99447.43394 + y * (376.9991397 + y)))));
240 }
241 
243 {
245  return vecCore::Blend<::ROOT::Double_v>(ax < 8, BesselJ1_Split1_Less(x), BesselJ1_Split1_More(ax, x));
246 }
247 
248 } // namespace TMath
249 
250 #endif
ROOT::Math::Cephes::gamma
double gamma(double x)
Definition: SpecFuncCephes.cxx:339
TMath::BesselJ1
Double_t BesselJ1(Double_t x)
Bessel function J0(x) for any real x.
Definition: TMath.cxx:1644
ROOT::v5::TFastFun::Sqrt
Double_t Sqrt(Double_t x)
Definition: TFormulaPrimitive_v5.cxx:289
TMath::BesselI0
Double_t BesselI0(Double_t x)
integer order modified Bessel function K_n(x)
Definition: TMath.cxx:1401
TMath::Cos
Double_t Cos(Double_t)
Definition: TMath.h:643
TMath::BesselJ0
Double_t BesselJ0(Double_t x)
modified Bessel function K_1(x)
Definition: TMath.cxx:1609
TMath::Gaus
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
TMath::Exp
Double_t Exp(Double_t x)
Definition: TMath.h:727
TMath::LaplaceDist
Double_t LaplaceDist(Double_t x, Double_t alpha=0, Double_t beta=1)
Computes the probability density function of Laplace distribution at point x, with location parameter...
Definition: TMath.cxx:2325
x
Double_t x[n]
Definition: legend1.C:17
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
Bool_t
bool Bool_t
Definition: RtypesCore.h:63
ROOT::Math::beta
double beta(double x, double y)
Calculates the beta function.
Definition: SpecFuncMathCore.cxx:111
VectorizedTMath.h
v
@ v
Definition: rootcling_impl.cxx:3635
TMath::LaplaceDistI
Double_t LaplaceDistI(Double_t x, Double_t alpha=0, Double_t beta=1)
Computes the distribution function of Laplace distribution at point x, with location parameter alpha ...
Definition: TMath.cxx:2341
TMath::Log2
Double_t Log2(Double_t x)
Definition: TMath.cxx:101
h
#define h(i)
Definition: RSha256.hxx:106
TMath::BreitWigner
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
TMath::Sin
Double_t Sin(Double_t)
Definition: TMath.h:639
y
Double_t y[n]
Definition: legend1.C:17
ROOT::Double_v
Double_t Double_v
Definition: Types.h:51
TMath::BesselI1
Double_t BesselI1(Double_t x)
modified Bessel function K_0(x)
Definition: TMath.cxx:1469
TMath::Freq
Double_t Freq(Double_t x)
Computation of the normal frequency function freq(x).
Definition: TMath.cxx:265
sigma
const Double_t sigma
Definition: h1analysisProxy.h:11
Double_t
double Double_t
Definition: RtypesCore.h:59
v3
@ v3
Definition: rootcling_impl.cxx:3639
v2
@ v2
Definition: rootcling_impl.cxx:3638
v4
@ v4
Definition: rootcling_impl.cxx:3640
TMath
TMath.
Definition: TMathBase.h:35
c1
return c1
Definition: legend1.C:41