lcd_even_closed_form.h

Go to the documentation of this file.

#ifndef LCD_EVEN_CLOSED_FORM_H
#define LCD_EVEN_CLOSED_FORM_H
 
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_expint.h>
 
#include <cassert>
#include <cmath>
#include <cstddef>
 
inline double lcd_ei(double x) {
  assert(x < 0.00);
 
  gsl_sf_result result;
  const int status = gsl_sf_expint_Ei_e(x, &result);
 
  if (status == GSL_EUNDRFLW || status == GSL_EDOM || status == GSL_EOVRFLW) return 0.00;
  assert(status == GSL_SUCCESS);
 
  return result.val;
}
 
inline double lcd_a_n(size_t N, double b) {
  assert(N >= 1);
  assert(b >= 0.00);
 
  const double bSqrd = b * b;
 
  size_t n = (N % 2 == 0) ? 2 : 1;
  double aN = (N % 2 == 0)
                  ? 0.50 * (bSqrd - std::log1p(bSqrd))
                  : 0.50 * (b * std::sqrt(1.00 + bSqrd) - std::asinh(b));
 
  for (n += 2; n <= N; n += 2) {
    const double dn = static_cast<double>(n);
    const double tail =
        bSqrd * std::pow(bSqrd / (1.00 + bSqrd), (dn - 2.00) / 2.00);
    aN = (dn * aN - tail) / (dn - 2.00);
  }
 
  return aN;
}
 
inline double lcd_c_repulsion(double b, double c) {
  assert(b > 0.00);
  assert(c >= 0.00);
 
  const double bSqrd = b * b;
  if (c <= 0.00) return 0.50 * bSqrd;  // coincident samples
 
  const double z = -c / (4.00 * bSqrd);
  return 0.50 * bSqrd * std::exp(z) + 0.125 * c * lcd_ei(z);
}
 
struct LcdDeltaBIntegrandParams {
  size_t k;      
  size_t q;      
  double cValue; 
  bool reduced;  
};
 
inline double lcd_delta_b_integrand(double b, void* params) {
  const LcdDeltaBIntegrandParams* p =
      static_cast<const LcdDeltaBIntegrandParams*>(params);
 
  const double u = 1.00 + 2.00 * b * b;
  double value = b * std::pow(b * b / u, static_cast<double>(p->k));
  for (size_t i = p->k; i < p->q; ++i) value /= u;
 
  const double exponent = -p->cValue / (2.00 * u);
  return value * (p->reduced ? std::expm1(exponent) : std::exp(exponent));
}
 
inline double lcd_delta_b0(size_t d, double u, double eULeading, double eU,
                           double e0, double dEi, double c) {
  assert(c > 0.00);
 
  if (d == 0) return 0.25 * u * eULeading - 0.25 * e0 + 0.125 * c * dEi;
  if (d == 1) return -0.25 * dEi;
 
  double coefficient = 0.50;
  double cPower = c;      // c^(d-j+1), j = d
  double uPower = 1.00;   // u^(j-2),   j = d ... built downwards below
  for (size_t j = d; j > 2; --j) uPower *= u;
 
  double sum = 0.00;
  for (size_t j = d;; --j) {
    sum += coefficient * (eU / (cPower * uPower) - e0 / cPower);
    if (j == 2) break;
    coefficient *= static_cast<double>(j - 2) * 2.00;
    cPower *= c;
    uPower /= u;
  }
 
  return sum;
}
 
inline double lcd_delta_bkk(size_t k, double bMax, double c);
 
inline double lcd_delta_bkk1(size_t k, double bMax, double c);
 
