Chombo + EB + MF  3.2
Factorial.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _FACTORIAL_H_
12 #define _FACTORIAL_H_
13 
14 #include "MayDay.H"
15 #include "RealVect.H"
16 #include "IndexTM.H"
17 #include "NamespaceHeader.H"
18 
19 /**
20  * These helper functions calculate factorials, binomial coeff, etc.
21  * for multinomials.
22  */
23 
24 /// Calculates factorial for an integer
25 inline Real factorial(const int n)
26 {
27  CH_assert(n >= 0);
28  Real nfact = 1.0;
29  for (int i = 2; i <= n; nfact*= i, i++);
30  return nfact;
31 }
32 /// computes x^p
33 inline Real POW(const Real& a_x, const int& a_p)
34 {
35  Real retval = 1;
36  for(int iexp = 0; iexp < Abs(a_p); iexp++)
37  {
38  if(a_p >= 0)
39  {
40  retval *= a_x;
41  }
42  else
43  {
44  retval /= a_x;
45  }
46  }
47  return retval;
48 }
49 
50 /// Calculates the binomial coefficient, "n choose k"
51 inline Real nCk(const int n, const int k)
52 {
53  CH_assert((n >= k)&&(k >= 0)&&(n >= 0));
54  Real nck = 1.0;
55  for (int i = n-k+1; i <= n; nck*= i, i++);
56  return nck / factorial(k);
57 }
58 
59 /// Calculates factorials for a multinomial
60 template <int Dim>
62 {
63  Real pfact = 1;
64  for(int idir = 0; idir < Dim; idir++)
65  {
66  CH_assert(p[idir] >= 0);
67  pfact *= factorial(p[idir]);
68  }
69  return pfact;
70 }
71 
72 /// Calculates the multinomial coefficient, "p choose k"
73 template <int Dim>
74 inline Real pCk(const IndexTM<int, Dim>& p, const IndexTM<int, Dim>& k)
75 {
76  //CH_assert((p >= IntVect::Zero)&&(k >= IntVect::Zero)&&(p >= k));
77 
78  Real pfact = 1;
79  for(int idir = 0; idir < Dim; idir++)
80  {
81  pfact *= nCk(p[idir],k[idir]);
82  }
83  return pfact;
84 }
85 
86 
87 /// calculate x^p
88 template <int Dim>
89 inline Real power(const IndexTM<Real, Dim>& a_x, const IndexTM<int, Dim>& a_p)
90 {
91  Real retval = 1;
92  for(int idir = 0; idir < Dim; idir++)
93  {
94  if((a_p[idir] >= 0))
95  {
96  for(int ipow = 0; ipow < a_p[idir]; ipow++)
97  {
98  retval *= a_x[idir];
99  }
100  }
101  else
102  {
103  for(int ipow = 0; ipow < a_p[idir]; ipow++)
104  {
105  retval /= a_x[idir];
106  }
107  }
108  }
109  return retval;
110 }
111 
112 
113 #include "NamespaceFooter.H"
114 #endif
#define CH_assert(cond)
Definition: CHArray.H:37
Real pfactorial(const IndexTM< int, Dim > p)
Calculates factorials for a multinomial.
Definition: Factorial.H:61
Real POW(const Real &a_x, const int &a_p)
computes x^p
Definition: Factorial.H:33
Real nCk(const int n, const int k)
Calculates the binomial coefficient, "n choose k".
Definition: Factorial.H:51
Definition: IndexTM.H:36
Real pCk(const IndexTM< int, Dim > &p, const IndexTM< int, Dim > &k)
Calculates the multinomial coefficient, "p choose k".
Definition: Factorial.H:74
double Real
Definition: REAL.H:33
T Abs(const T &a_a)
Definition: Misc.H:53
Real power(const IndexTM< Real, Dim > &a_x, const IndexTM< int, Dim > &a_p)
calculate x^p
Definition: Factorial.H:89
Real factorial(const int n)
Calculates factorial for an integer.
Definition: Factorial.H:25
constexpr int ipow(int M)
Definition: Misc.H:80