harmonic.cpp

Go to the documentation of this file.

/* -*- coding: latin-1 -*- */
/** \file common/sigmaodd/harmonic.cpp (February 14, 2018)
 *
 * GPLv3 --- Copyright (C) 2017, 2018 Olivier Pirson
 * http://www.opimedia.be/
 */

// \cond
#include <cassert>
#include <cfenv>
#include <cmath>
// \endcond

#include "divisors.hpp"
#include "harmonic.hpp"


namespace sigmaodd {

  /* *******************
   * Private constants *
   *********************/

  const double alpha = 0.425;
  const double alpha_bis = 0.7;
  const double beta = 0.40580406331047491619301581522449851;



  /* **********
   * Constant *
   ************/

  const double euler_gamma = 0.577215664901532860606512090082402431042;



  /* ***********
   * Functions *
   *************/
  double
  diff_half_harmonic_upper_bound(nat_type a, nat_type b) {
    assert(a != 0);
    assert(b != 0);

    const int round = fegetround();

    fesetround(FE_UPWARD);

    // 1/2 ln( (2a+1)^2 / (2(2b+1)) ) + gamma/2 + 1/(24 (a^2 + a + beta)) - 1/(48 (b^2 + b + alpha))
    const double bound = ((std::log(static_cast<double>(square(a*2 + 1))
                                    / static_cast<double>((b*2 + 1)*2))
                           + euler_gamma) / 2
                          - (0.5/(static_cast<double>(square(b) + b) + alpha)
                             - 1.0/(static_cast<double>(square(a) + a) + beta)) / 24.0);

    fesetround(round);

    return bound;
  }


  double
  diff_half_harmonic_upper_bound(nat_type n) {
    assert(n != 0);

    const int round = fegetround();

    fesetround(FE_UPWARD);

    // 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
    // 1/2 ln( (2n+1)^2 / (2n)) )    + gamma/2
    //   + 1/(24 (n^2 + n + beta)) - 1/(12 (n^2     + alpha_bis))   if n odd
    const double bound = (is_even(n)
                          ? ((std::log(static_cast<double>(square(n*2 + 1))
                                       / static_cast<double>((n + 1)*2))
                              + euler_gamma) / 2
                             - (1.0/(static_cast<double>(square(n + 1)) + alpha_bis)
                                - 0.5/(static_cast<double>(square(n) + n) + beta)) / 12)
                          : ((std::log(static_cast<double>(square(n*2 + 1))
                                       / static_cast<double>(n*2))
                              + euler_gamma) / 2
                             - (1.0/(static_cast<double>(square(n)) + alpha_bis)
                                - 0.5/(static_cast<double>(square(n) + n) + beta)) / 12));

    fesetround(round);

    return bound;
  }


  double
  diff_harmonic_upper_bound(nat_type a, nat_type b) {
    assert(a != 0);
    assert(b != 0);

    const int round = fegetround();

    fesetround(FE_UPWARD);

    // ln((2a + 1)/(2b + 1)) + 1/(24 (a^2 + a + beta)) - 1/(24 (b^2 + b + alpha))
    const double bound = (std::log(static_cast<double>(a*2 + 1)
                                   / static_cast<double>(b*2 + 1))
                          - (1.0/(static_cast<double>(square(b) + b) + alpha)
                             - 1.0/(static_cast<double>(square(a) + a) + beta)) / 24);

    fesetround(round);

    return bound;
  }


  double
  harmonic_lower_bound(nat_type n) {
    assert(n != 0);

    const int round = fegetround();

    fesetround(FE_DOWNWARD);

    // ln(n + 1/2) + gamma + 1/(24 (n^2 + n + alpha))
    const double bound = (std::log(static_cast<double>(n) + 0.5)
                          + euler_gamma
                          + 1.0/((static_cast<double>(square(n) + n) + alpha)*24));

    fesetround(round);

    return bound;
  }


  double
  harmonic_upper_bound(nat_type n) {
    assert(n != 0);

    const int round = fegetround();

    fesetround(FE_UPWARD);

    // ln(n + 1/2) + gamma + 1/(24 (n^2 + n + beta))
    const double bound = (std::log(static_cast<double>(n) + 0.5)
                          + euler_gamma
                          + 1.0/((static_cast<double>(square(n) + n) + beta)*24));

    fesetround(round);

    return bound;
  }


  nat_type
  sum_floor_n_harmonic_odd(nat_type n, nat_type to_n) {
    nat_type sum = 0;

    for (nat_type i = 1; i <= to_n; i += 2) {
      sum += static_cast<nat_type>(std::floor(static_cast<double>(n)/static_cast<double>(i)));
    }

    return sum;
  }

}  // namespace sigmaodd