Tuesday, September 29, 2015

Modular Inverse from 1 to N

We already learned how to find Modular Inverse for a particular number in a previous post, "Modular Multiplicative Inverse". Today we will look into finding Modular Inverse in a bulk.


Problem

Given $N$ and $M$ ( $N < M$ and $M$ is prime ), find modular inverse of all numbers between $1$ to $N$ with respect to $M$.

Since $M$ is prime and $N$ is less than $M$, we can be sure that Modular Inverse exists for all numbers. Why? Cause prime numbers are coprime to all numbers less than them.

We will look into two methods. Later one is better than the first one.

$O(NlogM)$ Solution

Using Fermat's little theorem, we can easily find Modular Inverse for a particular number.

$A^{-1} \:\%\: M = bigmod(A,M-2,M)$, where $bigmod()$ is a function from the post "Repeated Squaring Method for Modular Exponentiation". The function has complexity of $O(logM)$. Since we are trying to find inverse for all numbers from $1$ to $N$, we can find them in $O(NlogM)$ complexity by running a loop.
int inv[SIZE]; ///inv[x] contains value of (x^-1 % m)
for ( int i = 1; i <= n; i++ ) {
    inv[i] = bigmod ( i, m - 2, m );
}
But it's possible to do better. 

$O(N)$ Solution

This solution is derived using some clever manipulation of Modular Arithmetic.

Suppose we are trying to find the modular inverse for a number $a$, $a < M$, with respect to $M$. Now divide $M$ by $a$. This will be the starting point. 

$M = Q \times a + r$, (where $Q$ is the quotient and $r$ is the remainder) 
$M = \lfloor \frac{M}{a} \rfloor \times a + (M \:\%\: a )$

Now take modulo $M$ on both sides.

$0 \equiv \lfloor \frac{M}{a} \rfloor \times a + (M \:\%\: a ) \:\:\:\text{(mod M )}$
$  (M \:\%\: a ) \equiv -\lfloor \frac{M}{a} \rfloor \times a \:\:\:\text{(mod M )}$

Now divide both side by $a \times ( M \:\%\: a )$.

$$\frac{M \:\%\: a}{a \times ( M \:\%\: a )} \equiv \frac{-  \lfloor \frac{M}{a} \rfloor \times a } { a \times ( M \:\%\: a ) } \:\:\:\text{(mod M)}$$
$$\therefore a^{-1} \equiv - \lfloor \frac{M}{a} \rfloor \times ( M \:\%\: a )^{-1} \:\:\:\text{(mod M)}$$

The formula establishes a recurrence relation. The formula says that, in order to find the modular inverse of $a$, we need to find the modular inverse of $b = M \:\%\: a$ first. 

Since $b = M \:\%\: a$, we can say that its value lies between $0$ and $a-1$. But, $a$ and $M$ are coprime. So $a$ will never fully divide $M$. Hence we can ignore the possibility that $b$ will be $0$. So possible values of $b$ is between $1$ and $a-1$.

Therefore, if we have all modular inverse from $1$ to $a-1$ already calculated, then we can find the modular inverse of $a$ in $O(1)$.

Code

We can now formulate our code.
int inv[SIZE];
inv[1] = 1;
for ( int i = 2; i <= n; i++ ) {
    inv[i] = (-(m/i) * inv[m%i] ) % m;
    inv[i] = inv[i] + m;
}
In line $2$, we set the base case. Modular inverse of $1$ is always $1$. Then we start calculating inverse from $2$ to $N$. When $i=2$, all modular inverse from $1$ to $i-1=1$ is already calculated in array $inv[]$. So we can calculate it in $O(1)$ using the formula above at line $4$.

At line $5$, we make sure the modular inverse is non-negative.

Next, when $i=3$, all modular inverse from $1$ to $i-1=2$ is already calculated. This is process is repeated till we reach $N$.

Since we calculated each inverse in $O(1)$, the complexity of this code is $O(N)$.

Conclusion

I saw this code first time on CodeChef forum. I didn't know how it worked back then. I added it to my notebook and have been using it since then. Recently, while searching over the net for resources on Pollard Rho's algorithm, I stumbled on an article from Come On Code On which had the explanation. Thanks, fR0DDY, I have been looking for the proof.

Reference

  1. forthright48 - Modular Multiplicative Inverse
  2. forthright48 - Repeated Squaring Method for Modular Exponentiation
  3. Come On Code On - Modular Multiplicative Inverse

Saturday, September 26, 2015

Euler Phi Extension and Divisor Sum Theorem

Previously we learned about Euler Phi Function. Today we are going to look at two theorems related to Euler Phi that frequently appears in CPPS. I am not sure whether these theorems have any official names, so I just made them up. These allow easy references so I will be using these names from now on.


Euler Phi Extension Theorem

Theorem: Given a number $N$, let $d$ be a divisor of $N$. Then the number of pairs $\{a,N\}$, where $1 \leq a \leq N$ and $gcd(a,N) = d$, is $\phi(\frac{N}{d})$.

Proof

We will prove the theorem using Euler Phi Function and Arithmetic notion.

We need to find the numbe of pairs $\{a,N\}$ such that $gcd(a,N) = d$, where $1 \leq a \leq N$. 

Both $a$ and $N$ are divisible by $d$ and $d$ is the GCD. So, if we divide both $a$ and $N$ by $d$, then they will no longer have any common divisor.

$gcd(\frac{a}{d},\frac{N}{d}) = 1$,  where $1 \leq a \leq N$.

