Skip to main content

NT Part 2: Generating Primes, Prime Test, Prime Factorization

 

Generating primes fast is very important in some problems.
Let's cut to the chase and introduce Eratosthenes's Sieve.

The main idea is the following. Suppose we want to find all primes between 2 and 50.
Iterate from 2 to 50. We start with 2. Since it is not checked, it is a prime number. Now check all numbers that are multiple of $2$ except 2. Now we move on, to number 3. It's not checked, so it is a prime number.
Now check all numbers that are multiple of $3$, except 3. Now move on to 4. We see that this is checked - this is a multiple of 2! So 4 is not a prime. We continue doing this.

Here's the implementation.

#include <stdio.h>
int primechk[21000];
 
void preprocess(void)
{
    int i, j;
    for(i=2 ; i<=20000 ; i++)
    {
        primechk[i]=1;
    }
    for(i=2 ; i<=20000 ; i++)
    {
        if(primechk[i]==1)
        {
            for(j=2 ; i*j<=20000 ; j++)
            {
                primechk[i*j]=0;
            }
        }
    }
}
 
int main(void)
{
    preprocess();
    int i, cnt=0;
    for(i=2 ; i<=20000 ; i++)
    {
        if(primechk[i]==1)
        {
            cnt++;
            printf("%d is the %dth prime!\n",i,cnt);
        }
    }
}


Okay, so what is the time complexity? To get all primes in the interval $[1,n]$, the TC of this algorithm is $O(n \log \log n)$
Very very fast. To prove this, notice that the number of iterations are something like $O(n)+\sum_{p \le n} \frac{n}{p}$ where $p$ is a prime.

Well, Merten's Second Theorem states that $\sum_{p \le n} \frac{1}{p} = \log \log n + O(1)$ (natural logarithm, by the way) so this prove the TC.

Now we know how to generate prime numbers fast. How about primality testing?

Naive solutions first. Given an integer $n$, we can check numbers up to $\sqrt{n}$ to find if a number divides $n$. If there is such a number, $n$ is composite. If not, $n$ is a prime. This gives the solution in $O(\sqrt{n})$.

Here's a "solution" in $O(c\ln n)$ using the fast exponentiation we will talk about in the next section. It's called Miller-Rabin Primality Test.
I'll introduce the deterministic version. We probably (hopefully) won't see stuff like $n> 4 \cdot 10^9$ in contests.

Here's the sketch of the algorithm. Choose some set of $a$. We will run the algorithm with different $a$s, and the more $a$s we run this algorithm with, the more accurate this primality test is going to be.

Decompose $n-1$ as $2^s \cdot d$. Then check if the following holds.
$$a^d \not\equiv 1 \pmod{n} \text{ and }a^{2^rd} \not\equiv -1 \pmod{n} \text{ for all } r \in [0,s-1]$$If there is an $a$ that satisfies this, $n$ is composite. If not, $n$ is a prime.

For $n<4.7 \cdot 10^9$, we can just check for $a=2,7,61$ and be sure about it.
For $n<2^{64}$, we can check for $a=2,3,5,7,11,13,17,19,23,29,31,37$ and be confident.

Now let us look at prime factorization of numbers.

Assume that we generated prime numbers using the Eratosthenes's Sieve.
If $n$ is in the "prime-generated" range, we can actually do prime factorization in $O(\log n)$.
Make another array. While we are doing the Sieve, for composite numbers, put "the first prime that verified that this number is composite" and for prime numbers, put itself. This is easy to implement.

Then we can start with $n$, and continue to divide the prime numbers in the array.

#include <stdio.h>
int primechk[21000];
int fprime[21000]; 
void preprocess(void)
{
    int i, j;
    for(i=2 ; i<=20000 ; i++)
    {
        primechk[i]=1;
    }
    for(i=2 ; i<=20000 ; i++)
    {
        if(primechk[i]==1)
        {
            fprime[i]=i;
            for(j=2 ; i*j<=20000 ; j++)
            {
                primechk[i*j]=0;
                if(fprime[i*j]==0)
                {
                    fprime[i*j]=i;
                }
            }
        }
    }
}
 
int main(void)
{
    preprocess();
    int n;
    scanf("%d",&n);
    while(n!=1)
    {
        printf("%d\n",fprime[n]);
        n=n/fprime[n];
    }
}


If not, we can just check through all primes less than $\sqrt{n}$ and divide by those primes until we can't.
If all the primes multiply up to $n$, we are done. If not, there is exactly one prime more than $\sqrt{n}$ that divides $n$.

Another algorithm involving prime factorization is Pollard's rho algorithm - since the pseudocode is simple, I'll leave you the wikipedia link. https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm#Algorithm

Comments

Post a Comment

Popular posts from this blog

NT Part 1:GCD, LCM, Euclidean Algorithm

  The definitions of GCD and LCM are well-known, (and taught in middle school I think) I will skip the definitions. Also, since  ,  calculating GCD is equivalent to calculating LCM. Now, how do we calculate the GCD of two numbers? A naive solution would be iterating over all positive integers no more than  . This will get the GCD in  ,  very very slow. We can calculate the GCD of   in   using Euclidean Algorithm. This algorithm uses the easy-to-prove fact  ,  where   is the remainder when   is divided by  ,  or just a%b. We can now use the following code. #include <iostream> using namespace std ; int gcd ( int u, int v ) { return u % v == 0 ? v : gcd ( v,u % v ) ; } int main ( void ) { int x, y ; cin >> x >> y ; cout << gcd ( x,y ) ; } How do we prove that this algorithm is  ?  Well, let's suppose that we started with  . Then...
Post a Comment
Read more

NT Part 4: Mobius Function Part 1. The Introduction

  What is mobius function? This function is notated as  .  This function has lots of definitions. However, the main definition is the following.  is   if   or   is square-free and has even number of prime divisors.  is   if   is square-free and has odd number of prime divisors.  is   if   is not square-free.  also has a lot of interesting properties that make   so important.  for all  ,  and   is multiplicative. (A function   is multiplicative if   for all  . ) How do we calculate   fast? Again, we can tweak the Eratosthenes's Sieve a little bit. void preprocess ( void ) { int i, j ; for ( i = 1 ; i <= 111100 ; i ++ ) { mu [ i ] = 1 ; primechk [ i ] = 1 ; } primechk [ 1 ] = 0 ; for ( i = 2 ; i <= 111100 ; i ++ ) { if ( primechk [ i ] == ...
Post a Comment
Read more

NT Part 5: Mobius Function Part 2. Mobius Inversion, Dirichlet Convolution

  Theorem.  Let   be arithmetic function such that  .  Then,  . Proof. using   for all   and  . The Dirichlet Convolution of two function   is the following. Denote   as the constant function,  .  is a function which satisfies   and   for  .  is the identity function,  The mobius inversion formula and the basic property of mobius functions give  ,  and  Okay. So what? How can we use this formula to  improve our lives  solve problems? We can, change our "sum" or our answer using mobius inversion and mobius functions to calculate them fast. I will give 3 examples which changes the desired sum by using number theory (like mobius function) to calculate them fast. Example Problems Problem 1.  http://www.spoj.com/problems/LCMSUM/ Let   be the answer when  . Precalculate  ,  and add   to every  . ...
Post a Comment
Read more