MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExtendedMath.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include "NumericalVector.h"
27 #include "ExtendedMath.h"
28 #include "Quaternion.h"
29 #include <cmath>
30 #include "Logger.h"
32 //#include <limits>
33 //#include <iostream>
34 //#include <sys/stat.h>
35 #include <iomanip>
36 //#include <cmath>
37 //#include <sstream>
38 //#include <cstdlib>
39 //#include <limits>
40 //#include <string>
41 
42 
43 // sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...
45 {
46  Mdouble N = floor(x / (2.0 * constants::pi) + 0.5);
47  x -= N * 2.0 * constants::pi;
48  Mdouble Sum = 0;
49  Mdouble Power = x;
50  Mdouble Sign = 1;
51  const Mdouble x2 = x * x;
52  Mdouble Fact = 1.0;
53  for (unsigned int i = 1; i < 25; i += 2)
54  {
55  Sum += Sign * Power / Fact;
56  Power *= x2;
57  Fact *= (i + 1) * (i + 2);
58  Sign *= -1.0;
59  }
60  return Sum;
61 }
62 
63 // cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + ...
65 {
66  Mdouble N = floor(x / (2.0 * constants::pi) + 0.5);
67  x -= N * 2.0 * constants::pi;
68  Mdouble Sum = 1.0;
69  Mdouble Power = 1;
70  Mdouble Sign = 1;
71  const Mdouble x2 = x * x;
72  Mdouble Fact = 1.0;
73  for (unsigned int i = 2; i < 25; i += 2)
74  {
75  Power *= x2;
76  Fact *= i * (i - 1);
77  Sign *= -1.0;
78  Sum += Sign * Power / Fact;
79  }
80  return Sum;
81 }
82 
83 //from: http://www.codeproject.com/Tips/311714/Natural-Logarithms-and-Exponent
85 {
86  Mdouble X, P, Frac, I, L;
87  X = Exponent;
88  Frac = X;
89  P = (1.0 + X);
90  I = 1.0;
91 
92  do
93  {
94  I++;
95  Frac *= (X / I);
96  L = P;
97  P += Frac;
98  } while (L != P);
99 
100  return P;
101 }
102 
105 {
106  Mdouble N, P, L, R, A, E;
107  E = 2.71828182845905;
108  P = Power;
109  N = 0.0;
110 
111  // This speeds up the convergence by calculating the integral
112  while (P >= E)
113  {
114  P /= E;
115  N++;
116  }
117  N += (P / E);
118  P = Power;
119  do
120  {
121  A = N;
122  L = (P / (exp(N - 1.0)));
123  R = ((N - 1.0) * E);
124  N = ((L + R) / E);
125  } while (N < A);
126  //} while (N != A);
127 
128  return N;
129 }
130 
138 {
139  const Mdouble ep = 1e-5;
140 
141  if (gamma_in > 1.0 + ep)
142  {
143  return ((gamma_in - 1) * gamma(gamma_in - 1));
144  }
145  else
146  {
147 
148  if ((gamma_in - ep < 1.0) && (gamma_in + ep > 1.0))
149  return 1.0;
150  else if ((gamma_in - ep < 0.5) && (gamma_in + ep > 0.5))
151  return constants::sqrt_pi;
152  else
153  return std::numeric_limits<Mdouble>::quiet_NaN();
154  }
155 } //end func gamma
156 
165 {
166  return std::exp(std::lgamma(z) + std::lgamma(w) - std::lgamma(z + w));
167 }
168 
172 Mdouble mathsFunc::chi_squared(const Mdouble x, const unsigned int k)
173 {
174 
175  Mdouble prefactor = pow(2, k / 2.0) * gamma(k / 2.0);
176  Mdouble mainfactor = pow(x, k / 2.0 - 1) * exp(x / -2.0);
177 
178  return mainfactor / prefactor;
179 
180 }
181 
188 Mdouble mathsFunc::chi_squared_prob(const Mdouble x_max, const unsigned int k)
189 {
190 
191 //The current value was picked by tried were it stopped effect the 4 d.p.
192  const int num_steps_per_unit = 100;
193  Mdouble sum = 0;
194  Mdouble x = 0;
195  long int num_steps = static_cast<int>(num_steps_per_unit * x_max);
196 //Use trapezional rule, but ignoring the ends
197  for (int i = 0; i < num_steps; i++)
198  {
199  x = x_max / num_steps * (i + 0.5);
200  sum = sum + chi_squared(x, k);
201  }
202  return 1.0 - sum * x_max / num_steps;
203 
204 }
205 
207  Mdouble endCondition, Mdouble curVal)
208 {
209  if (std::abs(max - min) < endCondition)
210  {
211  return 0.5 * (min + max);
212  }
213  std::cout << "Min=" << min << " Max=" << max << " diff=" << max - min << std::endl;
214  Mdouble resphi = 2 - 0.5 * (1 + std::sqrt(5));
215  Mdouble x;
216  if (max - cur > cur - min)
217  {
218  x = cur + resphi * (max - cur);
219  }
220  else
221  {
222  x = cur - resphi * (cur - min);
223  }
224  if (std::isnan(curVal))
225  curVal = function(cur);
226  Mdouble xVal = function(x);
227  if (xVal < curVal)
228  {
229  if (max - cur > cur - min)
230  {
231  return goldenSectionSearch(function, cur, x, max, endCondition, xVal);
232  }
233  else
234  {
235  return goldenSectionSearch(function, min, x, cur, endCondition, xVal);
236  }
237  }
238  else
239  {
240  if (max - cur > cur - min)
241  {
242  return goldenSectionSearch(function, min, cur, x, endCondition, curVal);
243  }
244  else
245  {
246  return goldenSectionSearch(function, x, cur, max, endCondition, curVal);
247  }
248  }
249 }
250 
252 {
253  return std::abs(v1 - v2) <= absError;
254 }
255 
256 
257 bool mathsFunc::isEqual(Vec3D v1, Vec3D v2, Mdouble absError)
258 {
259  return isEqual(v1.X, v2.X, absError) && isEqual(v1.Y, v2.Y, absError) && isEqual(v1.Z, v2.Z, absError);
260 }
261 
263 {
264  return (isEqual(m1.XX, m2.XX, absError)
265  && isEqual(m1.XY, m2.XY, absError)
266  && isEqual(m1.XZ, m2.XZ, absError)
267  && isEqual(m1.YX, m2.YX, absError)
268  && isEqual(m1.YY, m2.YY, absError)
269  && isEqual(m1.YZ, m2.YZ, absError)
270  && isEqual(m1.ZX, m2.ZX, absError)
271  && isEqual(m1.ZY, m2.ZY, absError)
272  && isEqual(m1.ZZ, m2.ZZ, absError));
273 }
274 
276 {
277  return (isEqual(m1.XX, m2.XX, absError)
278  && isEqual(m1.XY, m2.XY, absError)
279  && isEqual(m1.XZ, m2.XZ, absError)
280  && isEqual(m1.YY, m2.YY, absError)
281  && isEqual(m1.YZ, m2.YZ, absError)
282  && isEqual(m1.ZZ, m2.ZZ, absError));
283 }
284 
285 bool mathsFunc::isEqual(Quaternion v1, Quaternion v2, double absError)
286 {
287  return isEqual(v1.getComponent(0), v2.getComponent(0), absError) &&
288  isEqual(v1.getComponent(1), v2.getComponent(1), absError) &&
289  isEqual(v1.getComponent(2), v2.getComponent(2), absError) &&
290  isEqual(v1.getComponent(3), v2.getComponent(3), absError);
291 }
292 
294 {
295  const Mdouble* p = coef;
296  Mdouble b0 = *p++;
297  Mdouble b1 = 0, b2;
298  int i = N - 1;
299 
300  logger.assert(i > 0, "i is greater than 0");
301  do
302  {
303  b2 = b1;
304  b1 = b0;
305  b0 = x * b1 - b2 + *p++;
306  } while (--i);
307 
308  return (0.5 * (b0 - b2));
309 }
310 
312 {
313  // Coefficients for [0..8]
314  const Mdouble A[] =
315  {
316  -4.415341646479339379501E-18,
317  3.330794518822238097831E-17,
318  -2.431279846547954693591E-16,
319  1.715391285555133030611E-15,
320  -1.168533287799345168081E-14,
321  7.676185498604935616881E-14,
322  -4.856446783111929460901E-13,
323  2.955052663129639834611E-12,
324  -1.726826291441555707231E-11,
325  9.675809035373236912241E-11,
326  -5.189795601635262906661E-10,
327  2.659823724682386650351E-9,
328  -1.300025009986248042121E-8,
329  6.046995022541918949321E-8,
330  -2.670793853940611733911E-7,
331  1.117387539120103718151E-6,
332  -4.416738358458750563591E-6,
333  1.644844807072889708931E-5,
334  -5.754195010082103703981E-5,
335  1.885028850958416557291E-4,
336  -5.763755745385823658851E-4,
337  1.639475616941335798421E-3,
338  -4.324309995050575944301E-3,
339  1.054646039459499831831E-2,
340  -2.373741480589946881561E-2,
341  4.930528423967070848781E-2,
342  -9.490109704804764442101E-2,
343  1.716209015222087753491E-1,
344  -3.046826723431983986831E-1,
345  6.767952744094760849951E-1
346  };
347 
348  // Coefficients for [8..infinity]
349  const Mdouble B[] =
350  {
351  -7.233180487874753954561E-18,
352  -4.830504485944182071261E-18,
353  4.465621420296759999011E-17,
354  3.461222867697461093101E-17,
355  -2.827623980516583484941E-16,
356  -3.425485619677219134621E-16,
357  1.772560133056526383601E-15,
358  3.811680669352622420751E-15,
359  -9.554846698828307648701E-15,
360  -4.150569347287222086631E-14,
361  1.540086217521409826911E-14,
362  3.852778382742142701141E-13,
363  7.180124451383666233671E-13,
364  -1.794178531506806117781E-12,
365  -1.321581184044771311881E-11,
366  -3.149916527963241364541E-11,
367  1.188914710784643834241E-11,
368  4.940602388224969589101E-10,
369  3.396232025708386345151E-9,
370  2.266668990498178064591E-8,
371  2.048918589469063741831E-7,
372  2.891370520834756482971E-6,
373  6.889758346916823984261E-5,
374  3.369116478255694089901E-3,
375  8.044904110141088316081E-1
376  };
377 
378  if (x < 0)
379  x = -x;
380 
381  if (x <= 8.0)
382  {
383  Mdouble y = (x / 2.0) - 2.0;
384  return (chebyshev(y, A, 30));
385  }
386 
387  return (chebyshev(32.0 / x - 2.0, B, 25) / sqrt(x));
388 
389 }
390 
392 {
393  if (x < 0)
394  x = -x;
395  return exp(x) * I0_exp(x);
396 }
397 
398 // This function computes all n order and positive order m associated Legendre polynomials using recursive formulations
399 // The polynomials are evaluated at x.
400 // They are organised as follows: P_0^0, P_1^0, P_1^1, P_2^0, P_2^1, ...
402 {
403  //Given n and m, we only have to compute P_n^(|m|)
404  //The function will return all these P values for theta
405  std::size_t nTerms = 0.5 * (n + 1) * (n + 2);
406  NumericalVector<> polynomials(nTerms);
407 
408  size_t location_current;
409  size_t location_previous;
410  Mdouble temp;
411 
412  polynomials(0) = 1; //P_0^0 = 1;
413  for (int l = 1; l <= n; l++)
414  {
415  //first compute P_l^l
416  location_current = 0.5 * (l + 1) * (l + 2) - 1;
417  location_previous = location_current - (l + 1);
418  polynomials(location_current) = -(2.0 * (l - 1.0) + 1.0) * std::sqrt(1.0 - x * x) *
419  polynomials(location_previous); // Recursive formula from wiki
420 
421  //second, compute P_l^(l-1) based on P_(l-1)^(l-1)
422  polynomials(location_current - 1) =
423  x * (2.0 * (l - 1) + 1) * polynomials(location_previous); // Recursive formula from wiki
424  }
425 
426  //thirdly, compute other values
427  for (int l = 2; l <= n; l++)
428  {
429  for (int m = (l - 2); m >= 0; m--)
430  {
431  location_current = (0.5 * (l + 1) * (l + 2) - 1) - l + m;
432  temp = polynomials(location_current + 2) +
433  2.0 * (m + 1) * x / sqrt(1 - x * x) * polynomials(location_current + 1);
434  polynomials(location_current) = temp / ((m - l) * (l + m + 1)); // variation on greengard eqn 3.34 from wiki
435  }
436  }
437 
438 
439  return polynomials;
440 }
441 
442 
443 // This function computes up to p order spherical harmonics as function of theta and phi
444 // They are organised as follows: Y_0^0, Y_1^-1, Y_1^0, Y_1^1, Y_2^-2, ...
446 {
447  std::size_t nTerms = 0.5 * (p + 1) * (2 * p + 2);
450 
451  //Compute the spherical harmonics
452  for (int n = 0; n <= p; n++)
453  {
454  for (int mt = -n; mt <= n; mt++)
455  {
456  Mdouble m = mt;
457  Mdouble m_abs = std::abs(mt);
458  std::size_t location_current = n * n + (m + n); //n^2 is begin of Y_n^-n
459  std::size_t location_polynomial = 0.5 * n * (n + 1) + m_abs;
460  int fact1 = mathsFunc::factorial(n - m_abs);
461  int fact2 = mathsFunc::factorial(n + m_abs);
462  Mdouble fact = 1.0 * fact1 / fact2;
463  std::complex<Mdouble> value =
464  std::sqrt(fact) * polynomials(location_polynomial) * std::exp(constants::i * m * phi);
465  Y(location_current) = value;
466  }
467  }
468  return Y;
469 }
470 
471 //Compute all squaredFactorials (see eqn 5.23 in a short course on fast multipole methods) up to order p
472 // They are organised as follows: A_0^0, A_1^-1, A_1^0, A_1^1, A_2^-2, ...
474 {
475  std::size_t nTerms = 0.5 * (p + 1) * (2 * p + 2);
476  NumericalVector<> squaredFactorials(nTerms);
477 
478  for (int n = 0; n <= p; n++)
479  {
480  for (int m = -n; m <= n; m++)
481  {
482  std::size_t location = n * n + (m + n); //(n^2 is begin of Y_n)
483  int fact1 = mathsFunc::factorial(n - m);
484  int fact2 = mathsFunc::factorial(n + m);
485  Mdouble fact = fact1 * fact2;
486  squaredFactorials(location) = std::pow(-1.0, n) / std::sqrt(fact);
487  }
488  }
489 
490  return squaredFactorials;
491 }
492 
493 
494 
Mdouble XY
Definition: Matrix.h:43
Implementation of a 3D quaternion (by Vitaliy).
Definition: Quaternion.h:62
Mdouble X
the vector components
Definition: Vector.h:65
Mdouble ZY
Definition: Matrix.h:43
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
Mdouble goldenSectionSearch(Mdouble(*function)(const Mdouble), Mdouble min, Mdouble cur, Mdouble max, Mdouble endCondition, Mdouble curVal=std::numeric_limits< Mdouble >::quiet_NaN())
This function performs a golden section search to find the location of the minimum of a function...
Mdouble getComponent(int index) const
Returns the requested component of this Quaternion.
Definition: Quaternion.cc:232
Mdouble XZ
Definition: Matrix.h:43
Mdouble exp(Mdouble Exponent)
Definition: ExtendedMath.cc:84
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
Mdouble ZX
Definition: Matrix.h:43
Mdouble I0(Mdouble x)
const Mdouble sqrt_pi
Definition: ExtendedMath.h:46
Mdouble log(Mdouble Power)
Mdouble ZZ
Definition: Matrix.h:43
bool isEqual(Mdouble v1, Mdouble v2, Mdouble absError)
Compares the difference of two Mdouble with an absolute error, useful in UnitTests.
Mdouble chi_squared_prob(Mdouble x, unsigned int k)
This is the function which actually gives the probability back using a chi squared test...
Mdouble YZ
Definition: Matrix.h:43
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
constexpr T factorial(const T t)
factorial function
Definition: ExtendedMath.h:153
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
Mdouble chi_squared(Mdouble x, unsigned int k)
This is a chi_squared function return the value x and degrees of freedom k.
const Mdouble pi
Definition: ExtendedMath.h:45
Mdouble I0_exp(Mdouble x)
Mdouble YX
Definition: Matrix.h:43
Mdouble chebyshev(Mdouble x, const Mdouble coef[], int N)
Namespace for evaluating the zeroth modified Bessel function of the first kind, I0(x), required in StatisticsPoint.hcc.
NumericalVector< std::complex< Mdouble > > sphericalHarmonics(int p, Mdouble theta, Mdouble phi)
Mdouble Y
Definition: Vector.h:65
Mdouble XX
all nine matrix elements
Definition: Matrix.h:43
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
NumericalVector associatedLegendrePolynomials(int n, Mdouble x)
Mdouble YY
Definition: Matrix.h:43
Implementation of a 3D matrix.
Definition: Matrix.h:37
Definition: Vector.h:49
Mdouble Z
Definition: Vector.h:65
Implementation of a 3D symmetric matrix.
NumericalVector computeSquaredFactorialValues(int p)
Mdouble XX
The six distinctive matrix elements.