harmonic.hpp

Go to the documentation of this file.

/* -*- coding: latin-1 -*- */
/** \file common/sigmaodd/harmonic.hpp (February 14, 2018)
 * \brief
 * Function to calculate harmonic number
 * H_n = 1 + 1/2 + 1/3 + 1/4 + ... + 1/n
 * and some variants and upper bounds.
 *
 * http://mathworld.wolfram.com/HarmonicNumber.html
 *
 * GPLv3 --- Copyright (C) 2017, 2018 Olivier Pirson
 * http://www.opimedia.be/
 */

#ifndef PROGS_SRC_COMMON_SIGMAODD_HARMONIC_HPP_
#define PROGS_SRC_COMMON_SIGMAODD_HARMONIC_HPP_

#include "helper.hpp"


namespace sigmaodd {

  /* **********
   * Constant *
   ************/
  /** \brief
   * The Euler-Mascheroni constant 0.577215664901532860606512090082402431042...
   * (rounded to an upper bound in double)
   * http://mathworld.wolfram.com/Euler-MascheroniConstant.html
   */
  extern const double euler_gamma;



  /* ************
   * Prototypes *
   **************/

  /** \brief
   * Return an upper bound of H_a - 1/2 H_b.
   *
   * Use Chen 2016 inequalities: https://link.springer.com/article/10.1007%2Fs11075-016-0116-9
   *
   * H_a - 1/2 H_b
   * < 1/2 ln( (2a+1)^2 / (2(2b+1)) ) + gamma/2 + 1/(24 (a^2 + a + beta)) - 1/(48 (b^2 + b + alpha))
   *
   * @param a >= 1
   * @param b >= 1
   *
   * @return an upper bound of H_a - 1/2 H_b
   */
  double
  diff_half_harmonic_upper_bound(nat_type a, nat_type b);


  /** \brief
   * Return an upper bound of H_n - 1/2 H_k with k = floor(n/2).
   *
   * Use Chen 2016 inequalities: https://link.springer.com/article/10.1007%2Fs11075-016-0116-9
   *
   * H_n - 1/2 H_k
   * < 1/2 ln( (2n+1)^2 / (2(n+1)) ) + gamma/2
   *   + 1/(24 (n^2 + n + beta)) - 1/(12 ((n+1)^2 + alpha_bis))   if n even
   *
   * H_n - 1/2 H_k
   * < 1/2 ln( (2n+1)^2 / (2n)) )    + gamma/2
   *   + 1/(24 (n^2 + n + beta)) - 1/(12 (n^2     + alpha_bis))   if n odd
   *
   * @param n >= 1
   *
   * @return an upper bound of H_n - 1/2 H_{n/2}
   */
  double
  diff_half_harmonic_upper_bound(nat_type n);


  /** \brief
   * Return an upper bound of H_a - H_b.
   *
   * Use Chen 2016 inequalities: https://link.springer.com/article/10.1007%2Fs11075-016-0116-9
   *
   * H_a - H_b < ln((2a + 1)/(2b + 1)) + 1/(24 (a^2 + a + beta)) - 1/(24 (b^2 + b + alpha))
   *
   * @param a >= 1
   * @param b >= 1
   *
   * @return an upper bound of H_a - H_b
   */
  double
  diff_harmonic_upper_bound(nat_type a, nat_type b);


  /** \brief
   * Return the harmonic number H_n = 1/1 + 1/2 + 1/3 + 1/4 + ... + 1/n.
   *
   * H_0 = 0
   *
   * http://mathworld.wolfram.com/HarmonicNumber.html
   *
   * @return the harmonic number H_n
   */
  constexpr
  double
  harmonic(nat_type n);


  /** \brief
   * Return 1/2 + 1/4 + 1/6 + 1/8 + ... + (1/n or 1/(n-1)).
   */
  constexpr
  double
  harmonic_even(nat_type n);


  /** \brief
   * Return 1/1 + 1/3 + 1/5 + 1/7 + ... + (1/n or 1/(n-1)).
   */
  constexpr
  double
  harmonic_odd(nat_type n);


  /** \brief
   * Return a lower bound of H_n.
   *
   * Use Chen 2016 inequalities: https://link.springer.com/article/10.1007%2Fs11075-016-0116-9
   *
   * ln(n + 1/2) + gamma + 1/(24 (n^2 + n + alpha)) < H_n
   *
   * @return an upper bound of H_n
   */
  double
  harmonic_lower_bound(nat_type n);


  /** \brief
   * Return an upper bound of H_n.
   *
   * Use Chen 2016 inequalities: https://link.springer.com/article/10.1007%2Fs11075-016-0116-9
   *
   * H_n <= ln(n + 1/2) + gamma + 1/(24 (n^2 + n + beta))
   *
   * @return an upper bound of H_n
   */
  double
  harmonic_upper_bound(nat_type n);


  /** \brief
   * Return floor(n/1) + floor(n/3) + floor(n/5) + floor(n/7) + ... + (n/to_n or floor(1/(to_n-1))).
   */
  nat_type
  sum_floor_n_harmonic_odd(nat_type n, nat_type to_n);

}  // namespace sigmaodd


#include "harmonic__inline.hpp"

#endif  // PROGS_SRC_COMMON_SIGMAODD_HARMONIC_HPP_