1 /**
2  *
3  *   Copyright (c) 2005-2021 by Pierre-Henri WUILLEMIN(_at_LIP6) & Christophe GONZALES(_at_AMU)
4  *   info_at_agrum_dot_org
5  *
6  *  This library is free software: you can redistribute it and/or modify
7  *  it under the terms of the GNU Lesser General Public License as published by
8  *  the Free Software Foundation, either version 3 of the License, or
9  *  (at your option) any later version.
10  *
11  *  This library is distributed in the hope that it will be useful,
12  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  *  GNU Lesser General Public License for more details.
15  *
16  *  You should have received a copy of the GNU Lesser General Public License
17  *  along with this library.  If not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 
22 /**
23  * @file
24  * @brief The class for computing Log2 (Gamma(x)).
25  *
26  * @author Christophe GONZALES(_at_AMU) and Pierre-Henri WUILLEMIN(_at_LIP6)
27  */
28 
29 namespace gum {
30 
gammaLog2(double x)31   ALWAYS_INLINE double GammaLog2::gammaLog2(double x) const {
32     if (x <= 0) GUM_ERROR(OutOfBounds, "log2(gamma()) should be called with a positive argument")
33 
34     // if x is small, use precomputed values
35     if (x < 50) {
36       if (x >= 0.01) {
37         if (_requires_precision_) {
38           const Idx index = int(x * 100);
39           return _small_values_[index]
40                + (_small_values_[index + 1] - _small_values_[index]) * double(x * 100 - index);
41         } else {
42           const Idx index = int(x * 100 + 0.5);
43           return _small_values_[index];
44         }
45       } else {
46         // for very small values of x, Gamma(x) is approximately equal to
47         // 1/x. Hence gammaLog2(x) is approximately equal to log2(1/x)
48         return std::log2(1.0 / x);
49       }
50     }
51 
52     // returns the approximation by the stirling formula
53     return (_log_sqrt_2pi_ + (x - 0.5f) * log(x) - x + log(1.0 + 1.0 / (12 * x))) * _inv_log2_;
54   }
55 
operator()56   ALWAYS_INLINE double GammaLog2::operator()(double x) const { return gammaLog2(x); }
57 
setPrecision(bool prec)58   INLINE void GammaLog2::setPrecision(bool prec) { _requires_precision_ = prec; }
59 
60 } /* namespace gum */
61