PraveenDhinwa's blog

By PraveenDhinwa, 11 years ago, In English

The blog is not yet complete. I just mistakenly clicked post button. So it is still in writing, I have to add formatting to it. Sorry for inconvience

Structure of the article. First I will explain what does sieve mean. Then I will give you some examples with corresponding java codes and finally some exercises :)

According to http://www.thefreedictionary.com/sieve , "sieve" means A utensil of wire mesh or closely perforated metal, used for straining, sifting, ricing, or puréeing. Similar to this definition sieve is a method of doing stuff, where you keep rejecting the things that are useless and reducing the space that you are currently looking at.

So much of thing in air, Let us now take an example.

  1. Finding primes upto N. You have to print all primes upto N.

Method1 : For all the numbers i from 1 to N, check if i is prime or not. If it is a prime, then print it.

_Subproblem_ :

        Checking whether a number K is prime.

  _Solution_ :

         1. For all numbers i from 2 to K-1, check if K is divisible by i (as every number is divisible by 1 and itself). If yes, then not a prime else the number is a prime.

            Complexity of this solution : O(K)

         2. Note that we do not need to check upto K-1, instead we  can very well check upto sqrt(K).

          Proof: Let us say a number K = a * b. Note that atleast one of a and b <= sqrt(K) otherwise product of them would exceed K. So check just upto sqrt(K).

         3. Either use some probabilistic method for finding whether a number is prime or not.
 More on this later. For now see link

                             http://en.wikipedia.org/wiki/Primality_test
                         http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=primalityTesting

Method 2 :

Now here comes the idea of sieve. So initially assume that all numbers are prime. Then you try to sieve/refine your search range by not looking at the entire numbers but at the reduced space. eg. When you find all the numbers which are divisible by 2, You do not look into those again, as they are already not prime. So those numbers are sieved out. Now try for 3,4, upto n.

  In other terms, You first try all the numbers which are divisible are 2 (except 2 itself),

Note that all those wont be primes. So you can remove those out of your consideration now. Now try the same for 3,4,...... N. Finally you will end up with only prime numbers.

For understanding the actual thing going on, see the code.

 So the code basically sets all the numbers upto N to be prime. Then for every number that is still prime, we set all of its multiples upto N to be non-prime.

import java.io.*; import java.util.*; import java.math.*;

public class Main { static boolean[] isPrime;

public static void sieveOptimized(int N) {
    isPrime = new boolean[N + 1];

    for (int i = 2; i <= N; i++)
        isPrime[i] = true;
    for (int i = 2; i * i <= N; i++) {
            if (isPrime[i]) {
            // For further optimization, You can do instead of j += i, j += (2 * i).
            // Proof is left to reader :)
            for (int j = i * i; j <= N; j += i) 
                isPrime[j] = false;
        }
    }
}


public static void sieve(int N) {
    isPrime = new boolean[N + 1];

    for (int i = 2; i <= N; i++)
        isPrime[i] = true;
    for (int i = 2; i <= N; i++) {
        if (isPrime[i]) {
            for (int j = i + i; j <= N; j += i) 
                isPrime[j] = false;
        }
    }
}

public static void main(String[] args) {
    Scanner sc = new Scanner(System.in);
    int limit = sc.nextInt();

    // call the sieve method
    sieveOptimized(limit);

    for (int i = 1; i <= limit; i++) {
        if (isPrime[i]) 
            System.out.printf("%d ",i);
    }
}

}

Now comes the complexity part: Complexity of this code is O(n * logn).

 _Proof_ : (This proof comes a lot of times in various algorithms, So pay attention).

 For all numbers i going from 2 to n, you need to check all the multiples of i. Note that number of multiples of i upto n are n / i. Hence Expression for the complexity will be written as n / 2 + n / 3 + .. + 1. Take n common out of expression. Now it becomes n * (1 / 2 + ..... + 1/n). 

 Now as the expression is definitely greater than n. So adding an n to the expression won't have any effect on the complexity, So add n to the expression and now it becomes n * (1 + 1 / 2 + ... + 1/ n).  The expression (1 + 1 / 2 + 1 / 3 + ... 1 / n) is harmonic sum and it's bounded by ln(n). Hence overall complexity is O(n * logn)

