<commit_before>///
/// @file  phi_vector.cpp
/// @brief The PhiCache class calculates the partial sieve function
///        (a.k.a. Legendre-sum) using the recursive formula:
///        phi(x, a) = phi(x, a - 1) - phi(x / primes[a], a - 1).
///        phi(x, a) counts the numbers <= x that are not divisible by
///        any of the first a primes. The algorithm used is an
///        optimized version of the algorithm described in Tomás
///        Oliveira e Silva's paper [1]. I have added 5 optimizations
///        to my implementation which significantly speed up the
///        calculation:
///
///        * Cache results of phi(x, a)
///        * Calculate phi(x, a) using formula [2] if a <= 6
///        * Calculate phi(x, a) using pi(x) lookup table
///        * Calculate all phi(x, a) = 1 upfront
///        * Stop recursion at c instead of 1
///
///       [1] <NAME>, Computing pi(x): the combinatorial
///           method, Revista do DETUA, vol. 4, no. 6, March 2006, p. 761.
///           http://sweet.ua.pt/tos/bib/5.4.pdf
///       [2] phi(x, a) = (x / pp) * φ(pp) + phi(x % pp, a)
///           with pp = 2 * 3 * ... * prime[a] 
///
/// Copyright (C) 2017 <NAME>, <<EMAIL>>
///
/// This file is distributed under the BSD License. See the COPYING
/// file in the top level directory.
///

#include <PiTable.hpp>
#include <primecount-internal.hpp>
#include <fast_div.hpp>
#include <imath.hpp>
#include <PhiTiny.hpp>
#include <min.hpp>

#include <stdint.h>
#include <algorithm>
#include <array>
#include <cstdlib>
#include <vector>
#include <limits>

#ifdef _OPENMP
  #include <omp.h>
#endif

using namespace std;
using namespace primecount;

namespace {

/// Cache phi(x, a) results if a < MAX_A
const int MAX_A = 100;

class PhiCache
{
public:
  PhiCache(vector<int64_t>& primes,
           PiTable& pi) :
    primes_(primes),
    pi_(pi)
  { }

  /// Calculate phi(x, a) using the recursive formula:
  /// phi(x, a) = phi(x, a - 1) - phi(x / primes_[a], a - 1)
  ///
  template <int SIGN>
  int64_t phi(int64_t x, int64_t a)
  {
    if (x <= primes_[a])
      return SIGN;
    else if (is_phi_tiny(a))
      return phi_tiny(x, a) * SIGN;
    else if (is_pix(x, a))
      return (pi_[x] - a + 1) * SIGN;
    else if (is_cached(x, a))
      return cache_[a][x] * SIGN;

    int64_t sqrtx = isqrt(x);
    int64_t pi_sqrtx = a;
    int64_t c = PhiTiny::get_c(sqrtx);
    int64_t sum = 0;

    if (sqrtx < pi_.size() && sqrtx < primes_[a])
      pi_sqrtx = pi_[sqrtx];

    // Move out of the loop the calculations where phi(x2, i) = 1
    // phi(x, a) = 1 if primes_[a] >= x
    // x2 = x / primes_[i + 1]
    // phi(x2, i) = 1 if primes_[i] >= x / primes_[i + 1]
    // phi(x2, i) = 1 if primes_[i] >= sqrt(x)
    // phi(x2, i) = 1 if i >= pi(sqrt(x))
    // \sum_{i = pi(sqrt(x))}^{a - 1} phi(x2, i) = a - pi(sqrt(x))
    //
    sum += (pi_sqrtx - a) * SIGN;
    sum += phi_tiny(x, c) * SIGN;

    for (int64_t i = c; i < pi_sqrtx; i++)
    {
      int64_t x2 = fast_div(x, primes_[i + 1]);

      if (is_pix(x2, i))
        sum += (pi_[x2] - i + 1) * -SIGN;
      else
        sum += phi<-SIGN>(x2, i);
    }

    update_cache(x, a, sum);

    return sum;
  }

private:
  using T = uint16_t;
  array<vector<T>, MAX_A> cache_;
  vector<int64_t>& primes_;
  PiTable& pi_;