We know that the possible values of $a$ lie in range $1 \leq a \leq N$. What about the possible values of $\frac{a}{d}$? $\frac{a}{d}$ must lie between $1  \leq \frac{a}{d} \leq \frac{N}{d}$ otherwise $a$ will cross its limit.

Therefore, $gcd(a,N) = d$, where $1 \leq a \leq N$ is same as, $gcd(\frac{a}{d},\frac{N}{d}) = 1$, where $1  \leq \frac{a}{d} \leq \frac{N}{d}$.

So all we need to do is find the value of  $gcd(\frac{a}{d},\frac{N}{d}) = 1$, where $1  \leq \frac{a}{d} \leq \frac{N}{d}$.

Let $N' = \frac{N}{d}$ and $a' = \frac{a}{d}$. How many pairs of $\{a',N'\}$ are there such that $gcd(a',N') = 1$ and $1 \leq a' \leq N'$? Isn't this what Euler Phi Function finds? The answer is $\phi(N') = \phi(\frac{N}{d})$.

Euler Phi Divisor Sum Theorem

Theorem: For a given integer $N$, the sum of Euler Phi of each of the divisors of $N$ equals to $N$, i.e, $$N = \sum_{d|N}\phi(d)$$

Proof

The proof is simple. I have broken down the proof in the following chunks for the ease of understanding.

Forming Array $A$

Imagine, we the following fractions in a list: $$\frac{1}{N}, \frac{2}{N}, \frac{3}{N}...\frac{N}{N}$$

Not very hard to imagine right? Let us convert this into an array of pairs. So now, we have the following array $A$: 

$$A = [ \{1,N\},\{2,N\},\{3,N\}...\{N,N\} ]$$

So we have an array of form $\{a,N\}$, where $a$ is between $1$ and $N$. There are exactly $N$ elements in the array.

Finding GCD of Pairs

Next, we find the GCD of each pair, $g$. What are the possible values of $g$? Since $g$ must divide both $a$ and $N$, $g$ must be a divisor of $N$. Therefore, we can conclude that, GCD of pair $\{a,N\}$ will be one of the divisors of $N$.

Let the divisors of $N$ be the following: $d_1, d_2, d_3...d_r$. So these are the only possible GCD.

Forming Parititions

Next, we form partitions $P_i$. Let us put all pairs which have $gcd(a,N) = d_i$ to partition $P_i$. Therefore, we will have $R$ partitions, where $R$ is the number of divisor of $N$. Note that each pair will belong to one partition only since a pair has a unique GCD. Therefore, $$N = \sum_{i=1}^{R}P_i$$

Size of Each Paritition

How many elements does partition $P_i$ contain? $P_i$ has all the pairs $\{a,N\}$ such that $gcd(a,N) = d_i$, $1 \leq a \leq N$. Using Euler Phi Extension Theorem from above, this value is $\phi(\frac{N}{d_i})$.

Wrapping it Up

We are almost done with the proof. Since $N = \sum_{i=1}^{R}P_i$, we can now write: $$N = \sum_{i=1}^{R}P_i = \sum_{i=1}^{R}\phi(\frac{N}{d_i})$$

But $d_i$ is just a divisor of $N$. So we can simplify and write:

$$N = \sum_{d|N}\phi(\frac{N}{d}) = \sum_{d|N}\phi(d)$$

Conclusion

These theorems may look so simple that you might think they are useless. Especially Euler Phi Divisor Sum Theorem, $N = \sum_{d|N} \phi(d)$. How is this useful at all? Hopefully, we will see one of its application on next post.

Reference

Wednesday, September 23, 2015

Modular Multiplicative Inverse

Problem

Given value of $A$ and $M$, find the value of $X$ such that $AX \equiv 1\:\text{(mod M)}$.

For example, if $A = 2$ and $M = 3$, then $X = 2$, since $2\times2 = 4 \equiv 1\:\text{(mod 3)}$.

We can rewrite the above equation to this:

$AX \equiv 1\:\text{(mod M)}$
$X \equiv \frac{1}{A}\:\text{(mod M)}$
$X \equiv A^{-1}\:\text{(mod M)}$

Hence, the value $X$ is known as Modular Multiplicative Inverse of $A$ with respect to $M$.

How to Find Modular Inverse?

First we have to determine whether Modular Inverse even exists for given $A$ and $M$ before we jump to finding the solution. Modular Inverse doesn't exist for every pair of given value.

Existence of Modular Inverse

Modular Inverse of $A$ with respect to $M$, that is, $X = A^{-1} \text{(mod M)}$ exists, if and only if $A$ and $M$ are coprime.

Why is that?

$AX \equiv 1 \:\text{(mod M)}$
$AX - 1 \equiv 0 \:\text{(mod M)}$

Therefore, $M$ divides $AX-1$. Since $M$ divides $AX-1$, then a divisor of $M$ will also divide$AX-1$. Now suppose, $A$ and $M$ are not coprime. Let $D$ be a number greater than $1$ which divides both $A$ and $M$. So, $D$ will divide $AX - 1$. Since $D$ already divides $A$, $D$ must divide $1$. But this is not possible. Therefore, the equation is unsolvable when $A$ and $M$ are not coprime.

From here on, we will assume that $A$ and $M$ are coprime unless state otherwise.

Using Fermat's Little Theorem

Recall Fermat's Little Theorem from a previous post, "Euler's Theorem and Fermat's Little Theorem". It stated that, if $A$ and $M$ are coprime and $M$ is a prime, then, $A^{M-1} \equiv 1 \text{(mod M)}$. We can use this equation to find the modular inverse.

