Computing Eulers Totient Function

Question

I am trying to find an efficient way to compute Euler's totient function.

What is wrong with this code? It doesn't seem to be working.

def isPrime(a):
    return not ( a < 2 or any(a % i == 0 for i in range(2, int(a ** 0.5) + 1)))

def phi(n):
    y = 1
    for i in range(2,n+1):
        if isPrime(i) is True and n % i  == 0 is True:
            y = y * (1 - 1/i)
        else:
            continue
    return int(y)

Morehead

1/i does not do what you think it does - try it. – Ginnifer

use python3 instead of python2 :-) – Temperament

Or put from __future__ import division at the top of your code to enable float division in Python 2. – Fastidious

Project Euler huh? I'd compile a list of primes beforehand or at least cache the ones you've found. – Gynecology

Answer 1

Here's a much faster, working way, based on this description on Wikipedia:

Thus if n is a positive integer, then φ(n) is the number of integers k in the range 1 ≤ k ≤ n for which gcd(n, k) = 1.

I'm not saying this is the fastest or cleanest, but it works.

from math import gcd

def phi(n):
    amount = 0        
    for k in range(1, n + 1):
        if gcd(n, k) == 1:
            amount += 1
    return amount

Ginnifer

Did you mean range(1, n+1)? – Catching

These properties comes handy If p is a prime number, ϕ(p)=p−1 as gcd(p,q)=1where 1<=q<p If p is a prime number and k>0, ϕ(p^k)=p^k−p^(k−1) If a and b are relatively prime, ϕ(ab)=ϕ(a).ϕ(b) e-maxx-eng.github.io/algebra/phi-function.html – Inmate

Answer 2

You have three different problems...

1. y needs to be equal to n as initial value, not 1
2. As some have mentioned in the comments, don't use integer division
3. n % i == 0 is True isn't doing what you think because of Python chaining the comparisons! Even if n % i equals 0 then 0 == 0 is True BUT 0 is True is False! Use parens or just get rid of comparing to True since that isn't necessary anyway.

Fixing those problems,

def phi(n):
    y = n
    for i in range(2,n+1):
        if isPrime(i) and n % i == 0:
            y *= 1 - 1.0/i
    return int(y)

Anchovy

Answer 3

Calculating gcd for every pair in range is not efficient and does not scales. You don't need to iterate throught all the range, if n is not a prime you can check for prime factors up to its square root, refer to https://mcmap.net/q/73369/-why-do-we-check-up-to-the-square-root-of-a-number-to-determine-if-the-number-is-prime. We must then update phi for every prime by phi = phi*(1 - 1/prime).

def totatives(n):
    phi = int(n > 1 and n)
    for p in range(2, int(n ** .5) + 1):
        if not n % p:
            phi -= phi // p
            while not n % p:
                n //= p
    #if n is > 1 it means it is prime
    if n > 1: phi -= phi // n 
    return phi

Indistinct

This could be made a lot faster by reworking the loop so it stops at the square root of the current n value rather than the initial n value. – Wentz

At least it beats the explicit primality tests and gcd computations all the other existing answers use, though. – Wentz

Answer 4

I'm working on a cryptographic library in python and this is what i'm using. gcd() is Euclid's method for calculating greatest common divisor, and phi() is the totient function.

def gcd(a, b):
    while b:
        a, b=b, a%b
    return a
def phi(a):
    b=a-1
    c=0
    while b:
        if not gcd(a,b)-1:
            c+=1
        b-=1
    return c

Nicknack

Answer 5

Most implementations mentioned by other users rely on calling a gcd() or isPrime() function. In the case you are going to use the phi() function many times, it pays of to calculated these values before hand. A way of doing this is by using a so called sieve algorithm.

https://mcmap.net/q/902937/-a-fast-prime-number-sieve-in-python This answer on stackoverflow provides us with a fast way of finding all primes below a given number.

Oke, now we can replace isPrime() with a search in our array.

Now the actual phi function:

Wikipedia gives us a clear example: https://en.wikipedia.org/wiki/Euler%27s_totient_function#Example

phi(36) = phi(2^2 * 3^2) = 36 * (1- 1/2) * (1- 1/3) = 30 * 1/2 * 2/3 = 12

In words, this says that the distinct prime factors of 36 are 2 and 3; half of the thirty-six integers from 1 to 36 are divisible by 2, leaving eighteen; a third of those are divisible by 3, leaving twelve numbers that are coprime to 36. And indeed there are twelve positive integers that are coprime with 36 and lower than 36: 1, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, and 35.

TL;DR

With other words: We have to find all the prime factors of our number and then multiply these prime factors together using foreach prime_factor: n *= 1 - 1/prime_factor.

import math

MAX = 10**5

# CREDIT TO https://mcmap.net/q/902937/-a-fast-prime-number-sieve-in-python
def sieve_for_primes_to(n): 
    size = n//2
    sieve = [1]*size
    limit = int(n**0.5)
    for i in range(1,limit):
        if sieve[i]:
            val = 2*i+1
            tmp = ((size-1) - i)//val 
            sieve[i+val::val] = [0]*tmp
    return [2] + [i*2+1 for i, v in enumerate(sieve) if v and i>0]

PRIMES = sieve_for_primes_to(MAX)
print("Primes generated")