_ Proof of harmonic sum:_ __ A simple Proof: Let us integrate 1 / x from 1 to n. (Note that we are doing integration, which means sum of area under the curve 1/x, which is greater than (1 + 1 / 2 + ... + 1 / n). Value of the integral can be found easily. In fact integration of 1/x dx is ln(x).

  1. Finding Sum of divisors of numbers upto N.

    Now you have to find sum of divisors of all numbers upto N. Here we are not just considering proper divisors(numbers other 1 and itself), we are considering all the divisors. Here you can do something like this.

    Here let us say divisorSum[i] denotes sum of divisors of i. Intially value of divisorSum[i] is equal to zero. Then for all the numbers i from 1 to n, We check all the multiples of i (let us say j) and add i to divisorSum[j].

    In other words, Start from 1 and for all the numbers which are multiples of 1, increment their sumDiviors by 1.

    Now do the same for 2,3, ... N. Note that for a number i, you are doing this adding operation upto N/i times. So the complexity calculation is same as before.

Look the beautifully commented code.

import java.io.*; import java.util.*; import java.math.*;

public class Main { public static void main(String[] args) { Scanner sc = new Scanner(System.in); int T = sc.nextInt();

while (T-- > 0) {
        int n = sc.nextInt();
        int divisorSum[] = new int [n + 1];

        // For every number i, You know that 2*i,3*i,4*i   upto k*i such that k*i<=n, will have i 
        // as one of it's divisors, so add that to divisorSum[j]

        for (int i = 1; i <= n; i++) {
            for (int j = i; j <= n; j += i) {
                divisorSum[j] += i;
            }
        }

        // Complexity of this code is O(n * logn)
        // Proof: Expression for the complexity can be written as n / 1 + n / 2 + ... + n / n
        // take n common
        // n * (1 + 1 / 2 + ..... + 1/n)
        // (1 + 1 / 2 + 1 / 3 + ... 1 / n) is harmonic sum and it's bounded by logn.
        // A simple Proof: Let us integrate 1 / x from 1 to n. 
        // (Note that we are doing integration, which means sum of area under the curve 1/x
        // which is greater than (1 + 1 / 2 + ... + 1 / n)
        // value of integration can be found easily
        // as integration of 1/x dx is ln(x)

        for (int i = 1; i <= n; i++) 
            System.out.printf("%d ", divisorSum[i]);

        System.out.printf("\n");
    }
}

}

  1. ** Finding No of divisors of numbers upto N.**

    This is also same as the previous example. Here instead of the storing sum in the array, store the number of divisors and for every multiple of i (say j), In the previous example, you were adding value i to divisorSum[j] , Here just increment the count of noOfDivisior[j] by one.
    
    Code is very easy and hence omitted. Complexity is also same.
  2. Sieve for finding euler phi function.

    I will denote the euler phi function for a number n by phi(n). phi(n) is defined as follows. It is count of how many number from 1 to n are coprime(having gcd value 1) to n. For example phi(6) is 2 as 1,5 are coprime to 6.

    Few properties of phi function :

    a). phi(p) = p - 1. Where p is prime. All the numbers from 1 to p - 1 are coprime to p.
    
    b). phi(a * b) = phi(a) * phi(b) where a and b are coprime.
    
    c). phi(p^k) = p^k - p^(k - 1). Note that here ^ denotes power. Here all the numbers from 1 to p^k are coprime to p^k except all the multiples of p, which are exactly p^(k -1).

    Method for finding:

    1. Simple : For all numbers from 1 to n, check if it is coprime to n or not, If yes add that to your answer.
    1. Let us say your number is n, which can be denoted as p1^k1 * p2^k2 ..... *p_m*k_m. Note that here p1, p2.... pm are prime. Basically n is written in it's prime representation.
    Then phi(n) would be [ p1^k1 - (p1^(k1-1) ) ] * ....... [p_m^k_m - (p_m^(k_m-1) )] . The expression for n can also be written as p1^k1 * p2^k2 * .... * p_m^k_m * (1 - 1/p1) * (1 - 1/p2) .... * (1 - 1/pm).

which is equal to n * (1 — 1/p1) * (1 — 1/p2) .... * (1 — 1/p_m).

See the code for details.

import java.io.*; import java.util.*; import java.math.*;

public class Main {

public static boolean isPrime (int n) {
    if (n < 2)
        return false;
    for (int i = 2; i * i <= n; i++)
        if (n % i == 0)
            return false;
    return true;
}

public static int eulerPhiDirect (int n) {
    int result = n;
    for (int i = 2; i <= n; i++) {
        if (isPrime(i))
            result -= result / i;
    }

    return result;
}

public static int eulerPhi (int n) {
    int result = n;
    // think that it like this, initially all numbers have gcd with n to be equal to 1.
    // Hence value of result is n
    // according to formulla  n * (1 - 1/p1) * (1 - 1/p2) .... * (1 - 1/pm). We will be calculating value 
    // of the product upto i. that is n * (1 - 1/p1) * ... (1 - 1/p_i)
    // So let us take example of p1. value of result after one iteration will be n - n / p1, which is precisly
    // n * (1 - 1/p1). 
    // Similarily by induction hypthesis we can say finally the result will be as required.

    for (int i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            result -= result / i;

            // By using while loop here, we are ensuing that all the numbers i will be prime.
            // because for every i, all it's multiples are gets removed.
            while (n % i == 0) {
                n /= i;
            }
        }
    }

