Friday, 22 March 2013

Fast Euler Totient Function in C++

Euler Totient Function or phi(n) gives the No. of numbers less than or equal to n which are relatively prime to n.Some Postulates of Euler totient Function-
phi(p) = p-1 ; if p is a prime ;
phi(p^n) = (p-1) * p^(n-1) ; still if p is a prime..
phi(a*b) = phi(a) * phi(b) whenever gcd(a,b) = 1
 Read ETF in detail from http://en.wikipedia.org/wiki/Euler%27s_totient_function

Now an easy but slow implementation will be to prime factorize N and then apply phi(p^n) formula to all the prime factors in prime factorization of N.
i.e
Let N = p1^e1 * p2^e2 *..... pn^en
Then
phi(N) = phi(p1^e1) * phi(p2^e2) * ... phi(pn^en)
phi(N) = ((p1 - 1) * (p1 ^ e1 - 1)) * ((p2 - 1) * (p2 ^ (e2 - 1))) * ((pn - 1) * (pn ^ (en - 1)))

Implementation of this one is -
 
#include<iostream>
#include<cmath>
using namespace std;

long long int phi(long long x)
 {
   long long int ret = 1,i,pow;
   for (i = 2; x != 1; i++) 
   {
     pow = 1;
     if(i>sqrt(x))break;
     while (!(x%i)) 
     {
       x /= i;
       pow *= i;
     }
     ret *= (pow - (pow/i));
    }
    if(x!=1)ret*=(x-1);
    return ret;
}

int main()
{
    int t;
    cin>>t;
    while(t--)
    {
           int n;
           cin>>n;
           cout<<phi(n)<<endl;   
    }
    return 0;
} 
 
Complexity of this Solution is O(T * sqrt(N)) as T can be maximum 10^3 and N <= 10^6
so overall complexity of this solution is O(10^6) which will not time out...
But there exists many better and faster implementation of the Euler Totient Function by using sieve

I will discuss two of them -
First implementation -
In first implementation using sieve we will keep one of the prime factors of the numbers from 1 to N in the array..
i.e prime[10] = 5, prime[14] = 7, prime[31] = 31 (prime array will contain one of the prime factor of the Number.)
Now when we will get a value N, we can factorize it in O(Number of prime factors)..
Now for a range [1,10^6] the maximum value Number of prime factors can get is <20 (as 2^20  = 10^6)
So after performing the sieve we can compute phi[n] in O( T * number of prime factors of N ).

Implementation of the above sieve -

#include<iostream>
#define N 1000000
#include<cmath>
#include<cstdio>
using namespace std;
int prime[N+1];
 
void sieve(){
    prime[0]=1;
    prime[1]=1;
    for(long long i=2;i<=N;i++){
        if(!prime[i]){
                        prime[i]=i;
            for(long long j=i;j*i<=N;j++){
                prime[j*i]=i;
            }
        }
    }
}

long long int primeFactorize(int n){
    long long int yy=1;
    for( ; n > 1 ; ){
        int p = prime[n] , e = 0 ;
        for( ; n % p == 0 ; n /= p , e++ ) ;  
        yy*= ((p-1)*(pow((double)p,e-1)));
    }
    return yy;
}


int main()
{
   sieve();
   int t;
   scanf("%d",&t);
   while(t--)
   {
    int n;
    scanf("%d",&n);
    long long int ll = primeFactorize(n);
    printf("%lld\n",ll);              
   } 
    
   return 0;
}
 
In the second implementation we are going to precompute all the values of phi[n].
It must be implemented when the number of testcases are high. Because after precomputation, the time required to answer the query becomes O(1).
But in precomputation the time required is O(nloglogn).
In this sieve we again perform standard sieve by combining the result of phi of earlier values relatively prime to current value.. and updating phi with the primes stored in an array..

Here is the implementation of the second Sieve -
 
#include<iostream>
#include<vector>
#include<cmath>
#include<cstdio>
using namespace std;

#define N 1000000
#define MAXN 1000000
int phi[MAXN + 1], prime[MAXN/10], sz=0;
vector<bool> mark(MAXN + 1);

int main()
{
    
   phi[1] = 1; 
   for (int i = 2; i <= MAXN; i++ ){
    if(!mark[i]){
        phi[i] = i-1;
        prime[sz++]= i;
    }
    for (int j=0; j<sz && prime[j]*i <= MAXN; j++ ){
        mark[prime[j]*i]=1;
        if(i%prime[j]==0){
            int ll = 0;int xx = i;
            while(xx%prime[j]==0)
            {
                           xx/=prime[j];
                           ll++;         
                       }
            int mm = 1;
            for(int k=0;k<ll;k++)mm*=prime[j];
            phi[i*prime[j]] = phi[xx]*mm*(prime[j]-1);
            break;
    }
        else phi[i*prime[j]] = phi[i]*(prime[j]-1 );
    }
}
   int t;
   scanf("%d",&t);
   while(t--)
   {
    int n;
    scanf("%d",&n);
    printf("%d\n",phi[n]);             
   } 
   return 0;    
} 
 
All of this will get more clear by reading
1. Multiplicative Function
2. Sieve of Erasthothenes

Questions to Try -
https://www.spoj.com/problems/ETF/
https://www.spoj.com/problems/DCEPCA03/