inline double lcd_delta_bkk_closed(size_t k, double bMax, double c,
                                   size_t shift, bool reduced = false) {
  assert(k >= 1);
  assert(bMax > 0.00);
  assert(c > 0.00);
  assert(!reduced || shift == 0);  // no reduction is defined for dBkk1
 
  const double u = 1.00 + 2.00 * bMax * bMax;
  const double exponent = -c / (2.00 * u);
  const double eU = std::exp(exponent);
  const double eULeading = reduced ? std::expm1(exponent) : eU;
  const double e0 = std::exp(-c / 2.00);
  const double dEi = lcd_ei(exponent) - lcd_ei(-c / 2.00);
 
  double sum = 0.00;
  double binomial = 1.00;  // binom(k, j)
  for (size_t j = 0; j <= k; ++j) {
    const double term =
        lcd_delta_b0(j + shift, u, eULeading, eU, e0, dEi, 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_delta_bkk_zero(size_t k, double bMax) {
  assert(k >= 1);
  assert(bMax > 0.00);
 
  const double u = 1.00 + 2.00 * bMax * bMax;
 
  double sum = (u - 1.00) - static_cast<double>(k) * std::log(u);
 
  double binomial = static_cast<double>(k);  // binom(k, 1), stepped up below
  for (size_t j = 2; j <= k; ++j) {
    binomial *= static_cast<double>(k - j + 1) / static_cast<double>(j);
    const double dj = static_cast<double>(j);
    const double term =
        binomial * (std::pow(u, 1.00 - dj) - 1.00) / (1.00 - dj);
    sum += (j % 2 == 0) ? term : -term;
  }
 
  return std::ldexp(sum, -static_cast<int>(k) - 2);
}
 
inline double lcd_delta_bkk_zero_minus_leading(size_t k, double bMax) {
  assert(k >= 1);
  assert(bMax > 0.00);
 
  const double u = 1.00 + 2.00 * bMax * bMax;
 
  double sum = -1.00 - static_cast<double>(k) * std::log(u);
 
  double binomial = static_cast<double>(k);
  for (size_t j = 2; j <= k; ++j) {
    binomial *= static_cast<double>(k - j + 1) / static_cast<double>(j);
    const double dj = static_cast<double>(j);
    const double term =
        binomial * (std::pow(u, 1.00 - dj) - 1.00) / (1.00 - dj);
    sum += (j % 2 == 0) ? term : -term;
  }
 
  return std::ldexp(sum, -static_cast<int>(k) - 2);
}
 
inline double lcd_delta_bkk_reduced(size_t k, double bMax, double c) {
  assert(k >= 1);
  assert(bMax > 0.00);
  assert(c >= 0.00);
 
  if (c <= 0.00) return 0.00;  // dBkk(0) - dBkk(0)
 
  // the expm1 form has already dropped pow(2,-k) * U/4; take off what is left
  // of dBkk(k, 0), which is only O(ln U)
  return lcd_delta_bkk_closed(k, bMax, c, 0, true) -
         lcd_delta_bkk_zero_minus_leading(k, bMax);
}
 
inline double lcd_delta_bkk_impl(size_t k, double bMax, double c,
                                 size_t shift) {
  assert(k >= 1);
  assert(bMax > 0.00);
  assert(c >= 0.00);
 
  return lcd_delta_bkk_closed(k, bMax, c, shift);
}
 
inline double lcd_delta_bkk(size_t k, double bMax, double c) {
  return lcd_delta_bkk_impl(k, bMax, c, 0);
}
 
inline double lcd_delta_bkk1(size_t k, double bMax, double c) {
  return lcd_delta_bkk_impl(k, bMax, c, 1);
}
 
 
inline double lcd_reported_offset(size_t N, double bMax, double weightSum) {
  assert(N >= 2);
  assert(N % 2 == 0);
  assert(bMax > 0.00);
 
  const size_t k = N / 2;
  const double bMaxSqrd = bMax * bMax;
 
  if (N == 2) {
    // K = W^2 bMax^2/2 + bMax^2/2 - W bMax^2 - ln(1+bMax^2)/2
    //     + W ln(1+2bMax^2)/2
    //
    // Two rearrangements:
    //  - the three bMax^2 terms collapse to (1-W)^2 bMax^2 / 2;
    //  - W ln(1+2b^2)/2 - ln(1+b^2)/2 is split as
    //    W/2 * ln((1+2b^2)/(1+b^2)) + (W-1)/2 * ln(1+b^2), whose first factor
    //    is log1p(b^2/(1+b^2)).
    const double weightDefect = 1.00 - weightSum;
    return 0.50 * bMaxSqrd * weightDefect * weightDefect +
           0.50 * weightSum * std::log1p(bMaxSqrd / (1.00 + bMaxSqrd)) -
           0.50 * weightDefect * std::log1p(bMaxSqrd);
  }
 
  return lcd_a_n(N, bMax) -
         std::ldexp(weightSum * lcd_delta_bkk_zero(k, bMax),
                    static_cast<int>(k) + 1) +
         0.50 * weightSum * weightSum * bMaxSqrd;
}
 
#endif  // LCD_EVEN_CLOSED_FORM_H