LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
BinomFact_omp.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 #include <cassert>
26 #include <cmath>
27 #include <vector>
28 
30 
31 namespace lsst {
32 namespace meas {
33 namespace algorithms {
34 namespace shapelet {
35 
36  double fact(int i)
37  {
38  // return i!
39  assert(i>=0);
40  static std::vector<double> f(10);
41  static bool first=true;
42 
43  if (first) {
44  f[0] = f[1] = 1.;
45  for(int j=2;j<10;j++) f[j] = f[j-1]*(double)j;
46  first = false;
47  }
48  if (i>=(int)f.size()) {
49 #ifdef _OPENMP
50 #pragma omp critical (fact)
51 #endif
52  {
53  for(int j=f.size();j<=i;j++)
54  f.push_back(f[j-1]*(double)j);
55  }
56  assert(i==(int)f.size()-1);
57  }
58  assert(i<(int)f.size());
59  return f[i];
60  }
61 
62  double sqrtfact(int i)
63  {
64  // return sqrt(i!)
65  static std::vector<double> f(10);
66  static bool first=true;
67  if (first) {
68  f[0] = f[1] = 1.;
69  for(int j=2;j<10;j++) f[j] = f[j-1]*sqrt((double)j);
70  first = false;
71  }
72  if (i>=(int)f.size()) {
73 #ifdef _OPENMP
74 #pragma omp critical (sqrtfact)
75 #endif
76  {
77  for(int j=f.size();j<=i;j++)
78  f.push_back(f[j-1]*sqrt((double)j));
79  }
80  }
81  assert(i<(int)f.size());
82  return f[i];
83  }
84 
85  double binom(int i,int j)
86  {
87  // return iCj, i!/(j!(i-j)!)
88  static std::vector<std::vector<double> > f(10);
89  static bool first=true;
90  if (first) {
91  f[0] = std::vector<double>(1,1.);
92  f[1] = std::vector<double>(2,1.);
93  for(int i1=2;i1<10;i1++) {
94  f[i1] = std::vector<double>(i1+1);
95  f[i1][0] = f[i1][i1] = 1.;
96  for(int j1=1;j1<i1;j1++) f[i1][j1] = f[i1-1][j1-1] + f[i1-1][j1];
97  }
98  first = false;
99  }
100  if (i>=(int)f.size()) {
101 #ifdef _OPENMP
102 #pragma omp critical (binom)
103 #endif
104  {
105  for(int i1=f.size();i1<=i;i1++) {
106  f.push_back(std::vector<double>(i1+1,1.));
107  for(int j1=1;j1<i1;j1++) f[i1][j1] = f[i1-1][j1-1] + f[i1-1][j1];
108  }
109  }
110  assert(i==(int)f.size()-1);
111  }
112  if (j<0 || j>i) return 0.;
113  assert(i<(int)f.size());
114  assert(j<(int)f[i].size());
115  return f[i][j];
116  }
117 
118  double sqrtn(int i)
119  {
120  // return sqrt(i)
121  static std::vector<double> f(10);
122  static bool first=true;
123  if (first) {
124  for(int j=0;j<10;j++) f[j] = std::sqrt((double)j);
125  first = false;
126  }
127  if (i>=(int)f.size()) {
128 #ifdef _OPENMP
129 #pragma omp critical (sqrtn)
130 #endif
131  {
132  for(int j=f.size();j<=i;j++)
133  f.push_back(std::sqrt((double)j));
134  }
135  }
136  assert(i<(int)f.size());
137  return f[i];
138  }
139 
140 
141 
142 }}}}