$A^{M-1} \equiv 1 \text{(mod M)}$ (Divide both side by $A$)
$A^{M-2} \equiv \frac{1}{A}\text{(mod M)}$
$A^{M-2} \equiv A^-1\text{(mod M)}$

Therefore, when $M$ is prime, we can find modular inverse by calculating the value of $A^{M-2}$. How do we calculate this? Using Modular Exponentiation.

This is the easiest method, but it doesn't work for non-prime $M$. But no worries since we have other ways to find the inverse.

Using Euler's Theorem

It is possible to use Euler's Theorem to find the modular inverse. We know that:

$A^{\phi(M)} \equiv 1 \text{(mod M)}$
$\therefore A^{\phi(M)-1} \equiv A^{-1} \text{(mod M)}$

This process works for any $M$ as long as it's coprime to $A$, but it is rarely used since we have to calculate Euler Phi value of $M$ which requires more processing. There is an easier way.

Using Extended Euclidean Algorithm

We are trying to solve the congruence, $AX \equiv 1 \text{(mod M)}$. We can convert this to an equation.

$AX \equiv 1 \text{(mod M)}$
$AX + MY = 1$

Here, both $X$ and $Y$ are unknown. This is a linear equation and we want to find integer solution for it. Which means, this is a Linear Diophantine Equation.

Linear Diophantine Equation can be solved using Extended Euclidean Algorithm. Just pass $\text{ext_gcd()}$ the value of $A$ and $M$ and it will provide you with values of $X$ and $Y$. We don't need $Y$ so we can discard it. Then we simply take the mod value of $X$ as the inverse value of $A$.

Code

$A$ and $M$ need to be coprime. Otherwise, no solution exists. The following codes do not check if $A$ and $M$ are coprime. The checking is left of the readers to implement.

When $M$ is Prime

We will use Fermat's Little Theorem here. Just call the $bigmod()$ function from where you need the value. 
int x = bigmod( a, m - 2, m ); ///(ax)%m = 1
Here $x$ is the modular inverse of $a$ which is passed to $bigmod()$ function.

When $M$ is not Prime

For this, we have to use a new function. 
int modInv ( int a, int m ) {
    int x, y;
    ext_gcd( a, m, &x, &y );

    ///Process x so that it is between 0 and m-1
    x %= m;
    if ( x < 0 ) x += m;
    return x;
}
I wrote this function since after using $\text{ext_gcd()}$ we need to process $x$ so that it's value is between $0$ and $M-1$. Instead of doing that manually, I decided to write a function.

So, if we want to find the modular inverse of $A$ with respect to $M$, then the result will be $X = modInv ( A, M )$.

Complexity

Repeated Squaring method has a complexity of $O(logP)$, so the first code has complexity of $O(logM)$, whereas Extended Euclidean has complexity of $O(log_{10}A+log_{10}B)$ so second code has complexity $O(log_{10}A + log_{10}M)$.

Why Do We Need Modular Inverse?

We need Modular Inverse to handle division during Modular Arithmetic. Suppose we are trying to find the value of the following equations:

$\frac{4}{2} \:\%\: 3$ - This is simple. We just simplify the equation and apply normal modular operation. That is, it becomes $\frac{4}{2} \:\%\: 3 = 2 \:\%\: 3 = 2$.

Then what happens when we try to do same with $\frac{12}{9}\:\%\:5$? First we simply. $\frac{12}{9}\:\%\:5 = \frac{4}{3}\:\%\:5$. Now we are facing an irreducible fraction. Should we simply perform the modular operation with numerator and denominator? That doesn't help since both of them are smaller than $5$.

This is where Modular Inverse comes to the rescue. Let us solve the equation $X \equiv 3^{-1}\:\text{(mod 5)}$. How do we find the value of $X$? We will see that on the later part of the post. For now, just assume that we know the value of $X$.

Now, we can rewrite the above equation in the following manner:

$\frac{12}{9}\:\%\:5$
$\frac{4}{3}\:\%\:5$
$(4 \times 3^{-1})\:\%\:5$
$( (4\:\%\:5) \times (3^{-1}\:\%\:5) ) \:\%\:5$
$\therefore 4X \:\%\:5$

So, now we can easily find the value of $\frac{A}{B} \:\%\: M$ by simply calculating the value of $(A \times B^{-1}) \:\%\: M$.

Conclusion

Modular Inverse is a small topic but look at the amount of background knowledge it requires to understand it! Euler's Theorem, Euler Phi, Modular Exponentiation, Linear Diophantine Equation, Extended Euclidian Algorithm and other small bits of information. We covered them all before, so we can proceed without any hitch.

Hopefully, you understood how Modular Inverse works. If not, make sure to revise the articles in the "Reference" section below.

Reference

  1. Wiki - Modular Multiplicative Inverse
  2. forthright48 - Euler's Theorem and Fermat's Little Theorem
  3. forthright48 - Modular Exponentiation
  4. forthright48 - Euler Phi
  5. forthright48 - Linear Diophantine Equation
  6. forthright48 - Extended Euclidean Algorithm

Monday, September 21, 2015

Repeated Squaring Method for Modular Exponentiation

Previously on Modular Exponentiation we learned about Divide and Conquer approach to finding the value of $B^P \:\%\: M $. In that article, which is recursive. I also mentioned about an iterative algorithm that finds the same value in same complexity, only faster due to the absence of recursion overhead. We will be looking into that faster algorithm on this post today.

Make sure you know about Bit Manipulation before proceeding.


Problem