    if (n > 1) {
        result -= result / n;
    }

    return result;
}

public static void main(String[] args) {
    Scanner sc = new Scanner(System.in);
    int T = sc.nextInt();

    while (T-- > 0) {
        int n = sc.nextInt();

        // call the eulerPhiDirect or eulerPhi method. 
        // culerPhi is more faster as it does not take sqrt(n) for checking for prime

        int ans = eulerPhiDirect (n);
        System.out.println(ans);
    }
}

}

Now Let us calculate value of sieve of all numbers from 1 to N.

Let us say eulerPhi[i] be the value of phi(i). Assign initially all the values of eulerPhi[i] to be i. Then for every prime p, for all multiples of p, we will multiply value of eulerPhi[i] by (1 - 1/p) as per the formula. multiplying eulerPhi[i] by (1 - 1/p) is exactly equal to eulerPhi[i] -= (eulerPhi[i] / p).

import java.io.*; import java.util.*; import java.math.*;

public class Main {

public static boolean isPrime (int n) {
    if (n < 2)
        return false;
    for (int i = 2; i * i <= n; i++)
        if (n % i == 0)
            return false;
    return true;
}

private static int[] eulerPhi;

public static void eulerSieve (int N) {
    eulerPhi = new int[N + 1];

    // set initial value of phi(i) = i
    for (int i = 1; i <= N; i++)
        eulerPhi[i] = i;

    // for every prime i, do as described in blog.
    // Note that we are using isPrime(i) function that takes sqrt(n) time. You are advised to write 
    // a seperate sieve of finding primes as described by me.
    // which will reduce the compleixty of this to n * log(n)
    // Proof of this is similar to previous ones. Left to reader.

    for (int i = 1; i <= N; i++) {
        if (isPrime(i))
            for (int j = i; j <= N; j += i)
                eulerPhi[j] -= eulerPhi[j] / i;
    }
}

public static void main(String[] args) {
    Scanner sc = new Scanner(System.in);
    int N = sc.nextInt();

    eulerSieve(N);

    for (int i = 1; i <= N; i++)
        System.out.printf("%d ",eulerPhi[i]);
}

}

  1. Sieve for finding inverse modulo m.

    Inverse of a number a with respect to m is defined as b * a = 1 mod m. Then b is called inverse of a modulo m.

    So first of all do you realize that it is not necessary to exist the inverse? Example for a = 4, b = 6.

    Let us play around the expressions. a * b = 1 mod m can be written as a * b = 1 + m * k. which can be written as a * b — m * k = 1. If the left side has a common factor of both a and m, means gcd(a,m) != 1, then note that right side won't be divisible by that number, Hence no solution of the equation when a and m has gcd non zero. Hence inverse will exist when a and m have non zero gcd.

    Now solving a * b + m * (-k) = 1. write the same as a * b + m * k1 = 1.

    So let us try to find solution of a * b + k * m = 1. This k is not equal to previous k, infact it is -k. It is equal to k1.

    So let us try to solve generic problem a * x + b * y = 1. where a and b can also be negative and gcd(a, b) = 1.

    Let us try a simpler version of the same problem which is solving b * x1 + (a % b) * y = 1;

    Now try to relate these equations. a % b = a — [ a / b] * b. where [] denotes floor function. This is same as removing all multiples of b from a, which is exactly equal to a % b.

    Now the equation turns into b * x1 + (a — [a/b] * b) * y1 = 1

    which is a * y1 + b * (x1 — [a/b] * y1) = 1.

    So this is recursive version of a * x + b * y = 1, where x is replaced with y1 and y is replaced with x1 — [a/b] * y1.

    Things getting really complex. (Note that this method is exactly similar to finding gcd of two numbers). Seeing the code will help you to understand more.)

    Complexity of this code is same as gcd as it has exactly the same recurrence relation as of that. Time complexity of gcd(m, n) is log (max(m , n)). Hence we can consider time complexity of this method around O(logn).

import java.io.*; import java.util.*; import java.math.*;

