<commit_before>///
/// @file  pi_deleglise_rivat3.cpp
/// @brief Implementation of the Lagarias-Miller-Odlyzko prime counting
///        algorithm with the improvements of Deleglise and Rivat.
///        This version uses compression (see FactorTable & PiTable)
///        to reduce the memory usage.
///
/// Copyright (C) 2014 <NAME>, <<EMAIL>>
///
/// This file is distributed under the BSD License. See the COPYING
/// file in the top level directory.
///

#include <PiTable.hpp>
#include <FactorTable.hpp>
#include <primecount-internal.hpp>
#include <primesieve.hpp>
#include <BitSieve.hpp>
#include <pmath.hpp>
#include <PhiTiny.hpp>
#include <tos_counters.hpp>

#include <stdint.h>
#include <algorithm>
#include <vector>

using namespace std;
using namespace primecount;

namespace {

/// Cross-off the multiples of prime in the sieve array.
/// For each element that is unmarked the first time update
/// the special counters tree data structure.
///
template <typename T1, typename T2>
void cross_off(int64_t prime, int64_t low, int64_t high, int64_t& next_multiple, T1& sieve, T2& counters)
{
  int64_t segment_size = sieve.size();
  int64_t k = next_multiple;

  for (; k < high; k += prime * 2)
  {
    if (sieve[k - low])
    {
      sieve.unset(k - low);
      cnt_update(counters, k - low, segment_size);
    }
  }
  next_multiple = k;
}

/// Calculate the contribution of the special leaves.
/// @see ../docs/computing-special-leaves.md
/// @pre y > 0 && c > 1
///
int64_t S2(int64_t x,
           int64_t y,
           int64_t z,
           int64_t c,
           vector<int32_t>& primes,
           FactorTable& factors)
{
  PiTable pi(y);
  int64_t pi_y = pi(y);
  int64_t pi_sqrty = pi(isqrt(y));
  int64_t pi_sqrtz = pi(min(isqrt(z), y));
  int64_t limit = z + 1;
  int64_t segment_size = next_power_of_2(isqrt(limit));
  int64_t S2_result = 0;

  BitSieve sieve(segment_size);
  vector<int32_t> counters(segment_size);
  vector<int64_t> next(primes.begin(), primes.begin() + pi_sqrtz + 1);
  vector<int64_t> phi(pi_sqrtz + 1, 0);

  // Segmented sieve of Eratosthenes
  for (int64_t low = 1; low < limit; low += segment_size)
  {
    // Current segment = interval [low, high[
    int64_t high = min(low + segment_size, limit);
    int64_t b = c + 1;

    // check if we need the sieve
    if (c < pi_sqrtz)
    {
      sieve.memset(low);

      // phi(y, i) nodes with i <= c do not contribute to S2, so we
      // simply sieve out the multiples of the first c primes
      for (int64_t i = 2; i <= c; i++)
      {
        int64_t k = next[i];
        for (int64_t prime = primes[i]; k < high; k += prime * 2)
          sieve.unset(k - low);
        next[i] = k;
      }

      // Initialize special tree data structure from sieve
      cnt_finit(sieve, counters, segment_size);
    }

    // For c + 1 <= b <= pi_sqrty
    // Find all special leaves: n = primes[b] * m, with mu[m] != 0 and primes[b] < lpf[m]
    // which satisfy: low <= (x / n) < high
    for (; b <= pi_sqrty; b++)
    {
      int64_t prime = primes[b];
      int64_t min_m = max(x / (prime * high), y / prime);
      int64_t max_m = min(x / (prime * low), y);

      if (prime >= max_m)
        goto next_segment;

      FactorTable::to_index(&min_m);
      FactorTable::to_index(&max_m);

      for (int64_t m = max_m; m > min_m; m--)
      {
        if (prime < factors.lpf(m))
        {
          int64_t n = prime * factors.get_number(m);
          int64_t count = cnt_query(counters, (x / n) - low);
          int64_t phi_xn = phi[b] + count;
          S2_result -= factors.mu(m) * phi_xn;
        }
      }

      phi[b] += cnt_query(counters, (high - 1) - low);
      cross_off(prime, low, high, next[b], sieve, counters);
    }

    // For pi_sqrty <= b < pi_y
    // Find all special leaves: n = primes[b] * primes[l]
    // which satisfy: low <= (x / n) < high
    for (; b < pi_y; b++)
    {
      int64_t prime = primes[b];
      int64_t l = pi(min(x / (prime * low), y));

      if (prime >= primes[l])
        goto next_segment;

      int64_t min_hard_leaf = max(x / (prime * high), y / prime);
      min_hard_leaf = in_between(prime, min_hard_leaf, y);
      int64_t min_trivial_leaf = min(x / (prime * prime), y);
      int64_t min_clustered_easy_leaf = min(isqrt(x / prime), y);
      int64_t min_sparse_easy_leaf = min(z / prime, y);

      min_trivial_leaf = max(min_hard_leaf, min_trivial_leaf);
      min_clustered_easy_leaf = max(min_hard_leaf, min_clustered_easy_leaf);
      min_sparse_easy_leaf = max(min_hard_leaf, min_sparse_easy_leaf);

      // Find all trivial leaves which satisfy:
      // phi(x / (primes[b] * primes[l]), b - 1) = 1
      if (primes[l] > min_trivial_leaf)
      {
        int64_t l_min = pi(min_trivial_leaf);
        S2_result += l - l_min;
        l = l_min;
      }

      // Find all clustered easy leaves which satisfy:
      // x / n <= y such that phi(x / n, b - 1) = pi(x / n) - b + 2
      // And phi(x / n, b - 1) == phi(x / m, b - 1)
      while (primes[l] > min_clustered_easy_leaf)
      {
        int64_t n = prime * primes[l];
        int64_t xn = x / n;
        int64_t phi_xn = pi(xn) - b + 2;
        int64_t m = prime * primes[b + phi_xn - 1];
        int64_t xm = max(x / m, min_clustered_easy_leaf);
        int64_t l2 = pi(xm);
        S2_result += phi_xn * (l - l2);
        l = l2;
      }

      // Find all sparse easy leaves which satisfy:
      // x / n <= y such that phi(x / n, b - 1) = pi(x / n) - b + 2
      for (; primes[l] > min_sparse_easy_leaf; l--)
      {
        int64_t n = prime * primes[l];
        int64_t xn = x / n;
        S2_result += pi(xn) - b + 2;
      }

      if (b <= pi_sqrtz)
      {
        // Find all hard leaves which satisfy:
        // low <= (x / n) < high
        for (; primes[l] > min_hard_leaf; l--)
        {
          int64_t n = prime * primes[l];
          int64_t xn = x / n;
          int64_t count = cnt_query(counters, xn - low);
          int64_t phi_xn = phi[b] + count;
          S2_result += phi_xn;
        }

        phi[b] += cnt_query(counters, (high - 1) - low);
        cross_off(prime, low, high, next[b], sieve, counters);
      }
    }

    next_segment:;
  }<commit_msg>Add assertions for easy leaves<commit_after>///
/// @file  pi_deleglise_rivat3.cpp
/// @brief Implementation of the Lagarias-Miller-Odlyzko prime counting
///        algorithm with the improvements of Deleglise and Rivat.
///        This version uses compression (see FactorTable & PiTable)
///        to reduce the memory usage.
///
/// Copyright (C) 2014 <NAME>, <<EMAIL>>
///
/// This file is distributed under the BSD License. See the COPYING
/// file in the top level directory.
///

#include <PiTable.hpp>
#include <FactorTable.hpp>
#include <primecount-internal.hpp>
#include <primesieve.hpp>
#include <BitSieve.hpp>
#include <pmath.hpp>
#include <PhiTiny.hpp>
#include <tos_counters.hpp>

#include <stdint.h>
#include <algorithm>
#include <vector>

using namespace std;
using namespace primecount;

namespace {

/// Cross-off the multiples of prime in the sieve array.
/// For each element that is unmarked the first time update
/// the special counters tree data structure.
///
template <typename T1, typename T2>
void cross_off(int64_t prime, int64_t low, int64_t high, int64_t& next_multiple, T1& sieve, T2& counters)
{
  int64_t segment_size = sieve.size();
  int64_t k = next_multiple;

  for (; k < high; k += prime * 2)
  {
    if (sieve[k - low])
    {
      sieve.unset(k - low);
      cnt_update(counters, k - low, segment_size);
    }
  }
  next_multiple = k;
}

/// Calculate the contribution of the special leaves.
/// @see ../docs/computing-special-leaves.md
/// @pre y > 0 && c > 1
///
int64_t S2(int64_t x,
           int64_t y,
           int64_t z,
           int64_t c,
           vector<int32_t>& primes,
           FactorTable& factors)
{
  PiTable pi(y);
  int64_t pi_y = pi(y);
  int64_t pi_sqrty = pi(isqrt(y));
  int64_t pi_sqrtz = pi(min(isqrt(z), y));
  int64_t limit = z + 1;
  int64_t segment_size = next_power_of_2(isqrt(limit));
  int64_t S2_result = 0;

  BitSieve sieve(segment_size);
  vector<int32_t> counters(segment_size);
  vector<int64_t> next(primes.begin(), primes.begin() + pi_sqrtz + 1);
  vector<int64_t> phi(pi_sqrtz + 1, 0);

  // Segmented sieve of Eratosthenes
  for (int64_t low = 1; low < limit; low += segment_size)
  {
    // Current segment = interval [low, high[
    int64_t high = min(low + segment_size, limit);
    int64_t b = c + 1;

    // check if we need the sieve
    if (c < pi_sqrtz)
    {
      sieve.memset(low);

      // phi(y, i) nodes with i <= c do not contribute to S2, so we
      // simply sieve out the multiples of the first c primes
      for (int64_t i = 2; i <= c; i++)
      {
        int64_t k = next[i];
        for (int64_t prime = primes[i]; k < high; k += prime * 2)
          sieve.unset(k - low);
        next[i] = k;
      }

      // Initialize special tree data structure from sieve
      cnt_finit(sieve, counters, segment_size);
    }

    // For c + 1 <= b <= pi_sqrty
    // Find all special leaves: n = primes[b] * m, with mu[m] != 0 and primes[b] < lpf[m]
    // which satisfy: low <= (x / n) < high
    for (; b <= pi_sqrty; b++)
    {
      int64_t prime = primes[b];
      int64_t min_m = max(x / (prime * high), y / prime);
      int64_t max_m = min(x / (prime * low), y);

      if (prime >= max_m)
        goto next_segment;

      FactorTable::to_index(&min_m);
      FactorTable::to_index(&max_m);

      for (int64_t m = max_m; m > min_m; m--)
      {
        if (prime < factors.lpf(m))
        {
          int64_t n = prime * factors.get_number(m);
          int64_t count = cnt_query(counters, (x / n) - low);
          int64_t phi_xn = phi[b] + count;
          S2_result -= factors.mu(m) * phi_xn;
        }
      }

      phi[b] += cnt_query(counters, (high - 1) - low);
      cross_off(prime, low, high, next[b], sieve, counters);
    }

    // For pi_sqrty <= b < pi_y
    // Find all special leaves: n = primes[b] * primes[l]
    // which satisfy: low <= (x / n) < high
    for (; b < pi_y; b++)
    {
      int64_t prime = primes[b];
      int64_t l = pi(min(x / (prime * low), y));

      if (prime >= primes[l])
        goto next_segment;

      int64_t min_hard_leaf = max(x / (prime * high), y / prime);
      min_hard_leaf = in_between(prime, min_hard_leaf, y);
      int64_t min_trivial_leaf = min(x / (prime * prime), y);
      int64_t min_clustered_easy_leaf = min(isqrt(x / prime), y);
      int64_t min_sparse_easy_leaf = min(z / prime, y);

      min_trivial_leaf = max(min_hard_leaf, min_trivial_leaf);
      min_clustered_easy_leaf = max(min_hard_leaf, min_clustered_easy_leaf);
      min_sparse_easy_leaf = max(min_hard_leaf, min_sparse_easy_leaf);

      // Find all trivial leaves which satisfy:
      // phi(x / (primes[b] * primes[l]), b - 1) = 1
      if (primes[l] > min_trivial_leaf)
      {
        int64_t l_min = pi(min_trivial_leaf);
        S2_result += l - l_min;
        l = l_min;
      }

      // Find all clustered easy leaves which satisfy:
      // x / n <= y such that phi(x / n, b - 1) = pi(x / n) - b + 2
      // And phi(x / n, b - 1) == phi(x / m, b - 1)
      while (primes[l] > min_clustered_easy_leaf)
      {
        int64_t n = prime * primes[l];
        int64_t xn = x / n;
        assert(xn < isquare(primes[b]));
        int64_t phi_xn = pi(xn) - b + 2;
        int64_t m = prime * primes[b + phi_xn - 1];
        int64_t xm = max(x / m, min_clustered_easy_leaf);
        int64_t l2 = pi(xm);
        S2_result += phi_xn * (l - l2);
        l = l2;
      }

      // Find all sparse easy leaves which satisfy:
      // x / n <= y such that phi(x / n, b - 1) = pi(x / n) - b + 2
      for (; primes[l] > min_sparse_easy_leaf; l--)
      {
        int64_t n = prime * primes[l];
        int64_t xn = x / n;
        assert(xn < isquare(primes[b]));
        S2_result += pi(xn) - b + 2;
      }

      if (b <= pi_sqrtz)
      {
        // Find all hard leaves which satisfy:
        // low <= (x / n) < high
        for (; primes[l] > min_hard_leaf; l--)
        {
          int64_t n = prime * primes[l];
          int64_t xn = x / n;
          int64_t count = cnt_query(counters, xn - low);
          int64_t phi_xn = phi[b] + count;
          S2_result += phi_xn;
        }

        phi[b] += cnt_query(counters, (high - 1) - low);
        cross_off(prime, low, high, next[b], sieve, counters);
      }
    }

    next_segment:;
  }
{'fraction_non_alphanumeric': 0.13802479008396643, 'fraction_numerical': 0.02463014794082367, 'mean_word_length': 2.1716966776566067, 'pattern_counts': {'":': 0, '<': 89, '<?xml version=': 0, '>': 59, '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': '1680151', 'n_tokens_mistral': 4980, 'n_tokens_neox': 4657, 'n_words': 1274}