Given three positive integers $B, P$ and $M$, find the value of $B^P \:\%\: M$.

For example, $B=2$, $P=5$ and $M=7$, then $B^P \:\%\: M = 2^5 \: \% \: 7 = 32 \: \% \: 7 = 4$.

Repeated Squaring Method

Repeated Squaring Method (RSM) takes the advantage of the fact that $A^x \times A^y = A^{x+y}$.

Now, we know that any number can be written as the sum of powers of $2$. Just convert the number to Binary Number System. Now for each position $i$ for which binary number has $1$ in it, add $2^i$ to the sum.

For example, $(19)_{10} = (10011)_2 = ( 2^4 + 2^1 + 2^0)_{10} = (16 + 2 + 1 )_{10}$

Therefore, in equation $B^P$, we can decompose $P$ to the sum of powers of $2$.

So let's say, $P = 19$, then $B^{19} = B^{2^4 + 2^1 + 2^0} = B^{16+2+1} = B^{16} \times B^2 \times B^1$.

This is the main concept for repeated squaring method. We decompose the value $P$ to binary, and then for each position $i$ (we start from 0 and loop till the highest position at binary form of $P$) for which binary of $P$ has $1$ in $i_{th}$ position, we multiply $B^{2^i}$ to result.

Code

Here is the code for RSM. The Explanation is below the code.
int bigmod ( int b, int p, int m ) {
    int res = 1 % m, x = b % m;
    while ( p ) {
        if ( p & 1 ) res = ( res * x ) % m;
        x = ( x * x ) % m;
        p >>= 1;
    }
    return res;
}

Explanation

At line $1$, we have the parameters. We simply send the value of $B$, $P$ and $M$ to this function, and it will return the required result.

At line $2$, we initiate some variables. $res$ is the variable that will hold our result. It contains the value $1$ initially. We will multiply $b^{2^i}$ with result to find $b^p$. $x$ is temporary variable that initially contains the value $b^{2^0} = b^1 = b$.

Now, from line $3$ the loop starts. This loop runs until $p$ becomes $0$. Huh? Why is that? Keep reading.

At line $4$ we first check whether the first bit of $p$ is on or not. If it is on, then it means that we have to multiply $b^{2^i}$ to our result. $x$ contains that value, so we multiply $x$ to $res$.

Now line $5$ and $6$ are crucial to the algorithm. Right now, $x$ contains the value of $b^{2^0}$ and we are just checking the $0_{th}$ position of $p$ at each step. We need to update our variables such that they keep working for positions other than $0$.

First, we update the value of $x$. $x$ contains the value of $b^{2^i}$. On next iteration, we will be working with position $i+1$. So we need to update $x$ to hold $b^{2^{i+1}}$.

$b^{2^{i+1}} = b^{2^i \times 2^1} = b ^ {2^i \times 2} = b^{2^i + 2^i} = b^{2^i} \times b^{2^i} = x \times x$.

Hence, new value of $x$ is $x \times x$. We make this update at line $5$.

Now, at each step we are checking the $0_{th}$ position of $p$. But next we need to check the $1_{st}$ position of $p$ in binary. Instead of checking $1_{st}$ position of $p$, what we can do is shift $p$ $1$ time towards right. Now, checking $0_{th}$ position of $p$ is same as checking $1_{st}$ position. We do this update at line $6$.

These two lines ensures that our algorithm works on each iteration. When $p$ becomes $0$, it means there is no more bit to check, so the loop ends.

Complexity

Since there cannot be more than $log_{2}(P)$ bits in $P$, the loop at line $3$ runs at most $log_{2}(P)$ times. So the complexity is $log_{2}(P)$.

Conclusion

RSM is significantly faster than D&C approach due to lack of recursion overhead. Hence, I always use this method when I have to find Modular Exponentiation.

The code may seem a little confusing, so feel free to ask questions.

When I first got my hands on this code, I had no idea how it worked. I found it in a forum with a title, "Faster Approach to Modular Exponentiation". Since then I have been using this code.

Resources

  1. forthrigth48 - Modular Exponentiation
  2. forthright48 - Bit Manipulation
  3. algorithmist - Repeated Squaring
  4. Wiki - Exponentiation by Squaring

Thursday, September 17, 2015

Euler's Theorem and Fermat's Little Theorem

We will be looking into two theorems at the same time today, Fermat's Little Theorem and Euler's Theorem. Euler's Theorem is just a generalized version of Fermat's Little Theorem, so they are quite similar to each other. We will focus on Euler's Theorem and its proof. Later we will use Euler's Theorem to prove Fermat's Little Theorem.

Euler's Theorem

Theorem - Euler's Theorem states that, if $a$ and $n$ are coprime, then $a^{\phi(n)} \equiv 1 (\text{mod n})$ - Wikipedia
Here $\phi(n)$ is Euler Phi Function. Read more about Phi Function on this post - Euler Totient or Phi Function.

Proof

Let us consider a set $A = \{b_1, b_2, b_3...,b_{\phi(n)} \}\:(\text{mod n})$, where $b_i$ is coprime to $n$ and distinct. Since there are $\phi(n)$ elements which are coprime to $n$, $A$ contains $\phi(n)$ integers.