  void update_cache(uint64_t x, uint64_t a, int64_t sum)
  {
    if (a < cache_.size() &&
        x <= numeric_limits<T>::max())
    {
      if (x >= cache_[a].size())
        cache_[a].resize(x + 1, 0);

      cache_[a][x] = (T) abs(sum);
    }
  }

  bool is_pix(int64_t x, int64_t a) const
  {
    return x < pi_.size() &&
           x < isquare(primes_[a + 1]);
  }

  bool is_cached(uint64_t x, uint64_t a) const
  {
    return a < cache_.size() && 
           x < cache_[a].size() && 
           cache_[a][x];
  }
};

} // namespace

namespace primecount {

/// Returns a vector with phi(x, i - 1) values such that
/// phi[i] = phi(x, i - 1) for 1 <= i <= a.
/// phi(x, a) counts the numbers <= x that are not
/// divisible by any of the first a primes.
///
vector<int64_t> phi_vector(int64_t x,
                           int64_t a,
                           vector<int64_t>& primes,
                           PiTable& pi, 
                           int threads)
{
  vector<int64_t> phi(a + 1, 0);

  if (x > 0 && a > 0)
  {
    phi[1] = x;
    PhiCache cache(primes, pi);

    int64_t thread_threshold = ipow(10ll, 10);
    threads = ideal_num_threads(threads, x, thread_threshold);

    #pragma omp parallel for num_threads(threads) schedule(dynamic, 16) firstprivate(cache)
    for (int64_t i = 2; i <= a; i++)
      phi[i] = cache.phi<-1>(x / primes[i - 1], i - 2);

    // calculate phi(x, a) using partial results
    for (int64_t i = 2; i <= a; i++)
      phi[i] += phi[i - 1];
  }

  return phi;
}

} // namespace<commit_msg>Fix phi(x, a) = 1 if (prime[a] >= x)<commit_after>///
/// @file  phi_vector.cpp
/// @brief The PhiCache class calculates the partial sieve function
///        (a.k.a. Legendre-sum) using the recursive formula:
///        phi(x, a) = phi(x, a - 1) - phi(x / primes[a], a - 1).
///        phi(x, a) counts the numbers <= x that are not divisible by
///        any of the first a primes. The algorithm used is an
///        optimized version of the algorithm described in Tomás
///        Oliveira e Silva's paper [1]. I have added 5 optimizations
///        to my implementation which significantly speed up the
///        calculation:
///
///        * Cache results of phi(x, a)
///        * Calculate phi(x, a) using formula [2] if a <= 6
///        * Calculate phi(x, a) using pi(x) lookup table
///        * Calculate all phi(x, a) = 1 upfront
///        * Stop recursion at c instead of 1
///
///       [1] <NAME>, Computing pi(x): the combinatorial
///           method, Revista do DETUA, vol. 4, no. 6, March 2006, p. 761.
///           http://sweet.ua.pt/tos/bib/5.4.pdf
///       [2] phi(x, a) = (x / pp) * φ(pp) + phi(x % pp, a)
///           with pp = 2 * 3 * ... * prime[a] 
///
/// Copyright (C) 2017 <NAME>, <<EMAIL>>
///
/// This file is distributed under the BSD License. See the COPYING
/// file in the top level directory.
///

#include <PiTable.hpp>
#include <primecount-internal.hpp>
#include <fast_div.hpp>
#include <imath.hpp>
#include <PhiTiny.hpp>
#include <min.hpp>

#include <stdint.h>
#include <algorithm>
#include <array>
#include <cstdlib>
#include <vector>
#include <limits>

#ifdef _OPENMP
  #include <omp.h>
#endif

using namespace std;
using namespace primecount;

namespace {

/// Cache phi(x, a) results if a < MAX_A
const int MAX_A = 100;

class PhiCache
{
public:
  PhiCache(vector<int64_t>& primes,
           PiTable& pi) :
    primes_(primes),
    pi_(pi)
  { }

