helper.cpp

Go to the documentation of this file.

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

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

#include <algorithm>
#include <fstream>
// \endcond

#include "../helper/helper.hpp"
#include "helper.hpp"


namespace sigmaodd {

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

  nat_type
  floor_eighth_root(nat_type n) {
    if (n >= 17878103347812890625ul) {  // >= (2^(64/8) - 1)^8 = 17878103347812890625
      return 255;  // = 2^(64/8) - 1
    }
    else {
      nat_type root
        = static_cast<nat_type>(std::floor(std::sqrt(std::sqrt(std::sqrt(static_cast<double>(n))))));

      // Correct possible rounding error due to floating point computation
      while (square(square(square(root))) <= n) {  // get the first value too big
        ++root;
      }
      do  {  // get the correct value
        --root;
      } while (square(square(square(root))) > n);

      assert(square(square(square(root))) <= n);
      assert(n < square(square(square(root + 1))));

      return root;
    }
  }


  nat_type
  floor_fourth_root(nat_type n) {
    if (n >= 18445618199572250625ul) {  // >= (2^(64/4) - 1)^4 = 18445618199572250625
      return 65535;  // = 2^(64/4) - 1
    }
    else {
      nat_type root = static_cast<nat_type>(std::floor(std::sqrt(std::sqrt(static_cast<double>(n)))));

      // Correct possible rounding error due to floating point computation
      while (square(square(root)) <= n) {  // get the first value too big
        ++root;
      }
      do  {  // get the correct value
        --root;
      } while (square(square(root)) > n);

      assert(square(square(root)) <= n);
      assert(n < square(square(root + 1)));

      return root;
    }
  }


  nat_type
  floor_square_root(nat_type n) {
    if (n >= 18446744065119617025ul) {  // >= (2^(64/2) - 1)^2 = 18446744065119617025
      return 4294967295;  // = 2^(64/2) - 1
    }
    else {
      nat_type sqrt_n = static_cast<nat_type>(floor(sqrt(static_cast<double>(n))));

      /* Correct possible rounding error due to floating point computation */
      while (square(sqrt_n) <= n) {  /* get the first value too big */
        ++sqrt_n;
      }
      do  {  /* get the correct value */
        --sqrt_n;
      } while (square(sqrt_n) > n);

      assert(square(sqrt_n) <= n);
      assert(n < square(sqrt_n + 1));

      return sqrt_n;
    }
  }


  std::vector<nat_type>
  load_bad_table(const std::string &filename,
                 const std::vector<nat_type> &bad_table) {
    std::vector<nat_type> new_bad_table = load_numbers(filename);

    if (!new_bad_table.empty()) {
      new_bad_table.insert(new_bad_table.begin(), bad_table.cbegin(), bad_table.cend());
    }

    sort(new_bad_table.begin(), new_bad_table.end());

    return new_bad_table;
  }


  std::vector<nat_type>
  load_numbers(const std::string &filename) {
    std::vector<nat_type> list;
    std::ifstream infile(filename, std::ios::in);

    if (infile.is_open()) {
      std::string line;

      while (getline(infile, line)) {
        const auto i = line.find_first_not_of("0123456789");
        const std::string s(i == std::string::npos
                            ? line
                            : (i > 0
                               ? line.substr(0, i)
                               : ""));

        if (!s.empty()) {
          const nat_type n = std::stoul(s);

          list.push_back(n);
        }
      }

      infile.close();
    }

    return list;
  }
}  // namespace sigmaodd