lbeta.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LBETA_HPP
#define STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LBETA_HPP
#ifdef STAN_OPENCL
 
#include <stan/math/opencl/stringify.hpp>
#include <string>
 
namespace stan {
namespace math {
namespace opencl_kernels {
 
// \cond
static constexpr const char* lbeta_device_function
    = "\n"
      "#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LBETA\n"
      "#define STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LBETA\n" STRINGIFY(
          // \endcond
          double stan_lbeta(double a, double b) {
            if (isnan(a) || isnan(b)) {
              return NAN;
            }
 
            double x;  // x is the smaller of the two
            double y;
            if (a < b) {
              x = a;
              y = b;
            } else {
              x = b;
              y = a;
            }
 
            // Special cases
            if (x == 0) {
              return INFINITY;
            }
            if (isinf(y)) {
              return -INFINITY;
            }
 
            // For large x or y, separate the lgamma values into Stirling
            // approximations and appropriate corrections. The Stirling
            // approximations allow for analytic simplification and the
            // corrections are added later.
            //
            // The overall approach is inspired by the code in R, where the
            // algorithm is credited to W. Fullerton of Los Alamos Scientific
            // Laboratory
            if (y < LGAMMA_STIRLING_DIFF_USEFUL) {
              // both small
              return lgamma(x) + lgamma(y) - lgamma(x + y);
            }
            double x_over_xy = x / (x + y);
            double log_xpy = log(x + y);
            if (x < LGAMMA_STIRLING_DIFF_USEFUL) {
              // y large, x small
              double stirling_diff
                  = lgamma_stirling_diff(y) - lgamma_stirling_diff(x + y);
              double stirling
                  = (y - 0.5) * log1p(-x_over_xy) + x * (1 - log_xpy);
              return stirling + lgamma(x) + stirling_diff;
            }
 
            // both large
            double stirling_diff = lgamma_stirling_diff(x)
                                   + lgamma_stirling_diff(y)
                                   - lgamma_stirling_diff(x + y);
            double stirling = (x - 0.5) * (log(x) - log_xpy)
                              + y * log1p(-x_over_xy)
                              + 0.5 * (M_LN2 + log(M_PI)) - 0.5 * log(y);
            return stirling + stirling_diff;
          }
          // \cond
          ) "\n#endif\n";  // NOLINT
// \endcond
 
}  // namespace opencl_kernels
}  // namespace math
}  // namespace stan
 
#endif
#endif