Code-Challenge: Farey sequence (II)

C - 0.1 Secs on Ideone

http://www.ideone.com/E3S2t

Explanation included in code.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

long long A[1000001]={0};
int main()
{ 
    int isprime[1001],d,n,k,i,e,p; 

    for (n=2;n<1001;n++)
        isprime[n]=1;

    //Sieve of Eratosthenes for Prime
    //Storing the smallest prime which divides n.
    //If A[n]=0 it means it is prime number.
    for(d=2;d<1001;d++)
    {
        if(isprime[d])
        {
            for (n=d*d;n<1001;n+=d)
            {
                isprime[n]=0;
                A[n]=d;
            }
            for (;n<=1000000;n+=d)
                A[n]=d;
        }
    }

    //Uses multiplicative property of
    //Euler Totient Function
    //Phi(n)=Phi(p1^k1)*Phi(p2^k2)...*Phi(pr^kr)
    //and Phi(pi^ki)=(pi-1)*(pi^(ki-1))
    A[1]=1;
    for(n=2;n<=1000000;n++)
    {
        if (A[n]==0)
           A[n]=n-1;
        else
        {
           p=A[n],k=n/p,e=1;
           while (k%p==0)
                 k/=p,e*=p;
           A[n]=A[k]*e*(p-1);
        }
    }
    //Number of terms in Farey Series
    //|F(n)| = 1 + Sigma(i,1,n,phi(i))
    A[1]=2;
    for(n=2;n<=1000000;n++)
        A[n]+=A[n-1];


    while (~scanf("%d",&i))
        printf("%lld\n",A[i]);
    return 0;
}

A little more explanation:

For all numbers we get a prime factor of it from the sieve(or 0 if it is a prime). Next we use the fact that ETF is multiplicative. That is if m and n are coprime then phi(m*n)=phi(m)*phi(n). Here we take out the multiple of prime factor out and hence the left part and the multiple part are co-prime. We already have the ETF for the left part since it is either smaller then current value or equal to 1. We only need to calculate the ETF for the multiple which we calculate using the formula phi(pi^ki)=(pi-1)*(pi^(ki-1)).

C - Less than a second for 999999

#include "stdio.h"
#include "time.h"
#define i64 unsigned long long
#define i32 unsigned int
#define PLIM 1100
#define FLIM 1000000

i32 primes[PLIM];
i64 f[FLIM];

int main(int argc, char *argv[]){
    i64 start=clock();
    i32 primesindex=0;
    f[0]=1;
    f[1]=2;
    i32 a,b,c;
    i32 notprime;
    //Make a list of primes
    for(a=2;a<PLIM;a++){
        notprime=0;
        for(b=0;b<primesindex;b++){
            if(!(a%primes[b])){
                notprime=1;
            }
        }
        if(!notprime){
            primes[primesindex]=a;
            primesindex++;
        }
    }
    i32 count,divided;
    i32 nextjump=4;
    i32 modlimit=2;
    i32 invalue;
    a=2;
    for(c=1;c<argc;c++){
        invalue=atoi(argv[c]);
        if(invalue<FLIM && invalue){
            //For each number from a to invalue find the totient by prime factorization
            for(;a<=invalue;a++){
                count=a;
                divided=a;
                b=0;
                while(primes[b]<=modlimit){
                    if(!(divided%primes[b])){
                        divided/=primes[b];
                        //Adjust the count when a prime factor is found
                        count=(count/primes[b])*(primes[b]-1);
                        //Discard duplicate prime factors
                        while(!(divided%primes[b])){
                            divided/=primes[b];
                        }
                    }
                    b++;
                }
                //Adjust for the remaining prime, if one is there
                if(divided>1){
                    count=(count/divided)*(divided-1);
                }
                //Summarize with the previous totients
                f[a]=f[a-1]+(i64)count;
                //Adjust the limit for prime search if needed
                if(a==nextjump){
                    modlimit++;
                    nextjump=modlimit*modlimit;
                }
            }
            //Output result
            printf("%I64u\n",f[invalue]);
        }
    }
    i64 end=clock();
    //printf("Runtime: %I64u",end-start);
    return 0;
}

This takes input from command line.

The computation is a simple sum of totients, it's only done once and only up to the biggest input.

Here is an Ideone version: http://www.ideone.com/jVbc0

C (0.2 sec in Ideone)

I thought of adding my approach too: (this is not as fast as Gaurav's)

 #include <stdio.h>
 #define N 1000000
 long long Phi[N+1];

 //Sieve ETF

int main(){
   int i,j;

 for(i=1;i<=N;i++)
    Phi[i]=i;

 for(i=2;i<=N;i++){
    if(Phi[i] == i)
       for(j=i;j<=N;j+=i)
           Phi[j] = (Phi[j]/i)*(i-1);
  }

int t,n;

 Phi[1]=2;
 for(i = 2; i < N; i++)
     Phi[i] += Phi[i-1];


for(;~scanf("%d",&n);)
    printf("%lld\n",Phi[n]);

return 0;
}

TESTING ...