class pair { public int x, y;

pair (int x,int y) {
    this.x = x;
    this.y = y;
}   

boolean isEquals (pair p) {
    if (this.x == p.x &amp;&amp; this.y == p.y) 
        return true;
    else 
        return false;
}

}

public class Main { public static int gcd (int a, int b) { if (b == 0) return a;
return gcd (b, a % b); }

public static pair solve (int a,int b) {
    if (b == 0) {
        // a * x + b * y = 1
        // here b = 0
        // hence a * x = 1
        // if a is not 1, then error else x = 1 and y = 0
        // Note that error wont be here, we will always find a which is not 1
        // as error case is already handle in solveThis function
        return new pair (1, 0);
    } else {
        // do the recursive call
        pair p = solve (b, a % b);
        int x1 = p.x;
        int y1 = p.y;

        int x = y1;
        int y = x1 - (a / b) * y1;

        return new pair (x, y);
    }
}

public static pair solveThis (int a, int b) {
    if (gcd (a, b) != 1)
        // (-1, -1) corresponds to error, that means no inverse exists in this case
        return new pair (-1, -1);
    else 
        return solve (a, b);
}

public static int modpow (long a, long b, long c) {
    long res = 1;
    while (b >= 0) {
        if (b % 2 == 1) {
            res = (res * a) % c;
        }
        a = (a * a) % c;
        b >>= 1;
    }

    return (int) res;
}

public static int findInverseModuloPrime (int a, int p) {
    return modpow (a, p - 2, p);
}

public static void main(String[] args) {
    Scanner sc = new Scanner(System.in);

    int a = sc.nextInt();
    int m = sc.nextInt();

    pair p = solveThis (a, m);

    if (p.isEquals(new pair(-1, -1)))
        System.out.printf("Error, inverse does not exist");
    else 
        System.out.printf("%d %d\n", p.x, p.y);
}

}

Now I will tell you another easier method . Generalized version of Fermat's theorem says that a ^ phi(m) = 1 mod m where gcd(a, m) = 1. Fir finding inverse a ^ (phi(m) — 1) = 1 / a (mod m) = (inverse(a)) mod m.

Hence for finding inverse(a) mod m, You can just find a ^ (phi(m) — 1) by modular exponention method. In case of m being prime, As phi(m) = m — 1. So just find a ^ (m — 2) % m. This is what precisely computed by modpow function in the above function.

Complexity :

Now Sieve for finding inverse modulo m:

You have to find inverse(i) mod m for i ranging from 1 to n. as complexity of modpow is log (n). and we are doing this for n numbers. Hence total complexity will be O(n * logn). Here I am going to describe a better method using sieve

Just use the identitiy inverse(i) = (-(m/i) * inverse(m % i) ) % m + m;

Proof: It is good one. I do not want reveal it now. Try yourself. If you come up with it, post it in the comments, I would check whether it is correct or not?

import java.io.*; import java.util.*; import java.math.*;

public class Main { public static void main(String[] args) { Scanner sc = new Scanner(System.in);

int n = sc.nextInt(); // denotes the range upto which you want to find the value of modInverse[i]
    int m = sc.nextInt();

    int modInverse[] = new int[n + 1];

    modInverse[1] = 1; // this is you know 1 * 1 mod m = 1

    for (int i = 2; i &lt;= n; i++) {
        modInverse[i] = (-(m/i) * modInverse[m % i]) % m + m;
    }

    for (int i = 2; i &lt;= n; i++) 
        System.out.printf("%d ", modInverse[i]);
}

}

Finally, Here are some few exercises for you to try

1. http://www.spoj.com/problems/GCDEX/

 I will keep adding the list if I found some problems.
  • Vote: I like it
  • +18
  • Vote: I do not like it

| Write comment?
»
11 years ago, # |
  Vote: I like it 0 Vote: I do not like it

very well written ..

»
11 years ago, # |
Rev. 3   Vote: I like it 0 Vote: I do not like it