Now, consider the set $B = \{ ab_1, ab_2, ab_3....ab_{\phi(n)} \} \:(\text{mod n})$. That is, $B$ is simply set $A$ where we multiplied $a$ with each element. Let $a$ be coprime to $n$. 
Lemma - Set $A$ and set $B$ contains the same integers.
We can prove the above lemma in three steps.
  1. $A$ and $B$ has the same number of elements
    Since $B$ is simply every element of $A$ multiplied with $a$, it contains the same number of elements as $A$. This is obvious.
  2. Every integer in $B$ is coprime to $n$
    An integer in $B$ is of form $a \times b_i$. We know that both $b_i$ and $a$ are coprime to $n$, so $ab_i$ is also coprime to $n$.
  3. $B$ contains distinct integers only
    Suppose $B$ does not contain distinct integers, then it would mean that there is such a $b_i$ and $b_j$ such that:

    $ab_i \equiv ab_j\:(\text{mod n})$
    $b_i \equiv b_j\:(\text{mod n})$

    But this is not possible since all elements of $A$ are distinct, that is, $b_i$ is never equal to $b_j$. Hence, $B$ contains distinct elements.
With these three steps, we claim that, since $B$ has the same number of elements as $A$ which are distinct and coprime to $n$, it has same elements as $A$.

Now, we can easily prove Euler's Theorem.

$ab_1 \times ab_2 \times ab_3...\times ab_{\phi(n)} \equiv b_1 \times b_2 \times b_3...\times b_{\phi(n)} \:(\text{mod n})$
$a^{\phi(n)} \times b_1 \times b_2 \times b_3...\times b_{\phi(n)} \equiv b_1 \times b_2 \times b_3...\times b_{\phi(n)}  \:(\text{mod n}) $
$$\therefore a^{\phi(n)} \equiv 1  \:(\text{mod n})$$

Fermat's Little Theorem

Fermat's Little Theorem is just a special case of Euler's Theorem.
Theorem - Fermat's Little Theorem states that, if $a$ and $p$ are coprime and $p$ is a prime, then $a^{p-1} \equiv 1 \:(\text{mod p})$ - Wikipedia
As you can see, Fermat's Little Theorem is just a special case of Euler's Theorem. In Euler's Theorem, we worked with any pair of value for $a$ and $n$ where they are coprime, here $n$ just needs to be prime.

We can use Euler's Theorem to prove Fermat's Little Theorem.

Let $a$ and $p$ be coprime and $p$ be prime, then using Euler's Theorem we can say that:

$a^{\phi(p)} \equiv 1\:(\text{mod p})$  (But we know that for any prime $p$, $\phi(p) = p-1$)
$a^{p-1} \equiv 1\:(\text{mod p})$

Conclusion

Both theorems have various applications. Finding Modular Inverse is a popular application of Euler's Theorem. It can also be used to reduce the cost of modular exponentiation. Fermat's Little Theorem is used in Fermat's Primality Test.

There are more applications but I think it's better to learn them as we go. Hopefully, I will be able to cover separate posts for each of the applications.

Reference

  1. Wiki - Euler's Theorem
  2. forthright48 - Euler Totient or Phi Function
  3. Wiki - Fermat's Little Theorem

Monday, September 7, 2015

Segmented Sieve of Eratosthenes

Problem

Given two integers $A$ and $B$, find number of primes inside the range of $A$ and $B$ inclusive. Here, $1 \leq A \leq B \leq 10^{12}$ and $B - A \leq 10^5$.

For example, $A = 11$ and $B = 19$, then answer is $4$ since there are $4$ primes within that range ($11$,$13$,$17$,$19$).

If limits of $A$ and $B$ were small enough ( $ \leq 10^8$ ), then we could solve this problem using the ordinary sieve. But here limits are huge, so we don't have enough memory or time to run normal sieve. But note that, $B - A \leq 10^5$. So even though we don't have memory/time to run sieve from $1$ to $N$, we have enough memory/time to cover $A$ to $B$.

$A$ to $B$ is a segment, and we are going to modify our algorithm for Sieve of Eratosthenes to cover this segment. Hence, the modified algorithm is called Segmented Sieve of Eratosthenes. 

Make sure you fully understand how "Normal" Sieve of Eratosthenes works.

How Normal Sieve Works

First, let us see the steps for "Normal" sieve. In order to make things simpler, we will be looking into a somewhat unoptimized sieve. You can implement your own optimizations later.

Suppose we want to find all primes between $1$ to $N$.
  1. First we define a new variable $sqrtn = \sqrt{N}$.
  2. We take all primes less than $sqrtn$. 
  3. For each prime $p$, we repeat the following steps.
    1. We start from $j = p \times p$.
    2. If $j \leq N$, we mark sieve at position $j$ to be not prime.
      Else, we break out of the loop.
    3. We increase the value of $j$ by $p$. And go back to step step $2$ of this loop.
  4. All positions in the sieve that are not marked are prime.
This is how the basic sieve works. We will now modify it to work on segments.

How Segmented Sieve Works

We will perform the same steps as normal sieve but just slightly modified.

Generate Primes Less Than $\sqrt{N}$

In the segmented sieve, what is he largest limit possible? $10^{12}$. So let $N = 10^{12}$

First of all, in the normal sieve we worked with primes less than $\sqrt{N}$ only. So, if we had to run sieve from $1$ to $N$, we would have required only primes less than $\sqrt{N} = 10^6$. So in order to run sieve on a segment between $1$ to $N$, we won't require primes greater than $\sqrt{N}$.

So, using normal sieve we will first generate all primes less than $\sqrt{N} = 10^6$.

Run on Segment

