bounds.cpp

Go to the documentation of this file.

/* -*- coding: latin-1 -*- */
/** \file common/bounds.cpp (January 5, 2018)
 * \brief
 * Some upper bounds for the sigma_odd problem.
 *
 * GPLv3 --- Copyright (C) 2017, 2018 Olivier Pirson
 * http://www.opimedia.be/
 */

// \cond
#include <cassert>
#include <cmath>
#include <cstdlib>

#include <iostream>
// \endcond

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


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

void
help_and_exit();



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

void
help_and_exit() {
  std::cerr << "Usage: bounds [options]" << std::endl
            << std::endl
            << "Options:" << std::endl
            << "  --first n   first odd number to check (1 by default)" << std::endl
            << "  --last  n   last odd number to check (101 by default)" << std::endl
            << "  --nb    n   number of odd 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();
  }


  sigmaodd::nat_type first = 1;
  sigmaodd::nat_type last = 1000;


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

    if (param == "--first") {
      first = std::max(1ul, helper::get_ulong(argc, argv, ++i, &help_and_exit));
      if (sigmaodd::is_even(first)) {
        ++first;
      }
    }
    else if (param == "--last") {
      last = helper::get_ulong(argc, argv, ++i, &help_and_exit);
      if (sigmaodd::is_even(last)) {
        --last;
      }
    }
    else if (param == "--nb") {
      last = first + helper::get_ulong(argc, argv, ++i, &help_and_exit)*2 - 1;
    }
    else {
      help_and_exit();
    }
  }


  // Print intern configuration
  helper::print_intern_config_compiler();


  // Print parameters
  std::cout << "First: " << first
            << "\tLast: " << last
            << "\tNb: " << (first <= last
                            ? (last - first)/2 + 1
                            : 0)
            << std::endl;


  // Print table legend
  std::cout << std::endl
            << "n\tsigma_odd\tvarsigma_odd\tUpper bound\tDeTemple\tDeTemple adapt\tSum_odd\tPow9_8\tSqrt8\tHalf_sqrt8"
            << std::endl;


  // Print data
  for (sigmaodd::nat_type n = first; n <= last; n += 2) {
    const sigmaodd::nat_type sum_odd_divisors = sigmaodd::sum_odd_divisors__factorize(n);
    const sigmaodd::nat_type varsum_odd_divisors = sigmaodd::divide_until_odd(sigmaodd::sum_odd_divisors__factorize(n));  // ???
    const sigmaodd::nat_type bound = sigmaodd::sum_odd_divisors_upper_bound(n);
    const sigmaodd::nat_type bound_DeTemple = sigmaodd::sum_odd_divisors_upper_bound__DeTemple(n);
    const sigmaodd::nat_type sum_odd = sigmaodd::sum_odd(n);

    // ??????
    sigmaodd::nat_type prime = sigmaodd::first_divisor(n);
    sigmaodd::nat_type n_divided;
    unsigned int alpha;
    sigmaodd::nat_type bound_DeTemple_adapt;

    assert(sigmaodd::is_odd(n));

    {
      const std::pair<sigmaodd::nat_type, unsigned int>
        result_alpha = (prime != 1
                        ? sigmaodd::divide_until_nb(n, prime)
                        : std::make_pair(n, 0u));

      if (n == 1) {
        bound_DeTemple_adapt = 1;
      }
      else {
        assert(prime > 1);

        if (prime == n) {  // n is prime
          assert(sigmaodd::is_prime(n));

          // ???????? bound_DeTemple_adapt = n + 1;
          bound_DeTemple_adapt = 0;
        }
        else if (prime <= sigmaodd::floor_square_root(result_alpha.first)) {
          assert(!sigmaodd::is_prime(n));
          // assert(prime <= sigmaodd::floor_square_root(result_alpha.first));

          n_divided = result_alpha.first;
          alpha = result_alpha.second;

          assert(sigmaodd::is_odd(n_divided));
          assert(1 <= alpha);

          bound_DeTemple_adapt = sigmaodd::divide_until_odd(sigmaodd::sum_geometric_progression(prime, alpha))*sigmaodd::sum_odd_divisors_upper_bound__DeTemple(n_divided, prime);
        }
        else {  // n_divided is prime
          assert(!sigmaodd::is_prime(n));
          // ???assert(sigmaodd::is_square(n_divided));

          // ???bound_DeTemple_adapt = sigmaodd::sum_odd_divisors__factorize(n);
          // ???bound_DeTemple_adapt = sigmaodd::sum_odd_divisors_upper_bound__DeTemple(n);
          bound_DeTemple_adapt = 0;
        }
      }
    }

    const sigmaodd::nat_type pow9_8
      = static_cast<sigmaodd::nat_type>(std::floor(std::pow(n, 9.0/8)*2));
    const sigmaodd::nat_type sqrt8
      = static_cast<sigmaodd::nat_type>(std::floor((std::pow(n - 1, 1.0/8) + 1.0/2)
                                                   * static_cast<double>(n)));
    const sigmaodd::nat_type half_sqrt8
      = static_cast<sigmaodd::nat_type>(std::floor((std::pow(n - 1, 1.0/8) + 1.0/2)
                                                   * static_cast<double>(n) / 2));

    std::cout << n
              << '\t' << sum_odd_divisors
              << '\t' << varsum_odd_divisors
              << '\t' << bound
              << '\t' << bound_DeTemple
              << '\t' << ((bound_DeTemple_adapt > 0)  // ??? && (n_divided())
                          ? std::to_string(bound_DeTemple_adapt)
                          : "")
              << '\t' << sum_odd
              << '\t' << pow9_8
              << '\t' << sqrt8
              << '\t' << half_sqrt8
              << std::endl;

    assert(varsum_odd_divisors <= sum_odd_divisors);

    assert(sum_odd_divisors <= bound);
    assert(sum_odd_divisors <= bound_DeTemple);
    if (bound_DeTemple_adapt > 0) {
      // ???assert(sum_odd_divisors <= bound_DeTemple_adapt);
      assert(varsum_odd_divisors <= bound_DeTemple_adapt);
    }
    assert(sum_odd_divisors <= pow9_8);
    if (n != 1) {
      assert(sum_odd_divisors <= sqrt8);
    }
    if (!sigmaodd:: is_square(n)) {
      assert(varsum_odd_divisors <= half_sqrt8);
    }

    assert(bound <= bound_DeTemple);
    assert(bound_DeTemple <= sum_odd);
    assert(bound_DeTemple_adapt <= sum_odd);

    assert(sqrt8 <= pow9_8);
    assert(half_sqrt8 <= sqrt8);
  }

  return EXIT_SUCCESS;
}