Home Problems Project Euler : Totient Chains

Project Euler : Totient Chains

0 94

Problem 214 of Project Euler requires one to find the sum of all the prime numbers upto 40000000 that have a totient chain of length 25.

I solved this problem by brute-force using memoization to prevent solving the same sub-problem over and over again. The values I memoized are:

  • totient_chain_length(n)
  • totient(n)
  • primefactors(n)

Also I precomputed all the primes till 40000000 so that I do not have to waste time doing that while computing. All the primes were stored in memory. The effect of such high memoization and prime number storage was that my process used almost 1.25 gb of memory.

One final thing to point out is the property of totient function for prime numbers.The totient value of a prime number n is n-1. There’s yet another property that could have further increased the performance – totient(n) * totient(m) = totient(m*n) if gcd(m,n) = 1.

All one needs to do is to compute the totient value of each number from 1 to 40000000. Then starting from n = 1 through n = 4000000, we need to update the value of totient_chain_length as totient_chain_length(n) = totient_chain_length(totient(n)) + 1.

Instead of computing all the values, I instead chose to compute and memoize them as required. Since only the chains starting with prime numbers are required, I got a generator of primes till 40000000. Iterating over the elements of the generator, I updated the totient_chain_length as mentioned above, calculating the missing values whenever required.

Here’s the main script that does this computation:

Created on Sep 28, 2012

@author: anuvrat
from utils.prime_utils import totient, primesfrom2to
from utils.memoize import Memoize

def totient_chain_length( n ):
    if n in range( 21 ): return ( 1, 1, 2, 3, 3, 4, 3, 4, 4, 4, 4, 5, 4, 5, 4, 5, 5, 6, 4, 5, 5 )[n]
    return 1 + totient_chain_length( totient( n ) )

total = 0
for prime in primesfrom2to( 40000000 ):
    if prime < 9000000: continue
    if totient_chain_length( prime - 1 ) == 24: total += prime

print( total )

The general purpose memoize class is:

Created on Sep 28, 2012

@author: anuvrat

class Memoize( object ):

    def __init__( self, f ):
        self.f = f
        self.cache = {}

    def __call__( self, *args ):
        if not args in self.cache: self.cache[args] = self.f( *args )
        return self.cache[args]

And the method used for computing the totient values is (note that it can be further optimized):

def totient( n ):
    if n == 0: return 1
    if is_probable_prime( n ): return n - 1

    tot = 1
    for p, exp in Counter( primefactors( n ) ).items():
        tot *= ( p - 1 ) * p ** ( exp - 1 )

    return tot

is_probable_prime actually is the Miller-Rabin test with a default trials value set to 5. However since I had already precomputed all the prime numbers till 40000000 it was a small matter of checking against the generated list.

Hi, my name's Anuvrat. I'm a Computer Science graduate from IIT Kharagpur. I have worked on the CART algorithm while I was at FICO and am now trying to make sense of the huge amount of data we have at Amazon Associates. In my free time I attempt Project Euler problems - the mathematical nature of those problems appeal to me. I am a big Chelsea fan - following them since the time Scholari was made the manager, only to be fired 8 months later, sadly. Following the ups and downs at that club has taught me quite a lot about life. After Chelsea, it's definitely Schumi and rock music that makes up my life.