Okay, now we can start our "Segmented" Sieve. We want to find primes between $A$ and $B$. 
  1. If $A$ is equal to $1$, then increase $A$ by $1$. That is, make $A = 2$. Since $1$ is not a prime, this does not change our answer.
  2. Define a new variable $sqrtn = \sqrt{B}$.
  3. Declare a new array of size $dif = \text{maximum difference of }(B - A) + 1$. Since it is given in our problem that $B-A \leq 10^5$, $dif = 10^5 + 1$ for this problem.

    Let the array be called $arr$. This array has index from $0$ to $dif-1$. Here $arr[0]$ represents the number $A$, $arr[1]$ represents $A+1$ and so on.
  4. Now, we will be working with all primes less than $sqrtn$. These primes are already generated using the normal sieve.
  5. For each prime $p$, we will repeat the following steps:
    1. We start from $j = p \times p$.
    2. But initial value of $j = p^2$ might be less than $A$. We want to mark positions between $A$ and $B$ only. So we will need to shift $j$ inside the segment.

      So, if $j < A$, then $j = ceil ( \frac{A}{p} ) \times p = \frac{A+p-1}{p} \times p$. This line makes $j$ the smallest multiple of $p$ which is bigger than $A$.
    3. If $j \leq B$, we mark sieve at position $j$ to be not prime.
      Else, we break out of the loop.

      But when marking, remember that our array $arr$ is shifted by $A$ positions. $arr[0]$ indicates position $A$ of normal sieve. So, we will mark position $j - A$ of $arr$ as not prime.
    4. Increase the value of $j$ by $p$. Repeat loop from step $3$.
  6. All positions in $arr$ which has not been marked are prime.
Step $1$ is important. Since we only mark multiples of prime as not prime in the pseducode above, $1$ which has no prime factor never gets marked. So we handle it by increassing value of $A$ by $1$ when $A = 1$.

Code

If we convert the above pseducode into C++, then it becomes something like this:
int arr[SIZE];

///Returns number of primes between segment [a,b]
int segmentedSieve ( int a, int b ) {
    if ( a == 1 ) a++;

    int sqrtn = sqrt ( b );

    memset ( arr, 0, sizeof arr ); ///Make all index of arr 0.

    for ( int i = 0; i < prime.size() && prime[i] <= sqrtn; i++ ) {
        int p = prime[i];
        int j = p * p;

        ///If j is smaller than a, then shift it inside of segment [a,b]
        if ( j < a ) j = ( ( a + p - 1 ) / p ) * p;

        for ( ; j <= b; j += p ) {
            arr[j-a] = 1; ///mark them as not prime
        }
    }

    int res = 0;
    for ( int i = a; i <= b; i++ ) {
        ///If it is not marked, then it is a prime
        if ( arr[i-a] == 0 ) res++;
    }
    return res;
}
In line $1$ we declare an array of $SIZE$. This array needs to be as large as maximum difference between $B-A+1$. Next in line $4$ we declare a function that finds number of primes between $a$ and $b$. Rest of the code is same as the psedocode above.

It is possible to optimize it further ( both for speed and memory ) but in expense of clarity. I am sure if readers understand the core concept behind this algorithm, then they will have no problem tweaking the code to suit their needs.

Conclusion

I first learned about Segmented Sieve from blog of +Zobayer Hasan. You can have a look at that post here. I wasn't really good at bit manipulation back then, so it looked really scary. Later I realized it's not as hard as it looks. Hopefully, you guys feel the same.

Leave a comment if you face any difficulty in understanding the post.

Reference

  1. forthright48 - Sieve of Eratosthenes
  2. zobayer - Segmented Sieve

Related Problems


Friday, September 4, 2015

Euler Totient or Phi Function

I have been meaning to write a post on Euler Phi for a while now, but I have been struggling with its proof. I heard it required Chinese Remainder Theorem, so I have been pushing this until I covered CRT. But recently, I found that CRT is not required and it can be proved much more easily. In fact, the proof is so simple and elegant that after reading it I went ahead and played MineCraft for 5 hours to celebrate.


Problem

Given an integer $N$, how many numbers less than or equal $N$ are there such that they are coprime to $N$? A number $X$ is coprime to $N$ if $gcd(X,N)=1$.

For example, if $N = 10$, then there are $4$ numbers, namely $\{1,3,7,9\}$, which are coprime to $10$.

This problem can be solved using Euler Phi Function, $\phi()$. Here is the definition from Wiki:
In number theory, Euler's totient function (or Euler's phi function), denoted as φ(n) or ϕ(n), is an arithmetic function that counts the positive integers less than or equal to n that are relatively prime to n. - Wiki
That's exactly what we need to find in order to solve the problem above. So, how does Euler Phi work?

Euler Phi Function

Before we go into its proof, let us first see the end result. Here is the formula using which we can find the value of the $\phi()$ function. If we are finding Euler Phi of $N = p_1^{a_1}p_2^{a_2}...p_k^{a_k}$, then:
$$\phi(n) = n \times \frac{p_1-1}{p_1}\times\frac{p_2-1}{p_2}...\times\frac{p_k-1}{p_k}$$
If you want you can skip the proof and just use the formula above to solve problems. That's what I have been doing all these years. But I highly recommend that you read and try to understand the proof. It's simple and I am sure someday the proof will help you out in an unexpected way.

Proof of Euler Phi Function

Even though the proof is simple, it has many steps. We will go step by step, and slowly you will find that the proof is unfolding in front of your eyes.

Base Case - $\phi(1)$

First, the base case. Phi function counts the number of positive numbers less than $N$ that are coprime to it. The keyword here is positive. Since the smallest positive number is $1$, we will start with this.

$\phi(1) = 1$, since $1$ itself is the only number which is coprime to it.

