Super speedy totient function The goal is simple: calculate the totient function for as many numbers as you can in 10 seconds and sum the numbers. 

You must print your result at the end and you must actually calculate it. No automated totient function is allowed, but bignum libraries are. You have to start at 1 and count up all the integers consecutively. You are not allowed to skip numbers.

Your score is how many numbers your program can calculate on your machine / how many my program can calculate on your machine. My code is a simple program in C++ (optimizations off), hopefully you can run it.

Important properties you could use!


if gcd(m,n) = 1, phi(mn) = phi(m) * phi(n)
if p is prime, phi(p) = p - 1 (for p &lt; 10^20)
if n is even, phi(2n) = 2 phi(n)
others listed in first link


My code

#include &lt;iostream&gt;
using namespace std;

int gcd(int a, int b)
{
    while (b != 0)
    {
        int c = a % b;
        a = b;
        b = c;
    }
    return a;
}

int phi(int n)
{
    int x = 0;
    for (int i=1; i&lt;=n; i++)
    {
        if (gcd(n, i) == 1)
            x++;
    }
    return x;
}

int main()
{
    unsigned int sum = 0;
    for (int i=1; i&lt;19000; i++) // Change this so it runs in 10 seconds
    {
        sum += phi(i);
    }
        cout &lt;&lt; sum &lt;&lt; endl;
        return 0;
}
Nimrod: ~38,667 (580,000,000/15,000)

This answer uses a pretty simple approach. The code employs a simple prime number sieve that stores the prime of the smallest prime power in each slot for composite numbers (zero for primes), then uses dynamic programming to construct the totient function over the same range, then sums the results. The program spends virtually all its time constructing the sieve, then calculates the totient function in a fraction of the time. It looks like it comes down to constructing an efficient sieve (with the slight twist that one also has to be able to extract a prime factor for composite numbers from the result and has to keep memory usage at a reasonable level).

Update: Improved performance by reducing memory footprint and improving cache behavior. It's possible to squeeze out 5%-10% more performance, but the increase in code complexity is not worth it. Ultimately, this algorithm primarily exercises a CPU's von Neumann bottleneck, and there are very few algorithmic tweaks that can get around that.

Also updated the performance divisor since the C++ code wasn't meant to be compiled with all optimizations on and nobody else did it. :)

Update 2: Optimized sieve operation for improved memory access. Now handling small primes in bulk via memcpy() (~5% speedup) and skipping multiples of 2, 3, and 5 when sieving bigger primes (~10% speedup).

C++ code: 9.9 seconds (with g++ 4.9)

Nimrod code: 9.9 seconds (with -d:release, gcc 4.9 backend)

proc handleSmallPrimes(sieve: var openarray[int32], m: int) =
  # Small primes are handled as a special case through what is ideally
  # the system's highly optimized memcpy() routine.
  let k = 2*3*5*7*11*13*17
  var sp = newSeq[int32](k div 2)
  for i in [3,5,7,11,13,17]:
    for j in countup(i, k, 2*i):
      sp[j div 2] = int32(i)
  for i in countup(0, sieve.high, len(sp)):
    if i + len(sp) &lt;= len(sieve):
      copyMem(addr(sieve[i]), addr(sp[0]), sizeof(int32)*len(sp))
    else:
      copyMem(addr(sieve[i]), addr(sp[0]), sizeof(int32)*(len(sieve)-i))
  # Fixing up the numbers for values that are actually prime.
  for i in [3,5,7,11,13,17]:
    sieve[i div 2] = 0

proc constructSieve(m: int): seq[int32] =
  result = newSeq[int32](m div 2 + 1)
  handleSmallPrimes(result, m)
  var i = 19
  # Having handled small primes, we only consider candidates for
  # composite numbers that are relatively prime with 31. This cuts
  # their number almost in half.
  let steps = [ 1, 7, 11, 13, 17, 19, 23, 29, 31 ]
  var isteps: array[8, int]
  while i * i &lt;= m:
    if result[i div 2] == 0:
      for j in 0..7: isteps[j] = i*(steps[j+1]-steps[j])
      var k = 1 # second entry in "steps mod 30" list.
      var j = 7*i
      while j &lt;= m:
        result[j div 2] = int32(i)
        j += isteps[k]
        k = (k + 1) and 7 # "mod 30" list has eight elements.
    i += 2

proc calculateAndSumTotients(sieve: var openarray[int32], n: int): int =
  result = 1
  for i in 2'i32..int32(n):
    var tot: int32
    if (i and 1) == 0:
      var m = i div 2
      var pp: int32 = 2
      while (m and 1) == 0:
        pp *= 2
        m = m div 2
      if m == 1:
        tot = pp div 2
      else:
        tot = (pp div 2) * sieve[m div 2]
    elif sieve[i div 2] == 0: # prime?
      tot = i - 1
      sieve[i div 2] = tot
    else:
      # find and extract the first prime power pp.
      # It's relatively prime with i/pp.
      var p = sieve[i div 2]
      var m = i div p
      var pp = p
      while m mod p == 0 and m != p:
        pp *= p
        m = m div p
      if m == p: # is i a prime power?
        tot = pp*(p-1)
      else:
        tot = sieve[pp div 2] * sieve[m div 2]
      sieve[i div 2] = tot
    result += tot

proc main(n: int) =
  var sieve = constructSieve(n)
  let totSum = calculateAndSumTotients(sieve, n)
  echo totSum

main(580_000_000)