pentagonal.cpp

Go to the documentation of this file.

/* -*- coding: latin-1 -*- */
/** \file common/sigmaodd/pentagonal.cpp (December 20, 2017)
 *
 * GPLv3 --- Copyright (C) 2017 Olivier Pirson
 * http://www.opimedia.be/
 */

// \cond
#include <cassert>

#include <map>
#include <vector>
// \endcond

#include "pentagonal.hpp"

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


namespace sigmaodd {

  /* *******************
   * Private constants *
   *********************/
  size_t sum_divisors__euler__cache_max_nb_ = 20000;

  size_t sum_divisors__euler__table_nb_ = 1024*64;



  /* *******************
   * Private variables *
   *********************/
  std::map<nat_type, nat_type> sum_divisors__euler__cache_;
  std::vector<nat_type> sum_divisors__euler__table_(sum_divisors__euler__table_nb_);



  /* *******************
   * Private prototype *
   *********************/

  nat_type
  sum_divisors_diff_(nat_type n, nat_type pi);



  /* ******************
   * Private function *
   ********************/

  nat_type
  sum_divisors_diff_(nat_type n, nat_type pi) {
    return (n > pi
            ? sum_divisors__euler(n - pi)
            : (n == pi
               ? n
               : 0));
  }



  /* ***********
   * Functions *
   *************/

  nat_type
  sum_divisors__euler(nat_type n) {
    assert(n != 0);

    if (n < sum_divisors__euler__table_.size()) {
      if (sum_divisors__euler__table_[n] != 0) {  // in the table
        return sum_divisors__euler__table_[n];
      }
    }
    else {
      const auto it = sum_divisors__euler__cache_.find(n);

      if (it != sum_divisors__euler__cache_.cend()) {  // in the cache
        return it->second;
      }
    }

    if (is_even(n)) {
      unsigned int alpha = 0;

      do {
        ++alpha;
        n >>= 1;
      } while (is_divide(2, n));

      return sum_divisors__euler(n)*sum_geometric_progression_strict(2, alpha);
    }

    for (unsigned int i = 0; i < 10; ++i) {
      const nat_type prime = odd_primes_table_by_index(i);

      if (is_divide(prime, n)) {
        unsigned int alpha = 0;

        do {
          ++alpha;
          n /= prime;
        } while (is_divide(prime, n));

        return sum_divisors__euler(n)*sum_geometric_progression_strict(prime, alpha);
      }
    }

    nat_type sum = 0;

    for (nat_type i = 1; ; ++i) {
      nat_type pi = pentagonal(i);
      nat_type pi_neg = pentagonal_neg(i);

      if ((n < pi) && (n < pi_neg)) {
        break;
      }

      if (is_odd(i)) {
        sum += sum_divisors_diff_(n, pi);
        sum += sum_divisors_diff_(n, pi_neg);
      }
      else {
        sum -= sum_divisors_diff_(n, pi);
        sum -= sum_divisors_diff_(n, pi_neg);
      }
    }

    if (n < sum_divisors__euler__table_.size()) {  // put in the table
      sum_divisors__euler__table_[n] = sum;
    }
    else {                                 // put in the cache
      if (sum_divisors__euler_cache_is_full()) {  // kill one item (at random)
        auto it = sum_divisors__euler__cache_.begin();

        std::advance(it,
                     static_cast<unsigned int>(std::rand())
                     % (sum_divisors__euler__cache_max_nb_ - 10));
        sum_divisors__euler__cache_.erase(it);
      }
      sum_divisors__euler__cache_[n] = sum;
    }

    return sum;
  }


  bool
  sum_divisors__euler_cache_is_full() {
    return (sum_divisors__euler__cache_.size() >= sum_divisors__euler__cache_max_nb_);
  }


  std::size_t
  sum_divisors__euler_cache_nb() {
    return sum_divisors__euler__cache_.size();
  }


  std::size_t
  sum_divisors__euler_table_nb() {
    return sum_divisors__euler__table_.size();
  }

}  // namespace sigmaodd