When $n$ is a Prime - $\phi(p)$

Next, we will consider the case when $n = p$. Here $p$ is any prime number. When $n$ is prime, it is coprime to all numbers less than $n$. Therefore, $\phi(n) = \phi(p) = p - 1$.

When $n$ is Power of Prime - $\phi(p^a)$

Next, we will consider $n$ where $n$ is a power of a single prime. In this case, how many numbers less than $n$ are coprime to it? Instead of counting that, we will count the inverse. How many numbers are there which are not coprime.

Since, $n = p^a$, we can be sure that $gcd(p,n) \neq 1$. Since both $n$ and $p$ are divisible by $p$. Therefore, the following numbers which are divisible by $p$ are not coprime to $n$,  $\{p, 2p, 3p$ $.... p^2, (p+1)p, (p+2)p$ $...(p^2)p, (p^2+1)p...(p^{a-1})p \}$. There are exactly $\frac{p^a}{p} = p^{a-1}$ numbers which are divisible by $p$. So, there are $n - p^{a-1}$ numbers which are coprime to $n$.

Hence, $\phi(n) = \phi(p^a) $ $ = n - \frac{n}{p} = p^a - \frac{p^a}{p} $ $= p^a ( 1 - \frac{1}{p} ) = p^a \times ( \frac{p - 1}{p} )$

It's starting to look like the equation above, right?

Assuming $\phi()$ is Multiplicative - $\phi( m \times n )$

This step is the most important step in the proof. This step claims that Euler Phi function is a multiplicative function. What does this mean? It means, if $m$ and $n$ are coprime, then $\phi( m \times n ) = \phi(m) \times \phi(n) $. Functions that satisfy this condition are called Multiplicative Functions.

So how do we prove that Euler Phi is multiplicative and how does Euler Phi being multiplicative helps us?

We will prove multiplicity of Euler Phi Function in the next section. In this section, we will assume it is multiplicative and see how it helps us calculating Euler Phi.

Let the prime factorization of $n$ be $p_1^{a_1}p_2^{a_2}...p_k^{a_k}$. Now, obviously $p_i$ nad $p_j$ are coprime to each other. Since $\phi$ funtion is multiplicative, we can simply rewrite the function as:

$\phi(n) = \phi(p_1^{a_1}p_2^{a_2}...p_k^{a_k})$
$\phi(n) = \phi(p_1^{a_1}) \times  \phi(p_2^{a_2}) ... \times  \phi(p_k^{a_k})$.

We can already calculate $\phi(p^a) = p^a \times \frac{p-1}{p}$. So our equationg becomes:

$\phi(n) = \phi(p_1^{a_1}) \times  \phi(p_2^{a_2}) ... \times  \phi(p_k^{a_k})$
$\phi(n) = p_1^{a_1} \times \frac{p_1 - 1}{p_1} \times  p_2^{a_2} \times \frac{p_2 - 1}{p_2}...\times  p_k^{a_k} \times \frac{p_k - 1}{p_k}$
$\phi(n) = ( p_1^{a_1} \times p_2^{a_2}... \times p_k^{a_k} ) \times \frac{p_1 - 1}{p_1} \times \frac{p_2 - 1}{p_2}... \times \frac{p_k - 1}{p_k}$
$$\therefore \phi(n) = n \times \frac{p_1 - 1}{p_1} \times \frac{p_2 - 1}{p_2}... \times \frac{p_k - 1}{p_k}$$
This is what we have been trying to prove. This equation was derived by assuming that Euler Phi Function is multiplicative. So all we need to do now is prove Euler Phi Function is multiplicative and we are done.

Proof for Multiplicity of Euler Phi Function

We are trying to prove the following theorem:
Theorem $1$: If $m$ and $n$ are coprime, then $\phi(m \times n ) = \phi ( m ) \times \phi(n)$
But in order to prove Theorem $1$, we will need to prove few other theorems first.

Theorem Related to Arithmetic Progression

Theorem $2$: In an arithmetic progression with difference of $m$, if we take $n$ terms and find their modulo by $n$, and if $n$ and $m$ are coprimes, then we will get the numbers from $0$ to $n-1$ in some order.
Umm, looks like Theorem $2$ is packed with too much information. Let me break it down.

Suppose you have an arithmetic progression (AP). Now, every arithmetic progression has two things. A starting element and a common difference. That is, arithmetic progressions are of the form $a + kb$ where $a$ is the starting element, $b$ is the common difference and $k$ is any number.

So take any arithmetic progression that has a common difference of $m$. Then take $n$ consecutive terms of that progression. So which terms will they be?

$\{ a + 0m, a + 1m, a + 2m, a + 3m ... a + (n-1)m\}$ There are exactly $n$ terms in this list.

Next, find their modulus by $n$. That is, find the remainder of each term after dividing by $n$.

Now, Theorem $2$ claims that, if $m$ and $n$ are coprime, then the list will contain a permutation of $0$ to $n-1$.

For example, let us try with $a = 1$, $m = 7$ and $n = 3$. So the $3$ terms in the list will be $(1, 8, 15)$. Now, if find modulus of each element, we get $(1,2,0)$.

I hope that now it's clear what the Theorem claims. Now we will look into the proof.

Proof of Theorem 2

