primecount 0.2.1

Rust wrapper for https://github.com/kimwalisch/primecount
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

///
/// @file   S2_easy_xy.cpp
/// @brief  Test the computation of the easy special leaves
///         S2_easy(x, y) used in the Lagarias-Miller-Odlyzko
///         and Deleglise-Rivat prime counting algorithms.
///
/// Copyright (C) 2017 Kim Walisch, <kim.walisch@gmail.com>
///
/// This file is distributed under the BSD License. See the COPYING
/// file in the top level directory.
///

#include <primecount-internal.hpp>
#include <PhiTiny.hpp>
#include <PiTable.hpp>
#include <generate.hpp>
#include <imath.hpp>
#include <S.hpp>

#include <stdint.h>
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <vector>
#include <random>

using std::min;
using std::max;
using namespace primecount;

void check(bool OK)
{
  std::cout << "   " << (OK ? "OK" : "ERROR") << "\n";
  if (!OK)
    std::exit(1);
}

int main()
{
  int64_t max_x = 100000;
  double max_alpha = get_alpha_deleglise_rivat(max_x);
  int64_t max_x13 = iroot<3>(max_x);
  int64_t max_y = (int64_t) (max_alpha * max_x13);
  int threads = 1;

  {
    auto primes = generate_primes<int64_t>(max_y);
    PiTable pi(max_y, threads);

    // test small x
    for (int64_t i = 1; i < max_x; i++)
    {
      int64_t x = i;
      double alpha = get_alpha_deleglise_rivat(x);
      int64_t x13 = iroot<3>(x);
      int64_t y = (int64_t) (alpha * x13);
      int64_t z = x / y;
      int64_t c = PhiTiny::get_c(y);
      int64_t pi_sqrty = pi[isqrt(y)];
      int64_t pi_x13 = pi[x13];
      int64_t s2_easy = 0;

      for (int64_t b = max(c, pi_sqrty) + 1; b <= pi_x13; b++)
      {
        int64_t min_trivial = min(x / (primes[b] * primes[b]), y);
        int64_t min_sparse = max(z / primes[b], primes[b]);
        int64_t l = pi[min_trivial];

        for (; primes[l] > min_sparse; l--)
          s2_easy += pi[x / (primes[b] * primes[l])] - b + 2;
      }

      std::cout << "S2_easy(" << x << ", " << y << ") = " << s2_easy;
      check(s2_easy == S2_easy(x, y, z, c, threads));
    }
  }

  max_x = 100000000;
  max_alpha = get_alpha_deleglise_rivat(max_x);
  max_x13 = iroot<3>(max_x);
  max_y = (int64_t) (max_alpha * max_x13);

  {
    auto primes = generate_primes<int64_t>(max_y);
    PiTable pi(max_y, threads);

    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_int_distribution<int64_t> dist(1, max_x);

    // test random x
    for (int64_t i = 0; i < 10000; i++)
    {
      int64_t x = dist(gen);
      double alpha = get_alpha_deleglise_rivat(x);
      int64_t x13 = iroot<3>(x);
      int64_t y = (int64_t) (alpha * x13);
      int64_t z = x / y;
      int64_t c = PhiTiny::get_c(y);
      int64_t pi_sqrty = pi[isqrt(y)];
      int64_t pi_x13 = pi[x13];
      int64_t s2_easy = 0;

      for (int64_t b = max(c, pi_sqrty) + 1; b <= pi_x13; b++)
      {
        int64_t min_trivial = min(x / (primes[b] * primes[b]), y);
        int64_t min_sparse = max(z / primes[b], primes[b]);
        int64_t l = pi[min_trivial];

        for (; primes[l] > min_sparse; l--)
          s2_easy += pi[x / (primes[b] * primes[l])] - b + 2;
      }

      std::cout << "S2_easy(" << x << ", " << y << ") = " << s2_easy;
      check(s2_easy == S2_easy(x, y, z, c, threads));
    }
  }

  std::cout << std::endl;
  std::cout << "All tests passed successfully!" << std::endl;

  return 0;
}