公式解説の方法では \(pq^3\) の値が64bit整数で収まらない可能性があり、やや面倒です。 代わりに、次のような解法を使えばより簡単な方針で解くことができます。

• まず、素数表から \(q\) にあたる整数を決め打つ。
• 次に、条件を満たす \(p\) の値を考える。これは、 \([1, \min(q-1, \frac{N}{q^3})]\) を満たす素数全てである。
• ここで、\([1, \min(q-1, \frac{N}{q^3})]\) を満たす素数の個数は、素数の個数の累積和を用いれば \(O(1)\) で計算できる。
• 以上の流れを \(N^{\frac{1}{3}}\) 以下の全ての素数 \(q\) に対して行う。

素数の整数性から、 \(\frac{N}{q^3}\)は端数を切り捨てても構いません(つまり、整数型で扱っても問題ありません)。 これにより、演算中の全ての値が \(N\) 以下に収まります。

実装例(C++):

#include<algorithm>
#include<cmath>
#include<iostream>
#include<vector>
/**
max以下の素数を昇順に列挙します。
計算量はO(max)です。
制約: 0≦max≦10^9
*/
std::vector<int> getPrimes(int max) {
  std::vector<int> primes;
  std::vector<int> lpf(max + 1, 0);
  for (int i = 2;i <= max;++ i) {
    if (lpf[i] == 0) {
      lpf[i] = i;
      primes.push_back(i);
    }
    for (int p : primes) {
      int mul = p * i;
      if (mul > max || p > lpf[i]) break; // ベルトラン＝チェビシェフの定理より、p*iはmaxの2倍以下
      lpf[mul] = p;
    }
  }
  return primes;
}

int main() {
  long long N;
  std::cin >> N;
  int cbrtN = std::cbrt(N);
  std::vector<int> primes = getPrimes(cbrtN);
  std::vector<int> cumsum(1, 0); // cumsum[i]はi以下の素数の個数を持つ
  cumsum.reserve(cbrtN + 1);
  for (int i = 1, pi = 0, p = 2;i <= cbrtN;++ i) {
    int last = cumsum[cumsum.size() - 1];
    if (i == p) {
      ++ last;
      ++ pi;
      p = pi >= (int)primes.size() ? cbrtN + 1 : primes[pi];
    }
    cumsum.push_back(last);
  }
  long long ans = 0;
  for (int q : primes) {
    int pMax = (int)std::min(q - 1LL, N / q / q / q);
    ans += cumsum[pMax];
  }
  std::cout << ans << std::endl;
}