What we need to prove is that the list after modulo operation has a permutation of numbers from $0$ to $n-1$. That means, all the numbers from $0$ to $n-1$ occurs in the list exactly once. There are three steps to this proof. Also, remember that $m$ and $n$ are coprime.
  1. There are exactly $n$ elements in the list.
    Well, since we took $n$ terms from the AP, this is obvious.
  2. Each element of the list has value between $0$ to $n-1$
    We performed modulo operations on each element by $n$. So this is also obvious.
  3. No remainder has the same value as another.
    Since there are $n$ values, and each value is between $0$ to $n-1$, if we can prove that each element is unique in the list, then our work is done.

    Suppose there are two numbers which have the same remainder. That means $a + pm$ has same remainder as $a + qm$, where $p$ and $q$ are two integer numbers such that $0 \leq p < q \leq n - 1$.

    Therefore, $( a + qm ) - ( a + pm ) \equiv 0\: ( mod\: n )$
    $(a+qm -a - pm) \equiv 0\: ( mod\: n )$
    $m ( q - p ) \equiv 0\: ( mod\: n )$

    That means, $n$ divides $m(q-p)$. But this is impossible. $n$ and $m$ are coprime and $q-p$ is smaller than $n$. So it is not possible for two numbers to have the same remainder.

Theorem Related to Remainder

Theorem $3$: If a number $x$ is coprime to $y$, then $(x \:\%\: y)$ will also be coprime to $y$.
The proof for this theorem is simple. Suppose $x$ and $y$ are coprime. Now, we can rewrite $x$ as $x = ky + r$, where $k$ is the quotient and $r$ is the remainder.

Theorem $3$ claims that $y$ and $r$ are coprime. What happens if this claim is false? Suppose they are not coprime. That means there is a number $d > 1$ which divides both $y$ and $r$. But then, $d$ also divides $ky + r = a$. So $d > 1$ divides $r, y$ and $x$, which is impossible cause $y$ and $x$ are coprime. There is no number greater than $1$ that can divide both of them. Hence, there is a contradiction. Hence, the theorem is proved.

Proof for Multiplicity of Euler Phi Function Continued

We are now ready to tackle proving Theorem $1$.

Suppose you have two numbers $m$ and $n$, and they are coprime. We want to show that $\phi(m\times n) = \phi (m ) \times \phi (n )$. 

What done $\phi (m\times n)$ gives us? It gives us the count of numbers which are coprime to $mn$. If a number $x$ is coprime to $mn$, then it is also coprime to $m$ and $n$ separately. So basically, we need to count the number of positive numbers less than or equal to $mn$ which are coprime to both $m$ and $n$.

Now, let us build a table of with $n$ rows and $m$ columns. Therefore, the table will look like the following:
123...m
1+m2+m3+m...2m
1+2m2+2m3+2m...3m
...............
1 + (n-1)m2 + (n-1)m3 + (n-1)m...mn
Now, notice that each column is an arithmetic progression with $n$ terms and has common difference of $m$. Also, $m$ and $n$ are coprime. This is exactly the same situation as Theorem $2$.

Now, how many numbers in each column are coprime to $n$? In order to figure this result out, we first need to consider what happens if we modulo all table values with $n$. Using theorem $2$, we know that each column will then contain a permutation of numbers from $0$ to $n-1$. Using theorem $3$, we know what if the remainder of a number is coprime to $n$ then the number itself will also be coprime. So, how many numbers between $0$ to $n-1$ is coprime to $n$? We can consider $0$ to be same as $n$ ( cause this is modular arithmetic), so it boils down to, how many numbers between $1$ to $n$ is coprime to $n$? Euler Phi Function calculates this values.

So, there are exactly $\phi(n)$ numbers which are coprime to $n$ in each column.

We need to find numbers that coprime to both $n$ and $m$. So, we cannot take $\phi(n)$ elements from every column, cause those elements may not be coprime to $m$. How do we decide which columns we should be taking?

Notice that, if we find the modulus of elements of the table by $m$, then each row has remainder between $0$ to $m-1$ occuring exactly once. If we consider $0$ to be $m$, then each row has values between $1$ to $m$.  That is the table becomes something like this:
123...m
123...m
123...m
...............
123...m
So, how many columns are there which are coprime to $m$? There are $\phi(m)$ columns which are coprime to $m$.

Now we just need to combine the two results from above. There are exactly $\phi(m)$ columns which are coprime to $m$ and in each column there are $\phi(n)$ values which are coprime to $n$. Therefore, there are $\phi(m) \times \phi(n)$ elements which are coprime to both $m$ and $n$.
$$\therefore \phi(m) \times \phi(n) = \phi(m \times n)$$

Code

Since we have to factorize $n$ in order to calculate $\phi(n)$, we can modify our $factorize()$ function from post "Prime Factorization of Integer Number" to handle Euler Phi.
int eulerPhi ( int n ) {
    int res = n;
    int sqrtn = sqrt ( n );
    for ( int i = 0; i < prime.size() && prime[i] <= sqrtn; i++ ) {
        if ( n % prime[i] == 0 ) {
            while ( n % prime[i] == 0 ) {
                n /= prime[i];
            }
            sqrtn = sqrt ( n );
            res /= prime[i];
            res *= prime[i] - 1;
        }
    }
    if ( n != 1 ) {
        res /= n;
        res *= n - 1;
    }
    return res;
}
I highlighted the lines that are different from $factorize()$ function. Notice that in line $10$ divided $res$ before multiplying in line $11$. This is an optimization that lowers the risk of overflowing.

Conclusion

That was a long post with lots of theorems and equations, but hopefully, they were easy to understand. Even though Theorem $2$ and $3$ were used as lemmas to prove Theorem $1$, they both are important by themselves. 

Leave a comment if you face any difficulty in understanding the post.

Reference

  1. Wiki - Euler Totient Function