Code-Challenge: Farey sequence (II)

Question

Challenge

In this task you would be given an integer N (less than 10^6), output the number of terms in the Farey sequence of order N.

The input N < 10⁶ is given in a single line, the inputs are terminated by EOF.

The challenge is to implement the fastest solution so that it can process a maximum of 10⁶-1 values as fast as possible.

Input

7
10
999
9999
79000
123000

Output

19
33
303793
30393487
1897046599
4598679951

Constraints

• You can use any language of your choice
• Take care about your fingers, do not use more than 4096 bytes of code.
• Fastest solution, based on the running time measured for N = 10⁶-1.

Answer 1

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)).

Answer 2

J, 0.2s, short code

It feels weird to write ungolfed J, but it actually is the language for the task, for once.

input =: ". ;. _2 stdin ''
phi =: 5 & p:
max =: >. / input
acc_totient =: +/ \ phi >: i. max
output =: >: (<: input) { acc_totient
echo output

Reads input, constructs accumulated sum of totients up to the max, and reads output from there.

Answer 3

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 ...

Answer 4

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

Answer 5

VB.NET, 0.11 s for 10^6

This is basically a totient computation using prime sieve. The time is measured 10 times and the average is computed.

Module Module1

Sub Main()
    Dim sw As New Stopwatch
    Dim M As Integer = Console.ReadLine()
    sw.Start()

    Dim is_prime(M) As Integer
    Dim phi(M) As Integer

    For i = 0 To M
        phi(i) = i
        is_prime(i) = True
    Next

    For i = 2 To M
        If is_prime(i) = True Then
            phi(i) = i - 1
            Dim j As Integer = 2 * i
            Do Until j > M
                is_prime(j) = False
                phi(j) = phi(j) - (phi(j) / i)
                j += i
            Loop
        End If
    Next

    Dim result As Long = 1
    For i = 1 To M
        result += phi(i)
    Next

    Console.WriteLine(result)
    sw.Stop()
    Console.WriteLine("Time Elasped : " & sw.ElapsedMilliseconds)
    Dim a = Console.ReadLine()
End Sub

End Module

Answer 6

Rust

Rust port of @fR0DDY's C answer.

use std::io::{self};

fn main() {
    let mut is_prime = [true; 1001];
    let mut a = [0i64; 1000001];

    // Sieve of Eratosthenes for Prime
    for d in 2..1001 {
        if is_prime[d] {
            for n in ((d * d)..1001).step_by(d) {
                is_prime[n] = false;
                a[n] = d as i64;
            }
            for n in ((d * d)..=1000000).step_by(d) {
                a[n] = d as i64;
            }
        }
    }

    // Euler's Totient Function computation
    a[1] = 1;
    for n in 2..=1000000 {
        if a[n] == 0 {
            a[n] = n as i64 - 1;
        } else {
            let p = a[n];
            let mut k = n / p as usize;
            let mut e = 1;
            while k % p as usize == 0 {
                k /= p as usize;
                e *= p as i64;
            }
            a[n] = a[k] * e * (p - 1);
        }
    }

    // Sum for Farey sequence length
    a[1] = 2;
    for n in 2..=1000000 {
        a[n] += a[n - 1];
    }

    // Reading input and printing results
    let stdin = io::stdin();
    let mut buffer = String::new();
    while stdin.read_line(&mut buffer).unwrap() > 0 {
        let i = buffer.trim().parse::<usize>().unwrap();
        println!("{}", a[i]);
        buffer.clear();
    }
}

Attempt This Online!