lcd_even_binf_closed_form.h

Go to the documentation of this file.

#ifndef LCD_EVEN_BINF_CLOSED_FORM_H
#define LCD_EVEN_BINF_CLOSED_FORM_H
 
#include <cassert>
#include <cmath>
#include <cstddef>
 
#include "lcd_even_closed_form.h"  // lcd_ei()
 
constexpr double lcdBinfEulerGamma = 0.5772156649015328606;
 
inline double lcd_binf_b0_zero(size_t d, double e0, double ei0, double c) {
  assert(c > 0.00);
 
  if (d == 0) return 0.25 * e0 + 0.125 * c * ei0;
  if (d == 1) return -0.25 * ei0;
 
  double term = 0.50 / c;  // m = d
  double sum = term;
  for (size_t m = d; m > 2; --m) {
    term *= 2.00 * static_cast<double>(m - 2) / c;  // step down to m - 1
    sum += term;
  }
 
  return e0 * sum;
}
 
inline double lcd_binf_b0_zero_deriv(size_t d, double e0, double ei0,
                                     double c) {
  assert(c > 0.00);
 
  if (d == 0) return 0.125 * ei0;
  if (d == 1) return -0.25 * e0 / c;
 
  double term = 0.50 / c;      // m = d
  double sum = term;           // T_d(c)
  double sumDeriv = -term / c; // T_d'(c), factor -(d - m + 1) = -1 at m = d
  for (size_t m = d; m > 2; --m) {
    term *= 2.00 * static_cast<double>(m - 2) / c;  // step down to m - 1
    sum += term;
    sumDeriv -= static_cast<double>(d - m + 2) * term / c;
  }
 
  return e0 * (sumDeriv - 0.50 * sum);
}
 
inline double lcd_binf_bkk_zero(size_t k, double c) {
  assert(k >= 1);
  assert(c > 0.00);
 
  const double e0 = std::exp(-0.50 * c);
  const double ei0 = lcd_ei(-0.50 * c);
 
  double sum = 0.00;
  double binomial = 1.00;  // binom(k, j)
  for (size_t j = 0; j <= k; ++j) {
    const double term = lcd_binf_b0_zero(j, e0, ei0, c);
    sum += (j % 2 == 0) ? binomial * term : -binomial * term;
    binomial *= static_cast<double>(k - j) / static_cast<double>(j + 1);
  }
 
  return std::ldexp(sum, -static_cast<int>(k));
}
 
inline double lcd_binf_bkk_zero_deriv(size_t k, double c) {
  assert(k >= 1);
  assert(c > 0.00);
 
  const double e0 = std::exp(-0.50 * c);
  const double ei0 = lcd_ei(-0.50 * c);
 
  double sum = 0.00;
  double binomial = 1.00;  // binom(k, j)
  for (size_t j = 0; j <= k; ++j) {
    const double term = lcd_binf_b0_zero_deriv(j, e0, ei0, c);
    sum += (j % 2 == 0) ? binomial * term : -binomial * term;
    binomial *= static_cast<double>(k - j) / static_cast<double>(j + 1);
  }
 
  return std::ldexp(sum, -static_cast<int>(k));
}
 
inline double lcd_binf_s(size_t k, double c) {
  assert(k >= 1);
  assert(c > 0.00);
 
  double sum = 0.00;
  double binomial = static_cast<double>(k);  // binom(k, 1), stepped up below
  double coefficient = 0.50;                 // (j-2)! * 2^(j-3), j = 2
  double cPower = c;                         // c^(j-1),          j = 2
 
  for (size_t j = 2; j <= k; ++j) {
    binomial *= static_cast<double>(k - j + 1) / static_cast<double>(j);
    const double term = binomial * coefficient / cPower;
    sum += (j % 2 == 0) ? term : -term;
 
    coefficient *= 2.00 * static_cast<double>(j - 1);
    cPower *= c;
  }
 
  return sum;
}
 
inline double lcd_binf_s_deriv(size_t k, double c) {
  assert(k >= 1);
  assert(c > 0.00);
 
  double sum = 0.00;
  double binomial = static_cast<double>(k);  // binom(k, 1), stepped up below
  double coefficient = 0.50;                 // (j-2)! * 2^(j-3), j = 2
  double cPower = c;                         // c^(j-1),          j = 2
 
  for (size_t j = 2; j <= k; ++j) {
    binomial *= static_cast<double>(k - j + 1) / static_cast<double>(j);
    const double term = binomial * coefficient *
                        (1.00 - static_cast<double>(j)) / (cPower * c);
    sum += (j % 2 == 0) ? term : -term;
 
    coefficient *= 2.00 * static_cast<double>(j - 1);
    cPower *= c;
  }
 
  return sum;
}
 
inline double lcd_binf_phi(size_t k, double c) {
  assert(k >= 1);
  assert(c > 0.00);
 
  return 0.125 * (c + 2.00 * static_cast<double>(k)) * std::log(0.25 * c) +
         lcd_binf_s(k, c) -
         std::ldexp(lcd_binf_bkk_zero(k, c), static_cast<int>(k));
}
 
inline double lcd_binf_phi_deriv(size_t k, double c) {
  assert(k >= 1);
  assert(c > 0.00);
 
  return 0.125 * std::log(0.25 * c) +
         (c + 2.00 * static_cast<double>(k)) / (8.00 * c) +
         lcd_binf_s_deriv(k, c) -
         std::ldexp(lcd_binf_bkk_zero_deriv(k, c), static_cast<int>(k));
}
 
inline double lcd_binf_constant(size_t k) {
  assert(k >= 1);
 
  double harmonic = 0.00;  // sum_{i=2..k} 1 / (2i)
  for (size_t i = 2; i <= k; ++i) harmonic += 0.50 / static_cast<double>(i);
 
  return static_cast<double>(k) * harmonic - 0.50 -
         0.50 * static_cast<double>(k) * lcdBinfEulerGamma;
}
 
#endif  // LCD_EVEN_BINF_CLOSED_FORM_H