I am facing problem in editing the post in the blog. I have to put br tag in each line to make sure that the line endings are shown as they are written by me. Could anybody suggest me some better way of putting line endings in codeforces blogs?

  • »
    »
    11 years ago, # ^ |
      Vote: I like it 0 Vote: I do not like it

    Leave a blank line to start a new paragraph.

  • »
    »
    11 years ago, # ^ |
    Rev. 2   Vote: I like it +3 Vote: I do not like it

    Perhaps article on markdown syntax may help you to use some other of its features to make your post more structured. For example, you probably find it useful to insert some

    Titles

    Subtitles

    Paragraph headers

    and outline the inline formulas like S = 3.14 * R ^ 2

    At the same time fixing some fragments where unintended styles were applied.

    • »
      »
      »
      11 years ago, # ^ |
        Vote: I like it 0 Vote: I do not like it

      Code part is not shown very well here. I tried to embed gist but it did not work :(

      Could you suggest some other alternative ?

      • »
        »
        »
        »
        11 years ago, # ^ |
        Rev. 2   Vote: I like it +3 Vote: I do not like it

        Code part is shown well enough, though it looks you tried to represent it in some improper way. As it is told in the mentioned article, code block is defined as any sequence of lines having 4 or more spaces at the beginning of each. Also, you should not use escape sequences (like <) in code blocks — they would be processed automatically:

        if (some < where(i) && !quit()) {
            heat(Gear.getEngine());
        }

        (and by the way I am not sure your nice article really need as much code snippets)

        UPD: Moreove, this post states that improved version of markdown is used, which allows syntax highliting.

      • »
        »
        »
        »
        10 years ago, # ^ |
          Vote: I like it 0 Vote: I do not like it

        To write some code:


        ~~~~~ #include <stdio.h> ... ~~~~~

        To write some little code


        `eulerPhi[i]`
»
11 years ago, # |
  Vote: I like it 0 Vote: I do not like it

Big thanks , very helpful :)

»
11 years ago, # |
  Vote: I like it 0 Vote: I do not like it

Hi, Thanks for your post it is helpful for me. I am having a doubt about calculating the number of numbers whose numbers of divisors are prime withing a range A — > B The sieve is the key but i cannot handle with the constrain that A ≤ B ≤ 10^12

How can i handle a good implementation using this property 0 ≤ B-A ≤ 10^6

  • »
    »
    11 years ago, # ^ |
      Vote: I like it +1 Vote: I do not like it

    Wait until CodeChef November challenge ends ;-)

  • »
    »
    11 years ago, # ^ |
      Vote: I like it 0 Vote: I do not like it

    I will not give a complete implementation right because there is a problem on CodeChef November that requires this idea. But I can give you a hint: in usual sieve you mark all multiples of prime i, but in this case you need just those multiples which lie in the range [A; B], think how you can fastly find the closest multiple of i which is <= A.

    I hope this helps.

    Cheers,

»
10 years ago, # |
  Vote: I like it +1 Vote: I do not like it

Very Well written.Too helpful.Interactive, informative ,optimized and bang on to the spot!!! In short excellent!

»
10 years ago, # |
  Vote: I like it 0 Vote: I do not like it

This link can be useful.

»
10 years ago, # |
Rev. 3   Vote: I like it +38 Vote: I do not like it

The pure euler-phi sieve:

int[] phi = new int[n];
for (int i = 0; i < n; i++) {
	phi[i] = i;	
}
for (int i = 1; i < n; i++) {
	for (int j = 2 * i; j < n; j += i){
		phi[j] -= phi[i];
	}
}

Proof:

And you can apply the same idea to the Mobius function.

int[] mu = new int[n];
mu[1] = 1;
for (int i = 1; i < n; i++) {
	for (int j = 2 * i; j < n; j += i) {
		mu[j] -= mu[i];
	}
}
  • »
    »
    10 years ago, # ^ |
      Vote: I like it +5 Vote: I do not like it

    Wow, nice beautiful trick. Not fastest asymptotically but still beautiful and short enough.

»
10 years ago, # |
Rev. 2   Vote: I like it 0 Vote: I do not like it

Hints :

  • With your code : you can use Block to view beautiful and clear.

Example :

public static int modpow (long a, long b, long c) {
    long res = 1;
    while (b >= 0) {
        if (b % 2 == 1) {
            res = (res * a) % c;
        }
        a = (a * a) % c;
        b >>= 1;
    }

    return (int) res;
}
  • If you want to express the important equations or note, you can use Inline.

    Example : a ^ (m — 2) % m. If you want to more beautiful,

    you can use insert Picture [Ctrl+P] :

»
10 years ago, # |
Rev. 2   Vote: I like it +12 Vote: I do not like it

PraveenDhinwa
the explanation for inverse modulo :
a * (m / a) + m % a = m
(a * (m / a) + m % a) mod m = m mod m, or
(a * (m / a) + m % a) mod m = 0, or
(- (m % a)) mod m = (a * (m / a)) mod m.
Dividing both sides by (a * (m % a)), we get
– inverse(a) mod m = ((m/a) * inverse(m % a)) mod m
inverse(a) mod m = (- (m/a) * inverse(m % a)) mod m