  /// Calculate phi(x, a) using the recursive formula:
  /// phi(x, a) = phi(x, a - 1) - phi(x / primes_[a], a - 1)
  ///
  template <int SIGN>
  int64_t phi(int64_t x, int64_t a)
  {
    if (x <= primes_[a])
      return SIGN;
    else if (is_phi_tiny(a))
      return phi_tiny(x, a) * SIGN;
    else if (is_pix(x, a))
      return (pi_[x] - a + 1) * SIGN;
    else if (is_cached(x, a))
      return cache_[a][x] * SIGN;

    int64_t sqrtx = isqrt(x);
    int64_t pi_sqrtx = a;
    int64_t c = PhiTiny::get_c(sqrtx);
    int64_t sum = 0;

    if (sqrtx < pi_.size() && sqrtx < primes_[a])
      pi_sqrtx = pi_[sqrtx];

    // Move out of the loop the calculations where phi(x2, i) = 1
    // phi(x, a) = 1 if primes_[a] >= x
    // x2 = x / primes_[i + 1]
    // phi(x2, i) = 1 if primes_[i] >= x / primes_[i + 1]
    // phi(x2, i) = 1 if primes_[i] >= sqrt(x)
    // phi(x2, i) = 1 if i >= pi(sqrt(x))
    // \sum_{i = pi(sqrt(x))}^{a - 1} phi(x2, i) = a - pi(sqrt(x))
    //
    sum += (pi_sqrtx - a) * SIGN;
    sum += phi_tiny(x, c) * SIGN;

    for (int64_t i = c; i < pi_sqrtx; i++)
    {
      int64_t x2 = fast_div(x, primes_[i + 1]);

      if (is_pix(x2, i))
        sum += (pi_[x2] - i + 1) * -SIGN;
      else
        sum += phi<-SIGN>(x2, i);
    }

    update_cache(x, a, sum);

    return sum;
  }

private:
  using T = uint16_t;
  array<vector<T>, MAX_A> cache_;
  vector<int64_t>& primes_;
  PiTable& pi_;

  void update_cache(uint64_t x, uint64_t a, int64_t sum)
  {
    if (a < cache_.size() &&
        x <= numeric_limits<T>::max())
    {
      if (x >= cache_[a].size())
        cache_[a].resize(x + 1, 0);

      cache_[a][x] = (T) abs(sum);
    }
  }

  bool is_pix(int64_t x, int64_t a) const
  {
    return x < pi_.size() &&
           x < isquare(primes_[a + 1]);
  }

  bool is_cached(uint64_t x, uint64_t a) const
  {
    return a < cache_.size() && 
           x < cache_[a].size() && 
           cache_[a][x];
  }
};

} // namespace

namespace primecount {

/// Returns a vector with phi(x, i - 1) values such that
/// phi[i] = phi(x, i - 1) for 1 <= i <= a.
/// phi(x, a) counts the numbers <= x that are not
/// divisible by any of the first a primes.
///
vector<int64_t> phi_vector(int64_t x,
                           int64_t a,
                           vector<int64_t>& primes,
                           PiTable& pi, 
                           int threads)
{
  int64_t size = a + 1;
  vector<int64_t> phi(size, x > 0);

  if (primes[a] > x)
    a = pi[x];

  if (x > 0 && a > 0)
  {
    phi[1] = x;
    int64_t thread_threshold = ipow(10ll, 10);
    threads = ideal_num_threads(threads, x, thread_threshold);
    PhiCache cache(primes, pi);

    #pragma omp parallel for num_threads(threads) schedule(dynamic, 16) firstprivate(cache)
    for (int64_t i = 2; i <= a; i++)
      phi[i] = cache.phi<-1>(x / primes[i - 1], i - 2);

    // calculate phi(x, a) using partial results
    for (int64_t i = 2; i <= a; i++)
      phi[i] += phi[i - 1];
  }

  return phi;
}

} // namespace
{'fraction_non_alphanumeric': 0.16772606110313718, 'fraction_numerical': 0.027168341193356572, 'mean_word_length': 2.073408947700063, 'pattern_counts': {'":': 0, '<': 95, '<?xml version=': 0, '>': 76, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '477282', 'n_tokens_mistral': 4017, 'n_tokens_neox': 3717, 'n_words': 1102}