table_harmonic.cpp

Go to the documentation of this file.

/* -*- coding: latin-1 -*- */
/** \file common/table_harmonic.cpp (January 5, 2018)
 * \brief
 * Print same information in the same format than
 * ../../miscellaneous/harmonic/harmonic.py,
 * to compare and check correction.
 *
 * GPLv3 --- Copyright (C) 2017, 2018 Olivier Pirson
 * http://www.opimedia.be/
 */

// \cond
#include <cstdlib>

#include <iomanip>
#include <iostream>
// \endcond

#include "helper/helper.hpp"
#include "sigmaodd/harmonic.hpp"
#include "sigmaodd/divisors.hpp"
#include "sigmaodd/primes.hpp"
#include "sigmaodd/sigmaodd.hpp"



/* ***********
 * Prototype *
 *************/

void
help_and_exit();



/* **********
 * Function *
 ************/

void
help_and_exit() {
  std::cerr << "Usage: table_harmonic [options]" << std::endl
            << std::endl
            << "Method options to calculate sum of divisors:" << std::endl
            << "  --bound       factorize method with bound shortcuts (by default)" << std::endl
            << "  --factorize   factorize number with partial prime numbers table" << std::endl
            << "  --naive       try each odd possible divisor)" << std::endl
            << std::endl
            << "Options:" << std::endl
            << "  --first n   first number to check (1 by default)" << std::endl
            << "  --last  n   last number to check (50 by default)" << std::endl
            << "  --nb    n   number of numbers to check" << std::endl;

  exit(EXIT_FAILURE);
}



/* ******
 * Main *
 ********/

int
main(int argc, const char* const argv[]) {
  // Load primes table
  if (!sigmaodd::read_primes_table()) {
    std::cerr << "! Impossible to load \"" << sigmaodd::prime_filename << '"' << std::endl
              << std::endl;
    help_and_exit();
  }


  enum Methods: unsigned int {factorize, factorize_bound, naive};

  const std::string method_strings[3] =  {"factorize",
                                          "factorize bound",
                                          "naive"};

  sigmaodd::nat_type first = 1;
  sigmaodd::nat_type last = 50;
  unsigned int method = factorize;


  // Read command line parameters
  for (unsigned int i = 1; i < static_cast<unsigned int>(argc); ++i) {
    const std::string param(argv[i]);

    if (param == "--bound") {
      method = factorize_bound;
    }
    else if (param == "--factorize") {
      method = factorize;
    }
    else if (param == "--first") {
      first = std::max(1ul, helper::get_ulong(argc, argv, ++i, &help_and_exit));
    }
    else if (param == "--last") {
      last = helper::get_ulong(argc, argv, ++i, &help_and_exit);
    }
    else if (param == "--naive") {
      method = naive;
    }
    else if (param == "--nb") {
      last = first + helper::get_ulong(argc, argv, ++i, &help_and_exit) - 1;
    }
    else {
      help_and_exit();
    }
  }


  // Print parameters (to stderr)
  std::cerr << "First: " << first
            << "\tLast: " << last
            << "\tNb: " << (first <= last
                            ? (last - first)/2 + 1
                            : 0)
            << "\tMethod: " << method_strings[method]
            << std::endl;


  // Print table legend
  std::cout << "n\tsqrt\t|s_odd\tbound\tGUB\tDeTemp.\t|sum_o.\t|Half bound\tdiff h.\tfl.harm._odd\tharmonic_odd"
            << std::endl;


  // Calculate and print data
  std::cout << std::fixed << std::setprecision(10);

  for (sigmaodd::nat_type n = first; n <= last; ++n) {
    const sigmaodd::nat_type sqrt_n = sigmaodd::floor_square_root(n);
    const sigmaodd::nat_type sqrt_n1 = sigmaodd::floor_square_root(n - 1);

    sigmaodd::nat_type exact;

    if (method == factorize) {
      exact = sigmaodd::sum_odd_divisors__factorize(n);
    }
    else if (method == factorize_bound) {
      exit(EXIT_FAILURE);
    }
    else if (method == naive) {
      exact = sigmaodd::sum_odd_divisors__naive(n);
    }
    else {
      exit(EXIT_FAILURE);
    }

    const sigmaodd::nat_type bound
      = (sigmaodd::sum_odd(sqrt_n) +
         sigmaodd::sum_floor_n_harmonic_odd(n, sigmaodd::floor_square_root(n - 1)));
    const sigmaodd::nat_type gub_bound
      = sigmaodd::sum_odd(sqrt_n)
      + static_cast<sigmaodd::nat_type>(std::floor(sigmaodd::harmonic_odd(sqrt_n1)
                                                   * static_cast<double>(n)));
    const sigmaodd::nat_type bound_detemple = sigmaodd::sum_odd_divisors_upper_bound__DeTemple(n);
    const sigmaodd::nat_type sum_odd_n = sigmaodd::sum_odd(n);

    const double diff_half = (sqrt_n1 == 0
                              ? 0
                              : sigmaodd::diff_half_harmonic_upper_bound(sqrt_n1));
    const double h_odd = sigmaodd::harmonic_odd(n);

    std::cout << n << '\t'
              << sqrt_n << '\t'
              << exact << '\t'
              << bound << '\t'
              << gub_bound << '\t'
              << bound_detemple << '\t'
              << sum_odd_n << '\t'
              << sigmaodd::diff_half_harmonic_upper_bound(n) << '\t'
              << sigmaodd::sum_floor_n_harmonic_odd(n, n) << '\t'
              << diff_half << '\t'
              << h_odd
              << std::endl;

    assert(exact <= bound);
    assert(exact <= gub_bound);
    assert(exact <= bound_detemple);

    assert(bound <= gub_bound);
    assert(gub_bound <= bound_detemple);

    assert(exact <= sum_odd_n);
    if ((n != 2) && (n != 4)) {
      assert(bound <= sum_odd_n);
    }
  }

  return EXIT_SUCCESS;
}