accumulate_eulerphis.cpp

/*
 * compute sigma[k=1..N] eulerphi(k)
 * http://oeis.org/A064018
 * http://stackoverflow.com/a/1134851
 */

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <omp.h>
#include <boost/multiprecision/cpp_int.hpp>

using namespace std;
using boost::multiprecision::uint128_t;

const uint64_t N = 1e9;
const uint64_t sqrtN = sqrt(N);
const size_t segment_size = 256*1024/8;
const int nthreads = 12;

vector<uint64_t> small_primes;

void init(void)
{
  vector<bool> v(sqrtN+1, true);

  for (size_t i = 2, root = sqrt(sqrtN); i <= root; i++) {
    if (v[i] == false) continue;
    for (size_t j = i*i; j <= sqrtN; j += i) v[j] = false;
  }
  for (size_t i = 2; i <= sqrtN; i++) {
    if (v[i] == true) small_primes.push_back(i);
  }
  small_primes.push_back(N+1); // sentinel
}

template <typename T> T align(T start, T p)
{
  return (start + p - 1) / p * p;
}

// segment: [start, stop]
uint128_t accumulate_eulerphis_segmented(size_t start, size_t stop)
{
  vector<uint64_t> phis(stop-start+1, 1);

  for (size_t i = 0; small_primes[i] <= stop; i++) {
    const auto& p = small_primes[i];
    for (auto power = p; power <= stop; power *= p) {
      for (auto j = align(start, power); j <= stop; j += power) phis[j-start] *= p;
    }
  }

  for (size_t i = 0, n = start; i < phis.size(); i++, n++) {
    phis[i] = phis[i] == n ? 1 : n/phis[i]-1;
  }

  for (size_t i = 0; small_primes[i] <= stop; i++) {
    const auto& p = small_primes[i];
    for (auto j = align(start, p); j <= stop; j += p) phis[j-start] *= p-1;
    for (auto power = p*p; power <= stop; power *= p) {
      for (auto j = align(start, power); j <= stop; j += power) phis[j-start] *= p;
    }
  }

  return accumulate(begin(phis), end(phis), uint128_t(0));
}

uint128_t accumulate_eulerphis(void)
{
  const int nsegments = (N + segment_size - 1) / segment_size;
  uint128_t accs[nthreads] = {};

# pragma omp parallel for num_threads(nthreads) schedule(dynamic)
  for (int i = 0; i < nsegments; i++) {
    const size_t start = segment_size * i + 1;
    const size_t stop  = min(N, start + segment_size - 1);
    accs[omp_get_thread_num()] += accumulate_eulerphis_segmented(start, stop);
  }

  return accumulate(begin(accs), end(accs), uint128_t(0));
}

int main(void)
{
  init();
  cout << "sigma[k=1.." << N << "] eulerphi(k) = " << accumulate_eulerphis() << endl;
}