Tag Archives: totient

Euler Problem 243 – Always the Totient

OK, so I had to look up the answer to this one. In my defence I did lots of research and determined that the resilience is related (once again) to Euler’s Totient function, specifically the relationship R(d) = phi(d)/(d-1). I figured that by using this relationship (and looking for R(d) lower than the ratio given) was the trick instead of having to brute force it. Nope. That is the brute force method. Ran for days and got nowhere.

So I broke down and looked up a solution. There aren’t very many for this one. Most were in interative languages like Ruby, which don’t help me. I did find one amazing solution is Haskell here, but for the life of me I don’t know why it works. I’ve adapted it and replaced as much as I can with clearer functions to help explain it. Perhaps not as efficient, but it still runs in less than a second.

Step One: Generate lists of primes and totients.
Step Two: Generate a list primorials consisting of tuples based on products of primes and totients.
Step Three: Use a recursive routine starting with the primorials list and a multiplier (m) of 1. If the current list head * m == the next primorial, then drop the head and loop. If the current list head * m is than the target then return it, otherwise increment m and loop.
Step Four: Shrug shoulders, open arms with palms upward, crinkle forehead, take a deep sign, and say, “WTF!”

import Data.List.Ordered

primes :: [Integer]
primes = 2 : sieve [3,5..] where
	sieve [] = []
	sieve (p:xs) = p : sieve (xs `minus` [p*p, p*p+2*p..])

totient :: Integer -> Integer
totient n = fromIntegral $ length [x | x <- [1..n], gcd x n == 1]

primorials :: [(Integer, Integer)]
primorials = [(product $ take x primes, product $ map totient $ take x primes) | x <- [1..]]

findRes :: [(Integer, Integer)] -> Integer -> Integer
findRes ((cpri,cphi):(npri,nphi):pps) m
	| cpri * m == npri = findRes ((npri,nphi):pps) 1
	| fromIntegral (cphi * m) / fromIntegral ((cpri * m) - 1) < ratio = cpri * m
	| otherwise = findRes ((cpri,cphi):(npri,nphi):pps) (m + 1)
	where ratio = 15499/94744

main = print $ findRes primorials 1

Euler Problem 72

Just add up the totient functions between 2 and 10E6.

#Euler72
class Algorithms:
    def totient(self, n):
        result = float(n)
        i = 2
        while i*i <= n:
            if n % i == 0:
                result -= result/i
            while n % i == 0:
                n /= i
            i += 1
        if n > 1:
            result -= result/n
        return result
class Problem72:
    def __init__(self, limit):
        self.limit = limit
        self.Euler = Algorithms()
    def Solution(self):
        count = 0
        for n in range (2, self.limit+1):
            count += int(self.Euler.totient(n))
            if n % 50000 == 0:
                print n
        return "Answer:"+str(count)
if __name__ == '__main__':
    P = Problem72(1000000)
    print P.Solution()

Euler Problem 70

So I could have kept calculating the totient function directly and left it to run for years. Instead, I googled totient and discovered an interesting little mathematical quirk that would enable me to significantly reduce my testing set.
If a number is a product of two prime numbers …

  • n = p1 * p2

then the totient of that number is …

  • phi (n) = (p1-1) * (p2-1)

The square root of 10E7 is about 3100 plus change. So I calculated all the prime numbers 30% above and below that number. Then tested each of those numbers whereby the products were < 10E7. Smallest n/phi(n) wins. Routine runs in about 5 minutes. By the way, I wrote this routine to leverage (and re-learn Pythonian OOP).

#Euler70
from math import sqrt
class Algorithms:
    def permutate(self, sequence):
        if not sequence:
            return [sequence]
        else:
            set = []
            for k in range(len(sequence)):
                part = sequence[:k] + sequence[k+1:]
                for m in self.permutate(part):
                    set.append(sequence[k:k+1] + m)
            return set 
    def primeGen (self, lowerLimit, upperLimit):
        isPrime = [False, True]
        i = 2
        while i < upperLimit:
            isPrime.append(True)
            i += 1
        p = 2
        while p * p < upperLimit:
            if isPrime[p] == True:
                m = p * p
                while m < upperLimit:
                    if m % p == 0:
                        isPrime[m] = False
                    m += p
            p += 1
        P = []
        count = 0
        for p in isPrime:
            if p and count >= lowerLimit :
                P.append(count)
            count += 1
        return P
class Problem70:
    def __init__(self, limit):
        self.limit = limit
        self.Euler = Algorithms()
    def Solution(self):
        primes = self.Euler.primeGen(.7 * sqrt(self.limit), 1.3 * sqrt(self.limit))
        print "Primes:"+str(primes)+"\n"
        i = 0
        minNDivPhi = 10
        minN = 0
        for p1 in primes:
          print "Testing:"+str(p1)
          i += 1
          for p2 in primes[i:]:
            n = p1 * p2
            if n > self.limit:
                break
            phi = (p1-1) * (p2-1)
            nDivPhi = n / float(phi)
            if nDivPhi < minNDivPhi:
                if str(int(n)) in self.Euler.permutate (str(int(phi))):
                    minNDivPhi = nDivPhi
                    minN = n
                    print "n:"+str(n)+" phi:"+str(phi)+" nDivPhi:"+str(nDivPhi)
        return "Value of n where nDivPhi is minimized:"+str(minN)
if __name__ == '__main__':
    P = Problem70(10000000)
    print P.Solution()

Euler Problem 69

I tried this first with a GCD function and found it too slow. In fact it would have taken nearly a month to run. So I looked for a more efficient algorithm. Surprisingly, I found this one in Java (ironic isn’t it?) then translated it to Python.

I must admit, working out Euler Problems is a bit easier with Python. I suppose manipulating lists etc. in a higher level scripting language requires less work. Java was a bit too system-level for this kind of work. So this routine is pretty slick. It runs through the entire set in under 90 seconds. Next.

#Euler69
def totient(n):
    result = float(n)
    i = 2
    while i*i <= n:
        if n % i == 0:
            result -= result/i
        while n % i == 0:
            n /= i
        i += 1
    if n > 1:
        result -= result/n
    return result
print "Crunching Phi's\n"
max_n = 0
max_n_div_phi = float(0)
for n in range(2,1000001):
    phi = totient(n)
    n_div_phi = float(n/phi)
    if n_div_phi > max_n_div_phi:
        max_n = n
        max_n_div_phi = n_div_phi
        print "n:"+str(n)+" phi:"+str(phi)+" n/phi:"+str(n_div_phi)
    if n % 100000 == 0:
        print "Testing "+str(n)+"s"
print "\nAnswer: Value of n with max_n_div_phi is "+str(max_n)