primes.cpp

Go to the documentation of this file.

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

// \cond
#include <cstdlib>

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

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


namespace sigmaodd {

  /* ***********
   * Constants *
   *************/

#ifdef PRIME16
  // See in ../../miscellaneous/offsets/
  const nat_type array_potential_prime_offsets_[array_potential_prime_offsets_nb_] = {
    10, 2, 4, 2, 4, 6, 2, 6, 4, 2,
    4, 6, 6, 2, 6, 4, 2, 6, 4, 6,
    8, 4, 2, 4, 2, 4, 8, 6, 4, 6,
    2, 4, 6, 2, 6, 6, 4, 2, 4, 6,
    2, 6, 4, 2, 4, 2, 10, 2};
#endif

  const std::string prime_filename = "big_data/prime28.bin";



  /* **********
   * Variable *
   ************/
  // http://primes.utm.edu/nthprime/
#ifdef PRIME16
  const prime_type array_odd_primes_[array_odd_primes_nb_ + 1] = {
#include "primes_table16.txt"
  };
#else
  prime_type* array_odd_primes_ = nullptr;
#endif



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

  std::vector<FactorExp>
  factorize(nat_type n) {
    assert(n != 0);
    assert(n <= MAX_POSSIBLE_N);
#ifdef PRIME16
    assert(n < 65536);
#endif

    std::vector<FactorExp> list;
    nat_type sqrt_n = floor_square_root(n);
    const prime_type* p = odd_primes_table_ptr() - 1;

    {
      std::pair<nat_type, unsigned int> n_exp = divide_until_odd_nb(n);

      if (n_exp.second > 0) {
        list.push_back({2, n_exp.second});
        n = n_exp.first;
      }
    }


    while (*(++p) <= sqrt_n) {
      assert(*p != 0);

      unsigned int alpha = 0;

      while (is_divide(*p, n)) {
        n /= *p;
        ++alpha;
      }

      if (alpha > 0) {
        list.push_back({*p, alpha});
        sqrt_n = floor_square_root(n);
      }
    }

    // Remain n_divided is prime
    if (n > 1) {
      list.push_back({n, 1});
    }

    return list;
  }


  nat_type
  factorization_to_n(std::vector<FactorExp> prime_exps) {
    nat_type n = 1;

    for (FactorExp prime_exp : prime_exps) {
      n *= pow_nat(prime_exp.factor, prime_exp.exp);
    }

    return n;
  }


  nat_type
  factorization_to_nu(std::vector<FactorExp> prime_exps) {
    nat_type nu = 1;

    for (FactorExp prime_exp : prime_exps) {
      nu *= prime_exp.exp + 1;
    }

    return nu;
  }


  nat_type
  factorization_to_nu_odd(std::vector<FactorExp> prime_exps) {
    nat_type nu = 1;

    for (FactorExp prime_exp : prime_exps) {
      if (prime_exp.factor != 2) {
        nu *= prime_exp.exp + 1;
      }
    }

    return nu;
  }


  nat_type
  factorization_to_sigma(std::vector<FactorExp> prime_exps) {
    nat_type sigma = 1;

    for (FactorExp prime_exp : prime_exps) {
      sigma *= sum_geometric_progression_strict(prime_exp.factor, prime_exp.exp);
    }

    return sigma;
  }


  nat_type
  factorization_to_sigma_odd(std::vector<FactorExp> prime_exps) {
    nat_type sigma = 1;

    for (FactorExp prime_exp : prime_exps) {
      if (prime_exp.factor != 2) {
        sigma *= sum_geometric_progression_strict(prime_exp.factor, prime_exp.exp);
      }
    }

    return sigma;
  }


  bool
  is_prime(nat_type n) {
    if (is_even(n) || (n <= 1)) {
      return (n == 2);
    }
    else if (n <= odd_primes_table_last()) {
      // Search prime present in the precalculated table
      return is_prime_in_odd_primes_table(n);
    }
    else {
      const nat_type sqrt_n = floor_square_root(n);

      // Search prime divisor present in the precalculated table
      for (unsigned int i = 0; i < odd_primes_table_nb(); ++i) {
        const nat_type prime = odd_primes_table_by_index(i);

        if (is_divide(prime, n)) {
          return (prime == n);
        }
        else if (prime > sqrt_n) {  // n is prime
          return true;
        }
      }

#ifdef PRIME16
      // Search other prime divisor by iteration on potential primes
      assert(odd_primes_table_last() % potential_prime_offsets_table_modulo() == 1);

      nat_type p = odd_primes_table_last();

      while (p < sqrt_n) {
        for (unsigned int i = 0; i < potential_prime_offsets_table_nb(); ++i) {
          p += potential_prime_offsets_table_by_index(i);

          if (is_divide(p, n)) {
            return (p == n);
          }
        }
      }
#endif

      return true;
    }
  }


  bool
  is_prime_in_odd_primes_table(nat_type n) {
    assert(sizeof(int32_t) == sizeof(prime_type));

    // Binary search in the table
    prime_type key = static_cast<prime_type>(n);

    return (nullptr
            != std::bsearch(&key, array_odd_primes_, odd_primes_table_nb(), sizeof(prime_type),
                            [](const void *ptr_a, const void *ptr_b) -> int {
                              const prime_type a = *static_cast<const prime_type*>(ptr_a);
                              const prime_type b = *static_cast<const prime_type*>(ptr_b);

                              return (a > b
                                      ? 1
                                      : (a < b
                                         ? -1
                                         : 0));
                            }));
  }


#ifdef PRIME16
  bool
  read_primes_table() {}
#else
  bool
  read_primes_table() {
    if (array_odd_primes_ != nullptr) {
      free(array_odd_primes_);
      array_odd_primes_ = nullptr;
    }

    if (sizeof(prime_type) != 4) {
      return false;
    }

    std::string filename = prime_filename;

    // Search the file from prime_filename
    for (unsigned int i = 0; i < 5; ++i) {
      if (helper::is_file_exists(filename)) {
        break;
      }

      filename = helper::concat_path("..", filename);
    }

    if (!helper::is_file_exists(filename)) {
      return false;
    }

    const std::size_t size = (odd_primes_table_nb() + 1)*4;

    array_odd_primes_ = reinterpret_cast<prime_type*>(malloc(size));

    if (array_odd_primes_ == nullptr) {
      return false;
    }

    std::ifstream infile(filename, std::ios::in | std::ios::binary);

    infile.read(reinterpret_cast<char*>(array_odd_primes_), 4);  // skip the first prime 2
    infile.read(reinterpret_cast<char*>(array_odd_primes_), size);  // read primes
    array_odd_primes_[odd_primes_table_nb()] = 0;

    infile.close();

    if (infile.fail()) {
      free(array_odd_primes_);
      array_odd_primes_ = nullptr;
    }

    return (array_odd_primes_ != nullptr);
  }
#endif

}  // namespace sigmaodd