def phi(n):
    original_n = n
    prime_factors = []
    prime_index = 0
    while n > 1: # As long as there are more factors to be found
        p = PRIMES[prime_index]
        if (n % p == 0): # is this prime a factor?
            prime_factors.append(p)
            while math.ceil(n / p) == math.floor(n / p): # as long as we can devide our current number by this factor and it gives back a integer remove it
                n = n // p

        prime_index += 1

    for v in prime_factors: # Now we have the prime factors, we do the same calculation as wikipedia
        original_n *= 1 - (1/v)

    return int(original_n)

print(phi(36)) # = phi(2**2 * 3**2) = 36 * (1- 1/2) * (1- 1/3) = 36 * 1/2 * 2/3 = 12

Linders

Answer 6

It looks like you're trying to use Euler's product formula, but you're not calculating the number of primes which divide a. You're calculating the number of elements relatively prime to a.

In addition, since 1 and i are both integers, so is the division, in this case you always get 0.

Spahi

Answer 7

With regards to efficiency, I haven't noticed anyone mention that gcd(k,n)=gcd(n-k,n). Using this fact can save roughly half the work needed for the methods involving the use of the gcd. Just start the count with 2 (because 1/n and (n-1)/k will always be irreducible) and add 2 each time the gcd is one.

Whipple

You're not considering the work required to calculate the work required to compare n with 2*k, as well as the work required to subtract n-k. – Recoverable

Answer 8

Here is a shorter implementation of orlp's answer.

from math import gcd
def phi(n): return sum([gcd(n, k)==1 for k in range(1, n+1)])

As others have already mentioned it leaves room for performance optimization.

Lakeishalakeland

Answer 9

Actually to calculate phi(any number say n)
We use the Formula
where p are the prime factors of n.

So, you have few mistakes in your code:
1.y should be equal to n
2. For 1/i actually 1 and i both are integers so their evaluation will also be an integer,thus it will lead to wrong results.

Here is the code with required corrections.

def phi(n):
    y = n
    for i in range(2,n+1):
        if isPrime(i) and n % i  == 0 :
            y -= y/i
        else:
            continue
    return int(y)

Pneumonectomy

Answer 10

By definition, the Euler Totient function φ(.) is given by

The python implementation to compute φ(.) is pretty straightforward:

def Euler_totient_def(n): # φ(n)
    return sum([gcd(i, n) == 1 for i in range(1, n)]

def gcd(a, b): # complexity O(lg a)
    while b:
        a, b = b, a % b
    return a

Also, by the fundamental theorem of arithmetic and by the property of the Totient function φ(.), we have,

and the corresponding implementation:

def sieve_Eratosthenes(n):      # precompute primes <= n: complexity O(n lglg n)
    P = np.ones(n+1, dtype=int)
    P[0:2] = 0
    for e in range(2, n+1):
        P[range(2*e, n+1, e)] = 0
    return [i for i in range(n+1) if P[i]]  
    
def factorize(n, P=None):       # prime factorization of n with precomputed table P: complexity O(lg n)
    if P is None:
        P = sieve_Eratosthenes(n)
    i = 0
    f = set([])
    while n > 1:
        if n % P[i] == 0:
            n //= P[i]
            f.add(P[i])
        else:
            i += 1
    return f

def Euler_totient_prime_factors(n, P): # φ(n)
    f = factorize(n, P)
    return int(n*np.prod([(1-1/p) for p in f])) 

Now let's compare the time to compute φ(.) with both the above implementations:

import numpy as np
from time import time   

n = 1234571
start = time()
print(Euler_totient_prime_factors(n, sieve_Eratosthenes(n)))
# 1089792
# prime factorization of n: 13 x 23 x 4129
print(time() - start)
# 23.920003175735474
start = time()
print(Euler_totient_def(n))
# 1089792
print(time() - start)
# 5.4549994468688965

The first implementation is much faster, as can be seen above. The precomputed prime table comes handy when we want to compute φ(.) for a bulk of numbers instead of a single number, as shown below:

def Euler_totients_prime_factors(a, N): # compute φ for numbers in list a
    P = sieve_Eratosthenes(N)
    return [Euler_totient_prime_factors(n, P) for n in a]

def Euler_totients_def(a):  # compute φ for numbers in list a
    return [Euler_totient_def(n) for n in a]

Now let's compute φ(.) for a list of numbers:

a = [12345, 123456, 1234561, 1234563, 1234567, 1234569, 1234571]
start = time()
print(Euler_totients_prime_factors(a, max(a)))
# [6576, 41088, 1228500, 704880, 1224720, 705456, 1089792]
# prime factorization 12345 = 3 x 5 x 823
#                     123456 = 2^6 x 3 x 643
#                     1234561 = 211 x 5851
#                     1234563 = 3 x 11^2 x 19 x 179
#                     1234567 = 127 x 9721
#                     1234569 = 3 x 7 x 58789
#                     1234571 = 13 x 23 x 4129
print(time() - start)
# 13.774996757507324

Ac can be seen from above, now the prime factorization with precomputed prime table computes φ values faster

start = time()
print(Euler_totients_def(a))
# [6576, 41088, 1228500, 704880, 1224720, 705456, 1089792]
print(time() - start)
# 18.693108558654785

Tape