Problem 214 says:
Let φ be Euler’s totient function, i.e. for a natural number n, φ(n) is the number of k, 1 ≤ k ≤ n, for which gcd(k,n) = 1.
By iterating φ, each positive integer generates a decreasing chain of numbers ending in 1.
E.g. if we start with 5 the sequence 5,4,2,1 is generated.
Here is a listing of all chains with length 4:
5,4,2,1
7,6,2,1
8,4,2,1
9,6,2,1
10,4,2,1
12,4,2,1
14,6,2,1
18,6,2,1Only two of these chains start with a prime, their sum is 12.
What is the sum of all primes less than 40000000 which generate a chain of length 25?
My first approach is to build a list of primes less than 40000000, and keeps a hashtable of all the totients that the program solves for, as well as a hashtable for storing all of the previous chains the program has come across. It works but it takes an age to give the answer, even when using .net 4′s task parallel library. So I’ll be reworking this and posting a better version.
Until then, here is the code. (C# requires .net 4)
using System;
using System.Collections;
using System.Collections.Generic;
using System.Threading;
using System.Threading.Tasks;
namespace Problem214
{
class MainClass
{
public static long[] primes = GetPrimesUpTo(40000000);
public static Hashtable knownTotients = new Hashtable();
public static Hashtable KnownChains = new Hashtable();
private static object lockObject = new object();
private static object lockObject2 = new object();
public static void Main (string[] args)
{
long sum = 0L;
Console.WriteLine("starting search");
Parallel.ForEach(primes, p =>
{
//Console.WriteLine(p);
if (Chain25((int)p))
{
lock (lockObject)
{
sum = sum + p;
Console.WriteLine("Current Sum :: "+sum.ToString());
}
}
});
Console.WriteLine(sum);
Console.ReadLine();
}
static bool Chain25 (int num)
{
if (ChainLen(num) == 25)
return true;
else
return false;
}
static int ChainLen (int num)
{
int c = chain(num,1);
KnownChains.Add(num, c);
return c;
}
static int chain(int num, int len)
{
if (KnownChains.Contains(num))
return len + (int)KnownChains[num];
if (len >= 25) return 60000;
int phin = eulerTotientWrap(num);
if (phin != 1) return chain(phin,len+1);
else return len++;
}
static bool[] GetPrimeSieve(long upTo)
{
long sieveSize = upTo + 1;
bool[] sieve = new bool[sieveSize];
Array.Clear(sieve, 0, (int)sieveSize);
sieve[0] = true;
sieve[1] = true;
long p, max = (long)Math.Sqrt(sieveSize) + 1;
for (long i = 2; i <= max; i++)
{
if (sieve[i]) continue;
p = i + i;
while (p < sieveSize) { sieve[p] = true; p += i; }
}
return sieve;
}
static long[] GetPrimesUpTo(long upTo)
{
if (upTo < 2) return null;
bool[] sieve = GetPrimeSieve(upTo);
long[] primes = new long[upTo + 1];
long index = 0;
for (long i = 2; i <= upTo; i++) if (!sieve[i]) primes[index++] = i;
Array.Resize(ref primes, (int)index);
return primes;
}
static int GetLastIndex(int n)
{
for (int i = 0; i< primes.Length;i++)
{
if (primes[i] >= n)
return i+1;
}
return primes.Length;
}
static int eulerTotientWrap(int n)
{
if (knownTotients.Contains(n))
return (int)knownTotients[n];
int t = eulerTotient(n);
lock (lockObject2)
{
if (!knownTotients.Contains(n))
knownTotients.Add(n, t);
}
return t;
}
static int eulerTotient(int n)
{
int numPrimes = GetLastIndex(n);
int totient = n;
int currentNum = n, temp, p, prevP = 0;
for (int i = 0; i < numPrimes; i++)
{
p = (int)primes[i];
if (p > currentNum) break;
temp = currentNum / p;
if (temp * p == currentNum)
{
currentNum = temp;
i--;
if (prevP != p) { prevP = p; totient -= (totient / p); }
}
}
return totient;
}
}
}