tag:blogger.com,1999:blog-8113711071749533092017-08-15T04:53:27.492+06:00forthright48Learning Never EndsMohammad Samiul Islamnoreply@blogger.comBlogger35125tag:blogger.com,1999:blog-811371107174953309.post-91332151576128007172015-09-29T22:20:00.000+06:002017-01-07T09:15:53.849+06:00Modular Inverse from 1 to NWe already learned how to find Modular Inverse for a particular number in a previous post, "<a href="http://forthright48.blogspot.com/2015/09/modular-multiplicative-inverse.html">Modular Multiplicative Inverse</a>". Today we will look into finding Modular Inverse in a bulk.<br /><hr /><br /><h2>Problem</h2>Given $N$ and $M$ ( $N < M$ and $M$ is prime ), find modular inverse of all numbers between $1$ to $N$ with respect to $M$.<br /><br />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.<br /><br />We will look into two methods. Later one is better than the first one.<br /><br /><h2>$O(NlogM)$ Solution</h2><div>Using Fermat's little theorem, we can easily find Modular Inverse for a particular number.</div><div><br /></div><div>$A^{-1} \:\%\: M = bigmod(A,M-2,M)$, where $bigmod()$ is a function from the post "<a href="http://forthright48.blogspot.com/2015/09/repeated-squaring-method-for-modular.html">Repeated Squaring Method for Modular Exponentiation</a>". 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.</div><pre class="brush:cpp">int inv[SIZE]; ///inv[x] contains value of (x^-1 % m)<br />for ( int i = 1; i <= n; i++ ) {<br /> inv[i] = bigmod ( i, m - 2, m );<br />}</pre><div>But it's possible to do better. </div><div><br /></div><h2>$O(N)$ Solution</h2><div>This solution is derived using some clever manipulation of Modular Arithmetic.</div><div><br /></div><div>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. </div><div><br /></div><div>$M = Q \times a + r$, <span style="background-color: #ffd966;">(where $Q$ is the quotient and $r$ is the remainder) </span></div><div>$M = \lfloor \frac{M}{a} \rfloor \times a + (M \:\%\: a )$</div><div><br /></div><div>Now take modulo $M$ on both sides.</div><div><br /></div><div>$0 \equiv \lfloor \frac{M}{a} \rfloor \times a + (M \:\%\: a ) \:\:\:\text{(mod M )}$</div><div>$ (M \:\%\: a ) \equiv -\lfloor \frac{M}{a} \rfloor \times a \:\:\:\text{(mod M )}$</div><div><br /></div><div>Now divide both side by $a \times ( M \:\%\: a )$.</div><div><br /></div><div>$$\frac{M \:\%\: a}{a \times ( M \:\%\: a )} \equiv \frac{- \lfloor \frac{M}{a} \rfloor \times a } { a \times ( M \:\%\: a ) } \:\:\:\text{(mod M)}$$</div><div>$$\therefore a^{-1} \equiv - \lfloor \frac{M}{a} \rfloor \times ( M \:\%\: a )^{-1} \:\:\:\text{(mod M)}$$</div><div><br /></div><div>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. </div><div><br /></div><div>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$.</div><div><br /></div><div>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)$.</div><div><br /></div><h2>Code</h2><div>We can now formulate our code.<br /><pre class="brush:cpp">int inv[SIZE];<br />inv[1] = 1;<br />for ( int i = 2; i <= n; i++ ) {<br /> inv[i] = (-(m/i) * inv[m%i] ) % m;<br /> inv[i] = inv[i] + m;<br />}<br /></pre>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$.<br /><br />At line $5$, we make sure the modular inverse is non-negative.<br /><br />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$.<br /><br />Since we calculated each inverse in $O(1)$, the complexity of this code is $O(N)$.<br /><br /><h2>Conclusion</h2>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 <a href="https://comeoncodeon.wordpress.com/2011/10/09/modular-multiplicative-inverse/">Come On Code On</a> which had the explanation. Thanks, <a href="https://comeoncodeon.wordpress.com/">fR0DDY</a>, I have been looking for the proof.<br /><br /><h2>Reference</h2><ol><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/09/modular-multiplicative-inverse.html">Modular Multiplicative Inverse</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/09/repeated-squaring-method-for-modular.html">Repeated Squaring Method for Modular Exponentiation</a></li><li>Come On Code On - <a href="https://comeoncodeon.wordpress.com/2011/10/09/modular-multiplicative-inverse/">Modular Multiplicative Inverse</a></li></ol></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com5tag:blogger.com,1999:blog-811371107174953309.post-46876640281923208592015-09-26T14:57:00.000+06:002015-09-26T15:00:23.154+06:00Euler Phi Extension and Divisor Sum TheoremPreviously we learned about <a href="http://forthright48.blogspot.com/2015/09/euler-totient-or-phi-function.html">Euler Phi Function</a>. 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.<br /><hr /><br /><h2>Euler Phi Extension Theorem</h2><blockquote class="tr_bq theorem"><b>Theorem</b>: 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})$.</blockquote><h3>Proof</h3><div>We will prove the theorem using Euler Phi Function and Arithmetic notion.</div><div><br />We need to find the numbe of pairs $\{a,N\}$ such that $gcd(a,N) = d$, where $1 \leq a \leq N$. </div><div><br /></div><div>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.</div><div><br /></div><div>$gcd(\frac{a}{d},\frac{N}{d}) = 1$, where $1 \leq a \leq N$.</div><div><br /></div><div>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.</div><div><br /></div><div>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}$.</div><div><br /></div><div>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}$.</div><div><br />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})$.</div><div><br /></div><div><h2>Euler Phi Divisor Sum Theorem</h2><blockquote class="tr_bq theorem"><b>Theorem</b>: 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)$$</blockquote><br /><h3>Proof</h3></div><div>The proof is simple. I have broken down the proof in the following chunks for the ease of understanding.</div><div><br /></div><h3>Forming Array $A$</h3><div>Imagine, we the following fractions in a list: $$\frac{1}{N}, \frac{2}{N}, \frac{3}{N}...\frac{N}{N}$$</div><div><br /></div><div>Not very hard to imagine right? Let us convert this into an array of pairs. So now, we have the following array $A$: </div><div><br /></div><div>$$A = [ \{1,N\},\{2,N\},\{3,N\}...\{N,N\} ]$$</div><div><br /></div><div>So we have an array of form $\{a,N\}$, where $a$ is between $1$ and $N$. There are exactly $N$ elements in the array.</div><div><br /></div><h3>Finding GCD of Pairs</h3><div>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$.</div><div><br /></div><div>Let the divisors of $N$ be the following: $d_1, d_2, d_3...d_r$. So these are the only possible GCD.</div><div><br /></div><h3>Forming Parititions</h3><div>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$$</div><div><br /></div><h3>Size of Each Paritition</h3><div>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})$.</div><div><br /></div><h3>Wrapping it Up</h3><div>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})$$</div><div><br /></div><div>But $d_i$ is just a divisor of $N$. So we can simplify and write:</div><div><br /></div><div>$$N = \sum_{d|N}\phi(\frac{N}{d}) = \sum_{d|N}\phi(d)$$</div><div><br /></div><div><h2>Conclusion</h2></div><div>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.<br /><br /><h2>Reference</h2><div><ol><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/09/euler-totient-or-phi-function.html">Euler Totient or Phi Function</a></li><li>Wiki - <a href="https://en.wikipedia.org/wiki/Euler%27s_totient_function#Divisor_sum">Divisor Sum</a></li></ol></div></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com1tag:blogger.com,1999:blog-811371107174953309.post-53243161863351745722015-09-23T17:39:00.000+06:002015-09-26T10:48:06.233+06:00Modular Multiplicative Inverse<h2>Problem</h2>Given value of $A$ and $M$, find the value of $X$ such that $AX \equiv 1\:\text{(mod M)}$.<br /><br />For example, if $A = 2$ and $M = 3$, then $X = 2$, since $2\times2 = 4 \equiv 1\:\text{(mod 3)}$.<br /><br />We can rewrite the above equation to this:<br /><br />$AX \equiv 1\:\text{(mod M)}$<br />$X \equiv \frac{1}{A}\:\text{(mod M)}$<br />$X \equiv A^{-1}\:\text{(mod M)}$<br /><br />Hence, the value $X$ is known as Modular Multiplicative Inverse of $A$ with respect to $M$.<br /><br /><div><h2>How to Find Modular Inverse?</h2></div><div>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.<br /><br /></div><h3>Existence of Modular Inverse</h3><div><span style="background-color: yellow;">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.</span></div><div><br /></div><div>Why is that?</div><div><br /></div><div>$AX \equiv 1 \:\text{(mod M)}$</div><div>$AX - 1 \equiv 0 \:\text{(mod M)}$</div><div><br /></div><div>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.<br /><br />From here on, we will assume that $A$ and $M$ are coprime unless state otherwise.<br /><br /><h3>Using Fermat's Little Theorem</h3></div><div>Recall Fermat's Little Theorem from a previous post, "<a href="http://forthright48.blogspot.com/2015/09/eulers-theorem-and-fermats-little.html">Euler's Theorem and Fermat's Little Theorem</a>". 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.<br /><br />$A^{M-1} \equiv 1 \text{(mod M)}$ <span style="background-color: #ffd966;">(Divide both side by $A$)</span><br />$A^{M-2} \equiv \frac{1}{A}\text{(mod M)}$<br />$A^{M-2} \equiv A^-1\text{(mod M)}$<br /><br />Therefore, when $M$ is prime, we can find modular inverse by calculating the value of $A^{M-2}$. How do we calculate this? Using <a href="http://forthright48.blogspot.com/2015/09/repeated-squaring-method-for-modular.html">Modular Exponentiation</a>.<br /><br />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.<br /><br /><h3>Using Euler's Theorem</h3>It is possible to use Euler's Theorem to find the modular inverse. We know that:<br /><br />$A^{\phi(M)} \equiv 1 \text{(mod M)}$<br />$\therefore A^{\phi(M)-1} \equiv A^{-1} \text{(mod M)}$<br /><br />This process works for any $M$ as long as it's coprime to $A$, but it is rarely used since we have to calculate <a href="http://forthright48.blogspot.com/2015/09/euler-totient-or-phi-function.html">Euler Phi</a> value of $M$ which requires more processing. There is an easier way.<br /><br /><h3>Using Extended Euclidean Algorithm</h3>We are trying to solve the congruence, $AX \equiv 1 \text{(mod M)}$. We can convert this to an equation.<br /><br />$AX \equiv 1 \text{(mod M)}$<br />$AX + MY = 1$<br /><br />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 <a href="http://forthright48.blogspot.com/2015/07/linear-diophantine-equation.html">Linear Diophantine Equation</a>.<br /><br />Linear Diophantine Equation can be solved using <a href="http://forthright48.blogspot.com/2015/07/extended-euclidean-algorithm.html">Extended Euclidean Algorithm</a>. 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$.<br /><br /><h2>Code</h2></div><div>$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.<br /><br /></div><h3>When $M$ is Prime</h3><div>We will use Fermat's Little Theorem here. Just call the $bigmod()$ function from where you need the value. </div><div><pre class="brush:cpp">int x = bigmod( a, m - 2, m ); ///(ax)%m = 1</pre></div>Here $x$ is the modular inverse of $a$ which is passed to $bigmod()$ function.<br /><br /><div><h3>When $M$ is not Prime</h3></div><div>For this, we have to use a new function. </div><div><pre class="brush:cpp">int modInv ( int a, int m ) {<br /> int x, y;<br /> ext_gcd( a, m, &x, &y );<br /><br /> ///Process x so that it is between 0 and m-1<br /> x %= m;<br /> if ( x < 0 ) x += m;<br /> return x;<br />}<br /></pre></div><div>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.<br /><br />So, if we want to find the modular inverse of $A$ with respect to $M$, then the result will be $X = modInv ( A, M )$.</div><div><br /></div><h2>Complexity</h2><div>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)$.<br /><br /></div><div><h2>Why Do We Need Modular Inverse?</h2>We need Modular Inverse to handle division during Modular Arithmetic. Suppose we are trying to find the value of the following equations:<br /><br />$\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$.<br /><br />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$.<br /><br />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$.<br /><br />Now, we can rewrite the above equation in the following manner:<br /><br />$\frac{12}{9}\:\%\:5$<br />$\frac{4}{3}\:\%\:5$<br />$(4 \times 3^{-1})\:\%\:5$<br />$( (4\:\%\:5) \times (3^{-1}\:\%\:5) ) \:\%\:5$<br />$\therefore 4X \:\%\:5$<br /><br /><div>So, now we can easily find the value of $\frac{A}{B} \:\%\: M$ by simply calculating the value of $(A \times B^{-1}) \:\%\: M$.<br /><br /><h2>Conclusion</h2>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.<br /><br />Hopefully, you understood how Modular Inverse works. If not, make sure to revise the articles in the "Reference" section below.<br /><br /><h2>Reference</h2></div><ol><li>Wiki - <a href="https://en.wikipedia.org/wiki/Modular_multiplicative_inverse">Modular Multiplicative Inverse</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/09/eulers-theorem-and-fermats-little.html">Euler's Theorem and Fermat's Little Theorem</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/09/repeated-squaring-method-for-modular.html">Modular Exponentiation</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/09/euler-totient-or-phi-function.html">Euler Phi</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/07/linear-diophantine-equation.html">Linear Diophantine Equation</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/07/extended-euclidean-algorithm.html">Extended Euclidean Algorithm</a></li></ol></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-42252168373449683112015-09-21T15:10:00.001+06:002015-09-26T10:48:13.633+06:00Repeated Squaring Method for Modular Exponentiation Previously on <a href="http://forthright48.blogspot.com/2015/08/modular-exponentiation.html">Modular Exponentiation</a> 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.<br /><br />Make sure you know about <a href="http://forthright48.blogspot.com/2015/08/bit-manipulation.html">Bit Manipulation</a> before proceeding.<br /><hr /><br /><h2>Problem</h2><div>Given three positive integers $B, P$ and $M$, find the value of $B^P \:\%\: M$.</div><div><br /></div><div>For example, $B=2$, $P=5$ and $M=7$, then $B^P \:\%\: M = 2^5 \: \% \: 7 = 32 \: \% \: 7 = 4$.<br /><br /><h2>Repeated Squaring Method</h2>Repeated Squaring Method (RSM) takes the advantage of the fact that $A^x \times A^y = A^{x+y}$.<br /><br />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.<br /><br />For example, $(19)_{10} = (10011)_2 = ( 2^4 + 2^1 + 2^0)_{10} = (16 + 2 + 1 )_{10}$<br /><br />Therefore, in equation $B^P$, we can decompose $P$ to the sum of powers of $2$.<br /><br />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$.<br /><br />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.<br /><br /><h2>Code</h2>Here is the code for RSM. The Explanation is below the code.</div><div><pre class="brush:cpp">int bigmod ( int b, int p, int m ) {<br /> int res = 1 % m, x = b % m;<br /> while ( p ) {<br /> if ( p & 1 ) res = ( res * x ) % m;<br /> x = ( x * x ) % m;<br /> p >>= 1;<br /> }<br /> return res;<br />}</pre></div><h3>Explanation</h3>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.<br /><br />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$.<br /><br />Now, from line $3$ the loop starts. This loop runs until $p$ becomes $0$. Huh? Why is that? Keep reading.<br /><br />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$.<br /><br />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$.<br /><br />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}}$.<br /><br />$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$.<br /><br />Hence, new value of $x$ is $x \times x$. We make this update at line $5$.<br /><br />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$.<br /><br />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.<br /><br /><h3>Complexity</h3>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)$.<br /><br /><h2>Conclusion</h2>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.<br /><br />The code may seem a little confusing, so feel free to ask questions.<br /><br />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.<br /><br /><h2>Resources</h2><ol><li>forthrigth48 - <a href="http://forthright48.blogspot.com/2015/08/modular-exponentiation.html">Modular Exponentiation</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/08/bit-manipulation.html">Bit Manipulation</a></li><li>algorithmist - <a href="http://www.algorithmist.com/index.php/Repeated_Squaring">Repeated Squaring</a></li><li>Wiki - <a href="https://en.wikipedia.org/wiki/Exponentiation_by_squaring">Exponentiation by Squaring</a></li></ol>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-27296711957405064362015-09-17T12:59:00.000+06:002015-09-17T13:01:38.310+06:00Euler's Theorem and Fermat's Little TheoremWe 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.<br /><hr /><h2>Euler's Theorem</h2><blockquote class="tr_bq theorem"><b>Theorem </b>- Euler's Theorem states that, if $a$ and $n$ are coprime, then $a^{\phi(n)} \equiv 1 (\text{mod n})$ - <a href="https://en.wikipedia.org/wiki/Euler%27s_theorem">Wikipedia</a></blockquote>Here $\phi(n)$ is Euler Phi Function. Read more about Phi Function on this post - <a href="http://forthright48.blogspot.com/2015/09/euler-totient-or-phi-function.html">Euler Totient or Phi Function</a>.<br /><br /><h3>Proof</h3><div>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.</div><div><br /></div><div>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$. </div><blockquote class="tr_bq lemma"><b>Lemma </b>- Set $A$ and set $B$ contains the same integers.</blockquote>We can prove the above lemma in three steps.<br /><div><ol><li><b>$A$ and $B$ has the same number of elements<br /></b>Since $B$ is simply every element of $A$ multiplied with $a$, it contains the same number of elements as $A$. This is obvious.</li><li><b>Every integer in $B$ is coprime to $n$<br /></b>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$.</li><li><b>$B$ contains distinct integers only</b><br />Suppose $B$ does not contain distinct integers, then it would mean that there is such a $b_i$ and $b_j$ such that:<br /><br />$ab_i \equiv ab_j\:(\text{mod n})$<br />$b_i \equiv b_j\:(\text{mod n})$<br /><br />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.</li></ol><div>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$.</div></div><div><br /></div><div>Now, we can easily prove Euler's Theorem.</div><div><br /></div><div>$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})$</div><div>$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}) $</div><div>$$\therefore a^{\phi(n)} \equiv 1 \:(\text{mod n})$$</div><div><br /></div><h2>Fermat's Little Theorem</h2><div>Fermat's Little Theorem is just a special case of Euler's Theorem.</div><blockquote class="tr_bq"><b>Theorem </b>- 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})$ - <a href="https://en.wikipedia.org/wiki/Fermat%27s_little_theorem">Wikipedia</a></blockquote>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.<br /><br />We can use Euler's Theorem to prove Fermat's Little Theorem.<br /><br />Let $a$ and $p$ be coprime and $p$ be prime, then using Euler's Theorem we can say that:<br /><br />$a^{\phi(p)} \equiv 1\:(\text{mod p})$ <span style="background-color: #ffd966;">(But we know that for any prime $p$, $\phi(p) = p-1$)</span><br />$a^{p-1} \equiv 1\:(\text{mod p})$<br /><br /><h2>Conclusion</h2>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.<br /><br />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.<br /><br /><h2>Reference</h2><ol><li>Wiki - <a href="https://en.wikipedia.org/wiki/Euler%27s_theorem">Euler's Theorem</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/09/euler-totient-or-phi-function.html">Euler Totient or Phi Function</a></li><li>Wiki -<a href="https://en.wikipedia.org/wiki/Fermat%27s_little_theorem"> Fermat's Little Theorem</a></li></ol>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-78255523635451854662015-09-07T14:56:00.000+06:002015-09-07T15:00:57.614+06:00Segmented Sieve of Eratosthenes<h2>Problem</h2><div>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$.</div><div><br /></div><div>For example, $A = 11$ and $B = 19$, then answer is $4$ since there are $4$ primes within that range ($11$,$13$,$17$,$19$).</div><div><br /></div><div>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$.</div><div><br /></div><div>$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. </div><div><br /></div><div>Make sure you fully understand how "Normal" <a href="http://forthright48.blogspot.com/2015/07/sieve-of-eratosthenes-generating-primes.html">Sieve of Eratosthenes</a> works.</div><div><br /></div><h2>How Normal Sieve Works</h2><div>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.</div><div><br /></div><div>Suppose we want to find all primes between $1$ to $N$.</div><div><ol><li>First we define a new variable $sqrtn = \sqrt{N}$.</li><li>We take all primes less than $sqrtn$. </li><li>For each prime $p$, we repeat the following steps.</li><ol><li>We start from $j = p \times p$.</li><li>If $j \leq N$, we mark sieve at position $j$ to be not prime.<br />Else, we break out of the loop.</li><li>We increase the value of $j$ by $p$. And go back to step step $2$ of this loop.</li></ol><li>All positions in the sieve that are not marked are prime.</li></ol><div>This is how the basic sieve works. We will now modify it to work on segments.</div></div><div><br /></div><h2>How Segmented Sieve Works</h2><div>We will perform the same steps as normal sieve but just slightly modified.</div><div><br /></div><h3>Generate Primes Less Than $\sqrt{N}$</h3><div>In the segmented sieve, what is he largest limit possible? $10^{12}$. So let $N = 10^{12}$</div><div><br /></div><div>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}$.</div><div><br /></div><div>So, using normal sieve we will first generate all primes less than $\sqrt{N} = 10^6$.</div><div><br /></div><h3>Run on Segment</h3><div>Okay, now we can start our "Segmented" Sieve. We want to find primes between $A$ and $B$. </div><div><ol><li>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.</li><li>Define a new variable $sqrtn = \sqrt{B}$.</li><li>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. <br /><br />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.</li><li>Now, we will be working with all primes less than $sqrtn$. These primes are already generated using the normal sieve.</li><li>For each prime $p$, we will repeat the following steps:</li><ol><li>We start from $j = p \times p$.</li><li>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.<br /><br />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$.<br /></li><li>If $j \leq B$, we mark sieve at position $j$ to be not prime.<br />Else, we break out of the loop.<br /><br />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.</li><li>Increase the value of $j$ by $p$. Repeat loop from step $3$.</li></ol><li>All positions in $arr$ which has not been marked are prime.</li></ol><div>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$.<br /><br /><h2>Code</h2><div>If we convert the above pseducode into C++, then it becomes something like this:</div><pre class="brush:cpp;">int arr[SIZE];<br /><br />///Returns number of primes between segment [a,b]<br />int segmentedSieve ( int a, int b ) {<br /> if ( a == 1 ) a++;<br /><br /> int sqrtn = sqrt ( b );<br /><br /> memset ( arr, 0, sizeof arr ); ///Make all index of arr 0.<br /><br /> for ( int i = 0; i < prime.size() && prime[i] <= sqrtn; i++ ) {<br /> int p = prime[i];<br /> int j = p * p;<br /><br /> ///If j is smaller than a, then shift it inside of segment [a,b]<br /> if ( j < a ) j = ( ( a + p - 1 ) / p ) * p;<br /><br /> for ( ; j <= b; j += p ) {<br /> arr[j-a] = 1; ///mark them as not prime<br /> }<br /> }<br /><br /> int res = 0;<br /> for ( int i = a; i <= b; i++ ) {<br /> ///If it is not marked, then it is a prime<br /> if ( arr[i-a] == 0 ) res++;<br /> }<br /> return res;<br />}<br /></pre>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.<br /><br />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.<br /><br /><h2>Conclusion</h2></div></div><div>I first learned about Segmented Sieve from blog of <a class="g-profile" href="https://plus.google.com/100568681605555500424" target="_blank">+Zobayer Hasan</a>. You can have a look at that post <a href="http://zobayer.blogspot.com/2009/09/segmented-sieve.html">here</a>. 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.<br /><br />Leave a comment if you face any difficulty in understanding the post.</div><div><br /></div><h2>Reference</h2><div><ol><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/07/sieve-of-eratosthenes-generating-primes.html">Sieve of Eratosthenes</a></li><li>zobayer - <a href="http://zobayer.blogspot.com/2009/09/segmented-sieve.html">Segmented Sieve</a></li></ol><h2>Related Problems</h2></div><div><ol><li><a href="http://www.spoj.com/problems/PRIME1/">SPOJ PRIME1 - Prime Generator</a></li><li><a href="http://lightoj.com/volume_showproblem.php?problem=1197">LOJ 1197 - Help Hanzo</a></li></ol></div><div><br /></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com1tag:blogger.com,1999:blog-811371107174953309.post-68915884867225433132015-09-04T17:38:00.000+06:002015-09-04T18:07:17.097+06:00Euler Totient or Phi FunctionI 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.<br /><hr /><br /><h2>Problem</h2>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$.<br /><br />For example, if $N = 10$, then there are $4$ numbers, namely $\{1,3,7,9\}$, which are coprime to $10$.<br /><br />This problem can be solved using Euler Phi Function, $\phi()$. Here is the definition from Wiki:<br /><blockquote>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. - <a href="https://en.wikipedia.org/wiki/Euler%27s_totient_function">Wiki</a></blockquote>That's exactly what we need to find in order to solve the problem above. So, how does Euler Phi work?<br /><br /><h2>Euler Phi Function</h2><div>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:</div><div>$$\phi(n) = n \times \frac{p_1-1}{p_1}\times\frac{p_2-1}{p_2}...\times\frac{p_k-1}{p_k}$$</div><div>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.</div><div><br /></div><h2>Proof of Euler Phi Function</h2><div>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.</div><div><br /></div><h3>Base Case - $\phi(1)$</h3><div>First, the base case. Phi function counts the number of <b>positive </b>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.</div><div><br /></div><div>$\phi(1) = 1$, since $1$ itself is the only number which is coprime to it.</div><div><br /></div><h3>When $n$ is a Prime - $\phi(p)$</h3><div>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$.</div><div><br /></div><h3>When $n$ is Power of Prime - $\phi(p^a)$</h3><div>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 <b>not </b>coprime.</div><div><br /></div><div>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$.</div><div><br />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} )$</div><div><br /></div><div>It's starting to look like the equation above, right?</div><div><br /></div><h3>Assuming $\phi()$ is Multiplicative - $\phi( m \times n )$</h3><div>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.</div><div><br /></div><div>So how do we prove that Euler Phi is multiplicative and how does Euler Phi being multiplicative helps us?</div><div><br /></div><div>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.</div><div><br /></div><div>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:</div><div><br /></div><div>$\phi(n) = \phi(p_1^{a_1}p_2^{a_2}...p_k^{a_k})$</div><div>$\phi(n) = \phi(p_1^{a_1}) \times \phi(p_2^{a_2}) ... \times \phi(p_k^{a_k})$.</div><div><br /></div><div>We can already calculate $\phi(p^a) = p^a \times \frac{p-1}{p}$. So our equationg becomes:</div><div><br /></div><div>$\phi(n) = \phi(p_1^{a_1}) \times \phi(p_2^{a_2}) ... \times \phi(p_k^{a_k})$</div><div>$\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}$</div><div>$\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}$</div><div>$$\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}$$</div><div>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.</div><div><br /></div><h2>Proof for Multiplicity of Euler Phi Function</h2><div>We are trying to prove the following theorem:</div><blockquote class="tr_bq theorem"><b>Theorem </b>$1$: If $m$ and $n$ are coprime, then $\phi(m \times n ) = \phi ( m ) \times \phi(n)$</blockquote>But in order to prove Theorem $1$, we will need to prove few other theorems first.<br /><br /><h3>Theorem Related to Arithmetic Progression</h3><blockquote class="tr_bq theorem"><b>Theorem </b>$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.</blockquote><div>Umm, looks like Theorem $2$ is packed with too much information. Let me break it down.<br /><br />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.<br /><br />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?<br /><br />$\{ a + 0m, a + 1m, a + 2m, a + 3m ... a + (n-1)m\}$ There are exactly $n$ terms in this list.<br /><br />Next, find their modulus by $n$. That is, find the remainder of each term after dividing by $n$.<br /><br />Now, Theorem $2$ claims that, if $m$ and $n$ are coprime, then the list will contain a permutation of $0$ to $n-1$.<br /><br />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)$.<br /><br />I hope that now it's clear what the Theorem claims. Now we will look into the proof.<br /><br /><h4><b>Proof of Theorem 2</b></h4><div>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.</div><div><ol><li><b>There are exactly $n$ elements in the list.</b> <br />Well, since we took $n$ terms from the AP, this is obvious.</li><li><b>Each element of the list has value between $0$ to $n-1$</b><br />We performed modulo operations on each element by $n$. So this is also obvious.</li><li><b>No remainder has the same value as another.</b><br />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.<br /><br />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$.<br /><br />Therefore, $( a + qm ) - ( a + pm ) \equiv 0\: ( mod\: n )$<br />$(a+qm -a - pm) \equiv 0\: ( mod\: n )$<br />$m ( q - p ) \equiv 0\: ( mod\: n )$<br /><br />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.</li></ol><div><br /></div><h3>Theorem Related to Remainder</h3></div></div><blockquote class="tr_bq theorem"><b>Theorem </b>$3$: If a number $x$ is coprime to $y$, then $(x \:\%\: y)$ will also be coprime to $y$.</blockquote><div>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.<br /><br />Theorem $3$ claims that $y$ and $r$ are coprime. What happens if this claim is false? Suppose they are <b>not </b>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.<br /><br /><h2>Proof for Multiplicity of Euler Phi Function Continued</h2></div><div>We are now ready to tackle proving Theorem $1$.</div><div><br /></div><div>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 )$. </div><div><br /></div><div>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$.</div><div><br /></div><div>Now, let us build a table of with $n$ rows and $m$ columns. Therefore, the table will look like the following:</div><div class="postTable"><table><tbody><tr><td>1</td><td>2</td><td>3</td><td>...</td><td>m</td></tr><tr><td>1+m</td><td>2+m</td><td>3+m</td><td>...</td><td>2m</td></tr><tr><td>1+2m</td><td>2+2m</td><td>3+2m</td><td>...</td><td>3m</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>1 + (n-1)m</td><td>2 + (n-1)m</td><td>3 + (n-1)m</td><td>...</td><td>mn</td></tr></tbody> </table></div>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$.<br /><br />Now, how many numbers in each column are coprime to $n$? Using theorem $2$, we know that each column has 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.<br /><br />So, there are exactly $\phi(n)$ numbers which are coprime to $n$ in each column.<br /><br />We need to find numbers that coprime to <b>both </b>$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?<br /><br />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:<br /><div class="postTable"><table><tbody><tr><td>1</td><td>2</td><td>3</td><td>...</td><td>m</td></tr><tr><td>1</td><td>2</td><td>3</td><td>...</td><td>m</td></tr><tr><td>1</td><td>2</td><td>3</td><td>...</td><td>m</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>1</td><td>2</td><td>3</td><td>...</td><td>m</td></tr></tbody></table></div>So, how many columns are there which are coprime to $m$? There are $\phi(m)$ columns which are coprime to $m$.<br /><br />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$.<br />$$\therefore \phi(m) \times \phi(n) = \phi(m \times n)$$<br /><h2>Code</h2>Since we have to factorize $n$ in order to calculate $\phi(n)$, we can modify our $factorize()$ function from post "<a href="http://forthright48.blogspot.com/2015/07/prime-factorization-of-integer.html">Prime Factorization of Integer Number</a>" to handle Euler Phi.<br /><pre class="brush:cpp;highlight:[1,2,10,11,15,16]">int eulerPhi ( int n ) {<br /> int res = n;<br /> int sqrtn = sqrt ( n );<br /> for ( int i = 0; i < prime.size() && prime[i] <= sqrtn; i++ ) {<br /> if ( n % prime[i] == 0 ) {<br /> while ( n % prime[i] == 0 ) {<br /> n /= prime[i];<br /> }<br /> sqrtn = sqrt ( n );<br /> res /= prime[i];<br /> res *= prime[i] - 1;<br /> }<br /> }<br /> if ( n != 1 ) {<br /> res /= n;<br /> res *= n - 1;<br /> }<br /> return res;<br />}<br /></pre>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.<br /><br /><h2>Conclusion</h2><div>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. </div><div><br /></div><div>Leave a comment if you face any difficulty in understanding the post.<br /><br /></div><h2>Reference</h2><ol><li>Wiki - <a href="https://en.wikipedia.org/wiki/Euler%27s_totient_function">Euler Totient Function</a></li></ol><br />Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com2tag:blogger.com,1999:blog-811371107174953309.post-73915137558748769762015-08-31T21:11:00.000+06:002015-09-02T14:58:53.233+06:00Bit ManipulationOn this post, we will look into how to use "<a href="http://forthright48.blogspot.com/2015/08/bitwise-operators.html">Bitwise Operators</a>" to manipulate bits of a number. Bit manipulation is a useful tool that can sometimes perform complicated tasks in a simple manner.<br /><br />We will work on integer values and assume each number has $32$ bits in them. If a number has only $0$ or $1$ as its digit, then assume that it is in binary form. All position of bits is $0$ indexed.<br /><br />Since binary values can only be $0$ or $1$, we sometimes refer them like light bulbs. A bit being "Off" means it has value $0$ and being "On" means it has value $1$. Bit $1$ is also referred as "Set" and bit $0$ as "Reset".<br /><br />We need to have a strong grasp of how AND, OR, Negation, XOR and Shift Operators work to understand how bit manipulation works. If you have forgotten about them, then make sure you revise.<br /><br /><h2>Checking Bit at Position X</h2><div>Given a number $N$, find the value of its bit at position $X$. For example, if $N=12=(1100)_2$ and $X=2$, then value of bit is $1$. </div><div><br /></div><div>So how do we find it? We can find this value in three steps:</div><div><ol><li>Let, $mask = 1 << X$</li><li>Let, $res = N \: \& \: mask$</li><li>If $res$ is $0$, then bit is $0$ at that position, else bit is $1$.</li></ol></div><div>Let us see it in action with the example above, $N = 12$ and $X = 2$.</div><div><br /></div><div>$1100$</div><div>$\underline{0100} \:\&$ <span style="background-color: #ffd966;">( $1 << 2$ is just $1$ shifted to left twice )</span></div><div>$0100$</div><div><br /></div><div>Now, what happened here? In step $1$, when we left shifted $1$ $X$ times, we created a number $mask$ which has bit $1$ only at position $X$. This is a useful application of left shift operator. Using this, we can move $1$ bit anywhere we want. </div><div><br /></div><div>Next, at step $2$, we find the result of performing AND operator between $N$ and $mask$. Since $mask$ has $0$ in all positions except $X$, the result will have $0$ in all places other than $X$. That's because performing AND with $0$ will always give $0$. Now, what about position $X$. Will the result at position $X$ have $0$ too? That depends on the bit of $N$ at position $X$. Since, the $X$ bit in the mask is $1$, the only way result will have $0$ at that position if $N$ has $0$ in that position.</div><div><br /></div><div>So, when all position of $res$ is $0$, that is $res$ has value $0$, we can say that $N$ has $0$ bit in position $X$, otherwise it doesn't.</div><div><br /></div><h2>Set Bit at Position $X$</h2><div>Given a value $N$, turn the bit at position $X$ of $N$ to $1$. For example, $N=12=1100$ and $X=0$, then $N$ will become $13=1101$. </div><div><br /></div><div>This can be done in two steps:</div><div><ol><li>Let, $mask = 1 << X$</li><li>$N = N \:|\: mask$</li></ol><div>$1100$</div><div>$\underline{0001}\:|$</div><div>$1101$</div><div><br /></div><div>In step $1$, we shift a $1$ bit to position $X$. Then we perform OR between $N$ and $mask$. Since $mask$ has 0 in all places except $X$, in result all those places remain unchanged. But $mask$ has $1$ in position $X$, so in result position $X$ will be $1$. </div></div><div><br /></div><div><h2>Reset Bit at Position $X$</h2></div><div><div>Given a value $N$, turn the bit at position $X$ of $N$ to $0$. For example, $N=12=1100$ and $X=3$, then $N$ will become $4=0100$. </div></div><div><br /></div><div>This cane be done in three steps:</div><div><ol><li>Let, $mask = 1 << X$</li><li>$mask = \sim mask$</li><li>$N = N \:\&\: mask$</li></ol><div>$1100$</div></div><div>$\underline{0111}\:\&$ <span style="background-color: #ffd966;">( We first got $1000$ from step $1$. Then used negation from step 2.)</span></div>$0100$<br /><div><br /></div><div>In step $1$ we move $1$ bit to position $X$. This gives us a number with $1$ in position $X$ and $0$ in all other places. Next we negate the $mask$. This flips all bits of $mask$. So now we have $0$ in position $X$ and $1$ in all other places. Now when we perform AND between $mask$ and $N$, it forces the bit at the $X$ position to $0$ and all other bits stay intact.</div><div><br /></div><h2>Toggle Bit at Position X</h2><div>Given a value $N$, toggle the bit at position $X$, i.e, if bit at $X$ is $0$ then make it $1$ and if it is already $1$, make it $0$. </div><div><br /></div><div>This can be done in two steps:</div><div><ol><li>Let, $mask = 1 << X$</li><li>$N = N \:\wedge\: mask$</li></ol><div>First we shift bit $1$ to our desired position $X$ in step $1$. When XOR is performed between a bit and $0$, the bit remains unchanged. But when XOR performed between a bit and $1$, it toggles. So all position except $X$ toggles during step 2, since all position in mask except $X$ is $0$.</div></div><div><br /></div><h2>Coding Tips</h2><div>When coding these in C++, make sure you use lots of brackets. Bitwise operators have lower precedence than $==$ and $!=$ operator which often causes trouble. See <a href="http://en.cppreference.com/w/cpp/language/operator_precedence">here</a>. </div><div><br />What would be the output of this code:</div><div><pre class="brush:cpp">if ( 0 & 6 != 7 ) printf ( "True\n" );<br />else printf ( "False\n" );<br /></pre></div><div>You might expect things to follow this order: $0 \:\&\: 6 = 0$ and $0 \:!= 7$, so "True" will be output. But due to precedence of operators, the output comes out "False". First $6\: != 7$ is checked. This is true, so $1$ is returned. Next $0 \:\&\:1$ is performed which is $0$.<br /><br />The best way to avoid these kind of mishaps is to use brackets whenever you are using bitwise operators. The correct way to write the code above is:<br /><pre class="brush:cpp">if ( (0 & 6) != 7 ) printf ( "True\n" );<br />else printf ( "False\n" );<br /></pre>This prints "True" which is our desired result.<br /><br /></div><h2>Conclusion</h2><div>These are the basics of bit manipulation and by no means the end of it. There are lots of other tricks that are yet to be seen, such as "Check Odd or Even", "Check Power of Two", "Right Most Bit". We will look into them some other day.<br /><br /><h2>Reference</h2><div><ol><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/08/bitwise-operators.html">Bitwise Operators</a></li><li>CPPReference - <a href="http://en.cppreference.com/w/cpp/language/operator_precedence">C++ Operator Precedence</a></li></ol></div><br /></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-54947300137314188832015-08-27T17:07:00.000+06:002015-08-27T17:08:48.879+06:00Bitwise OperatorsWe know about arithmetic operators. The operators $+,-,/$ and $\times$ adds, subtracts, divides and multiplies respectively. We also have another operator $\%$ which finds the modulus.<br /><div><br /></div><div>Today, we going to look at $6$ more operators called "Bitwise Operators". Why are they called "Bitwise Operators"? That's because they work using the binary numerals (bits, which are the individual digits of the binary number) of their operands. Why do we have such operators? That's because in computers all information is stored as strings of bits, that is, binary numbers. Having operators that work directly on them is pretty useful.</div><div><br /></div><div>We need to have a good idea how Binary Number System works in order to understand how these operators work. Read more on number system in [1]<a href="http://forthright48.blogspot.com/2015/08/introduction-to-number-systems.html">Introduction to Number Systems</a>. Use the $decimalToBase()$ function from [1] to convert the decimal numbers to binary and see how they are affected.</div><div><br /></div><h2>Bitwise AND ($\&$) Operator</h2><div>The $\&$ operator is a binary operator. It takes two operands and returns single integer as the result. Here is how it affects the bits.</div><div><br /></div><div>$0 \: \& \: 0 = 0$</div><div>$0 \: \& \: 1 = 0$</div><div>$1 \: \& \: 0 = 0$</div><div>$1 \: \& \: 1 = 1$</div><div><br /></div><div>It takes two bits and returns another bit. The $\&$ operator will take two bits $a, b$ and return $1$ only if both $a$ <b>AND </b>$b$ are $1$. Otherwise, it will return $0$.</div><div><br /></div><div>But that's only for bits. What happens when we perform $\&$ operation on two integer number?</div><div><br /></div><div>For example, what is the result of $A \: \& \: B$ when $A = 12$ and $ B =10$?</div><div><br /></div><div>Since $\&$ operator works on bits of binary numbers we have to convert $A$ and $B$ to binary numbers.</div><div><br /></div><div>$A = (12)_{10} = (1100)_2$</div><div>$B = (10)_{10} = (1010)_2$</div><div><br /></div><div>We know how to perform $\&$ on individual bits, but how do we perform $\&$ operation on strings of bits? Simple, take each position of the string and perform and operations using bits of that position.</div><div><br /></div><div>$1100$</div><div>$\underline{1010}\: \&$</div><div>$1000$</div><div><br /></div><div>Therefore, $12 \: \& \: 10 = 8$.</div><div><br /></div><div>In C++, it's equivalent code is:</div><pre class="brush:cpp;">printf ("%d\n", 12 & 10 );</pre><div><br /><div></div><h2>Bitwise OR ($|$) Operator</h2><div>The $|$ operator is also a binary operator. It takes two operands and returns single integer as the result. Here is how it affects the bits:<br /><br /></div><div></div><div><div>$0 \: | \: 0 = 0$</div><div>$0 \: | \: 1 = 1$</div><div>$1 \: | \: 0 = 1$</div><div>$1 \: | \: 1 = 1$<br /><br /></div></div><div></div><div>The $|$ operator takes two bits $a, b$ and return $1$ if $a$ <b>OR</b> $b$ is $1$. Therefore, it return $0$ only when both $a$ and $b$ are $0$.<br /><br /></div><div></div><div>What is the value of $A\:|\:B$ if $A=12$ and $B=10$? Same as before, convert them into binary numbers and apply $|$ operator on both bits of each position.<br /><br /></div><div></div><div><div>$1100$</div><div>$\underline{1010}\: |$</div><div>$1110$<br /><br /></div></div><div></div><div>Therefore $12\:|\:10 = 14$.</div><pre class="brush:cpp;">printf ( "%d\n", 12 | 10 );</pre><div><br /></div><h2>Bitwise XOR ($\wedge$) Operator</h2><div>Another binary operator that takes two integers as operands and returns integer. Here is how it affects two bits:<br /><br /></div><div></div><div><div>$0 \: \wedge \: 0 = 0$</div><div>$0 \: \wedge \: 1 = 1$</div><div>$1 \: \wedge \: 0 = 1$</div><div>$1 \: \wedge \: 1 = 0$<br /><br /></div></div><div></div><div>XOR stands for Exclusive-OR. This operator returns $1$ only when both operand bits are not same. Otherwise, it returns $0$.<br /><br /></div><div></div><div>What is the value of $A\:\wedge\:B$ if $A=12$ and $B=10$?<br /><br /></div><div></div><div><div>$1100$</div><div>$\underline{1010}\: \wedge$</div><div>$0110$<br /><br /></div></div><div></div><div>Therefore, $12\:\wedge\:10 = 6$.<br /><br /></div><div></div><div>In mathematics, XOR is represented with $\oplus $, but I used $\wedge$ cause in C++ XOR is performed with $\wedge$.</div><pre>printf ( "%d\n", 12 ^ 10 );</pre><br /><h2>Bitwise Negation ($\sim$) Operator</h2>This is a unary operator. It works on a single integer and flips all its bits. Here is how it affects individual bits:<br /><br />$\sim\:0 = 1$<br />$\sim\: 1 = 0$<br /><br />What is the value of $\sim \: A$ if $A = 12$?<br /><br />$\sim \: (1100)_2 = (0011)_2 = (3)_{10}$<br /><br />But this will not work in code cause $12$ in C++ is not $1100$, it is $0000...1100$. Each integer is $32$ bits long. So when each of the bits of the integer is flipped it becomes $11111...0011$. If you don't take unsigned int, the value will even come out negative.<br /><pre class="brush:cpp">printf ( "%d\n", ~12 );</pre><br /><h2>Bitwise Left Shift ($<<$) Operator</h2>The left shift operator is a binary operator. It takes two integers $a$ and $b$ and shifts the bits of $a$ towards <b>LEFT</b> $b$ times and adds $b$ zeroes to the end of $a$ in its binary system.<br /><br />For example, $(13)_{10} << 3 = (1101)_2 << 3 = (110100)_2$.<br /><br />Shifting the bits of a number $A$ left once is same as multiplying it by $2$. Shifting it left three times is same as multiplying the number by $2^3$.<br /><br />Therefore, the value of $A << B = A \times 2^B$.<br /><pre class="brush:cpp">printf ( "%d\n", 1 << 3 );</pre><br /><h2>Bitwise Right Shift ($>>$) Operator</h2>The $>>$ Operator does opposite of $<<$ operator. It takes two integer $a$ and $b$ and shifts the bits of $a$ towards <b>RIGHT </b>$b$ times. The rightmost $b$ bits are lost and $b$ zeroes are added to the left end.<br /><br />For example, $(13)_{10} >> 3 = (1101)_2 >> 3 = (1)_2$.<br /><br />Shifting the bits of a number $A$ right once is same as dividing it by $2$. Shifting it right three times is same as dividing the number by $2^3$.<br /><br />Therefore, the value of $A >> B = \lfloor \frac{A}{2^B} \rfloor$.<br /><pre class="brush:cpp;">printf ( "%d\n", 31 >> 3 );</pre><br /><h2>Tips and Tricks</h2><div><ol><li>When using $<<$ operator, careful about overflow. If $A << B$ does not fit into $int$, make sure you type cast $A$ into $long\:long$. Typecasting $B$ into $long\:long$ does not work.</li><li>$A \:\&\:B \leq MIN(A,B)$ </li><li>$A\:|\:B \geq MAX(A,B)$</li></ol><div><br /></div></div><hr /><br />That's all about bitwise operations. These operators will come useful during "Bits Manipulation". We will look into it on next post.<br /><br /><h2>Resource</h2><div><ol><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/08/introduction-to-number-systems.html">Introduction to Number Systems</a></li><li>Wiki - <a href="https://en.wikipedia.org/wiki/Bitwise_operation">Bitwise Operation</a></li></ol></div></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com2tag:blogger.com,1999:blog-811371107174953309.post-19454745132883150242015-08-27T12:46:00.000+06:002015-08-27T12:50:27.033+06:00SPOJ LCMSUM - LCM Sum<h2>Problem</h2><div>Problem Link - <a href="http://www.spoj.com/problems/LCMSUM/">SPOJ LCMSUM</a></div><div><br /></div><div>Given $n$, calculate the sum $LCM(1,n) + LCM(2,n) + ... + LCM(n,n)$, where $LCM(i,n)$ denotes the Least Common Multiple of the integers $i$ and $n$.</div><div><br /></div><h2>Solution</h2>I recently solved this problem and found the solution really interesting. You can find the formula that solves this problem on <a href="https://oeis.org/A051193">OEIS</a>. I will show you the derivation here.<br /><div><br /></div><div>In order to solve this problem, you need to know about Euler Phi Function, finding Divisor using Sieve and some properties of LCM and GCD.</div><div><br /></div><h3>Define SUM</h3><div>Let us define a variable SUM which we need to find.</div><div><br /></div><div>$SUM = LCM(1,n) + LCM(2,n) + ... + LCM(n,n)$ <span style="background-color: #ffd966;">(take $LCM(n,n)$ to the other side for now )</span></div>$SUM - LCM(n,n) = LCM(1,n) + LCM(2,n) + ... + LCM(n-1,n)$<br /><div><br /></div><div>We know that $LCM(n,n) = n$.</div><div><br /></div><div>$SUM - n = LCM(1,n) + LCM(2,n) + ... + LCM(n-1,n)$ <span style="background-color: #ffd966;">($eq 1$)</span></div><div><br /><h3>Reverse and Add</h3><div>$SUM - n = LCM(1,n) + LCM(2,n) + ... + LCM(n-1,n)$ <span style="background-color: #ffd966;">(Reverse $eq 1$ to get $eq 2$)</span></div><div>$SUM - n = LCM(n-1,n) + LCM(n-2,n) + ... + LCM(1,n)$ <span style="background-color: #ffd966;">( $eq 2$ )</span></div><br />No what will we get if we add $eq 1$ with $eq 2$? Well, we need to do more work to find that.<br /><div><br /></div><h3>Sum of $LCM(a,n) + LCM(n-a,n)$</h3><div>$x = LCM(a,n) + LCM(n-a,n)$</div><div>$x = \frac{an}{gcd(a,n)} + \frac{(n-a)n}{gcd(n-a,n)}$ <span style="background-color: #ffd966;">($eq 3$)</span></div><div><br /></div><div>Arghh. Now we need to prove that $gcd(a,n)$ is equal to $gcd(n-a,n)$.<br /><br /></div><div>If $c$ divides $a$ and $b$, then, $c$ will divide $a+b$ and $a-b$. This is common property of division. So, if $g = gcd(a,n)$ divides $a$ and $n$, then it will also divide $n-a$. Hence, $gcd(a,n)=gcd(n-a,n)$.<br /><br />So $eq 3$ becomes:<br /><br />$x = \frac{an}{gcd(a,n)} + \frac{(n-a)n}{gcd(a,n)}$<br />$x = \frac{an + n^2 -an}{gcd(a,n)}$.<br />$x = \frac{n^2}{gcd(a,n)}$.<br /><br />Now, we can continue adding $eq 1$ with $eq 2$.<br /><br /></div><h3>Reverse and Add Continued</h3><div><div><div>$SUM - n = LCM(1,n) + LCM(2,n) + ... + LCM(n-1,n)$ <span style="background-color: #ffd966;">($eq 1$)</span></div><div>$SUM - n = LCM(n-1,n) + LCM(n-2,n) + ... + LCM(1,n)$ <span style="background-color: #ffd966;">( $eq 2$ Add them )</span></div></div>$2(SUM-n)= \frac{n^2}{gcd(1,n)} + \frac{n^2}{gcd(2,n)} + ... \frac{n^2}{gcd(n-1,n)}$</div><div>$2(SUM-n ) = \sum_{i=1}^{n-1}\frac{n^2}{gcd(i,n)}$</div><div>$2(SUM-n ) = n\sum_{i=1}^{n-1}\frac{n}{gcd(i,n)}$ <span style="background-color: #ffd966;">( take $n$ common )</span></div><div><br /></div><h3>Group Similar GCD</h3><div>What are the possible values of $g = gcd(i,n)$? Since $g$ must divide $n$, $g$ needs to be a divisor of $n$. So we can list the possible values of $gcd(i,n)$ by finding the divisors of $n$.</div><div><br /></div><div>Let $Z = n\sum_{i=1}^{n-1}\frac{n}{gcd(i,n)}$. Now, every time $gcd(i,n) = d$, where $d$ is a divisior of $n$, $n \times \frac{n}{d}$ is added to $Z$. Therefore, we just need to find, for each divisor $d$, how many times $n \times \frac{n}{d}$ is added to $Z$.</div><div><br /></div><div>How many values $i$ can we put such that $gcd(i,n) = d$? There are $\phi(\frac{n}{d})$ possible values of $i$ for which we get $gcd(i,n)=d$. Therefore:</div><div><br /></div><div>$2(SUM-n ) = n\sum_{i=1}^{n-1}\frac{n}{gcd(i,n)}$</div><div>$2(SUM-n ) = n\sum_{d\:|\:n, d \neq n } \phi(\frac{n}{d}) \times \frac{n}{d}$ <span style="background-color: #ffd966;">(but, $\frac{n}{d}$ is also a divisor of $n$ )</span></div><div>$2(SUM-n ) = n\sum_{d\:|\:n, d \neq 1 } \phi(d) \times d$ <span style="background-color: #ffd966;">( when $d = 1$, we get $\phi(1)\times 1$ )</span> </div><div>$2(SUM-n ) = n( \sum_{d\:|\:n } ( \phi(d) \times d ) -1 )$ </div><div>$2(SUM-n ) = n \sum_{d\:|\:n } (\phi(d) \times d ) - n$</div><div>$2SUM - 2n + n = n\sum_{d\:|\:n } \phi(d) \times d $</div><div>$2SUM - n = n\sum_{d\:|\:n } \phi(d) \times d $<br />$$2SUM = n( \sum_{d\:|\:n } \phi(d) \times d )+ n \\2SUM = n ( \sum_{d\:|\:n } ( \phi(d) \times d ) + 1 )\\ \therefore SUM = \frac{n}{2} ( \sum_{d\:|\:n } ( \phi(d) \times d ) + 1 )$$<br />Using this formula we can solve the problem.</div><div><br /></div><h2>Code</h2></div><pre class="brush: cpp;">#include <bits/stdc++.h><br /><br />#define FOR(i,x,y) for(vlong i = (x) ; i <= (y) ; ++i)<br /><br />using namespace std;<br />typedef long long vlong;<br /><br />vlong res[1000010];<br />vlong phi[1000010];<br /><br />void precal( int n ) {<br /> ///Calculate phi from 1 to n using sieve<br /> FOR(i,1,n) phi[i] = i;<br /> FOR(i,2,n) {<br /> if ( phi[i] == i ) {<br /> for ( int j = i; j <= n; j += i ) {<br /> phi[j] /= i;<br /> phi[j] *= i - 1;<br /> }<br /> }<br /> }<br /><br /> ///Calculate partial result using sieve<br /> ///For each divisor d of n, add phi(d)*d to result array<br /> FOR(i,1,n){<br /> for ( int j = i; j <= n; j += i ) {<br /> res[j] += ( i * phi[i] );<br /> }<br /> }<br />}<br /><br />int main () {<br /> precal( 1000000 );<br /><br /> int kase;<br /> scanf ( "%d", &kase );<br /><br /> while ( kase-- ) {<br /> vlong n;<br /> scanf ( "%lld", &n );<br /><br /> ///We already have partial result in res[n]<br /> vlong ans = res[n] + 1;<br /> ans *= n;<br /> ans /= 2;<br /><br /> printf ( "%lld\n", ans );<br /> }<br /><br /> return 0;<br />}<br /></pre>We nee to precalculate partial result using a sieve for all values of $N$. Precalculation has a complexity of $O(N\: lnN)$. After pre-calculation, for each $N$ we can answer in $O(1)$.<br /><br /><h2>Reference</h2><div><ol><li>OEIS - <a href="https://oeis.org/A051193">A051193</a></li></ol></div><br />Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com8tag:blogger.com,1999:blog-811371107174953309.post-49413725002649047392015-08-25T22:51:00.000+06:002015-08-25T22:51:59.047+06:00Contest Analysis: BUET Inter-University Programming Contest - 2011, UVa 12424 - 12432Last week (2015-08-21) we practiced with this set. Some of the problems didn't have proper constraints, which caused some headaches. This contest originally followed Google Code Jam format with two subproblems ( small and large ) for each problem. Each of the sub-problems had a different constraint on them. But when they uploaded it online, they removed both constraints from some problems, making things confusing. We ended up guessing the constraints when solving them.<br /><br />During contest time, we managed to solve $B, E$, and $H$. I was stuck with $F$. Should have moved on. $C$ and $D$ were far easier than $F$. Bad decision from my end.<br /><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=3855">UVa 12424 - Answering Queries on a Tree</a></h2><div><b>Complexity: $O(logN)$</b><br /><b>Category: Segment Tree, DFS, LCA</b><br /><br />At first it seemed like a problem on Heavy Light Decomposition. But after a moment I realized it can be done without HLD.<br /><br />Run a dfs from any node and make the tree rooted tree. Now consider the color $X$. For every node that have color $X$, the edge from its parent will have cost $1$. Since the root can also have color $X$, we need to add a dummy node $0$ and make that root. Now, run another dfs and calculate the distance from the root to each node. We repeat this for all $10$ color and build $10$ trees. Let us call tree with color of node $X$, $T_X$.<br /><br />Now, whenever we have an update to change node $u$ from color $X$ to color $Y$, it means that the edge entering $u$ in tree $T_X$ will now have a cost $0$. Since its cost decreased, the distance from root of all the nodes in the subtree of $u$ will be lowered by $1$. We can apply this change using preorder traversal dfs + segment tree. Do the opposite for $T_Y$. Increase the value of all nodes under subtree $u$ in $T_Y$.<br /><br />And for query $u$ and $v$, find the distance between these two nodes in each of the $10$ trees. This can be done in $O(logN)$ using LCA + Segment Tree.<br /><br />Code: <a href="http://ideone.com/UEU9J4">http://ideone.com/UEU9J4</a><br /><br /></div><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=3856">UVa 12425 - Best Friend</a></h2><div><b>Complexity: $O(\sqrt{N} + \sqrt[3]{N} \times \text{Complexity of } \phi (N) )$</b></div><div><b>Category: Number Theory, Euler Phi, GCD</b></div><div><br /></div><div>So for given integers $N$ and $X$, we have to find out how many number $I$ are there such that $gcd(N, I) \leq X$.</div><div><br /></div><div>In order to calculate $gcd(N,I) \leq X$, I first need to be able to calculate how many numbers are there such that $gcd(N,I) = Y$. I started with a small value of $Y$ first. So the first thing I asked myself, how many numbers are there such that $gcd(N, I) = 1$? The answer is $\phi (N)$. Then I asked how many are there such that $gcd(N, I) = 2$? </div><div><br /></div><div>This took a moment, but I soon realized, if $gcd(N,I)$ is equal $2$, then if I divide both $N$ and $I$ by $2$, then $gcd(\frac{N}{2},\frac{I}{2})$ will be equal to $1$. Now, all I had to calculate how many $I$ are there such that $gcd(\frac{N}{2},\frac{I}{2}) = 1$. The answer is $\phi (\frac{N}{2})$.</div><div><br /></div><div>Therefore, there are $\phi (\frac{N}{Y})$ number of $I$ such that $gcd(N,I) = Y$.</div><div><br /></div><div>Great. Now we can calculate number of $I$ such that $gcd(N,I) = Y$. Now all we need to do is find the sum of all number of $I$ for which $gcd(N, I) = Y$ and $Y \leq X$. GCD of $N$ and $I$ can be as high as $10^{12}$, so running a loop won't do.</div><div><br /></div><div>Next I thought, what are the possible values of $g = gcd(N,I)$? Since $g$ is a number that divides $N$ and $I$, $g$ must be a divisor of $N$. So I simply calculated all the divisors of $N$ using a $\sqrt{N}$ loop. Next for each divisor $d$, I calculated $\phi ( \frac{N}{d})$ which is the number of ways to chose $I$ such that $gcd(N,I) = d$.</div><div><br />Now, for each query $X$, all I needed to do was find sum of $\phi ( \frac{N}{d} ) $, where $d$ is a divisor of $N$ and $d \leq X$. Since all values of $\phi ( \frac{N}{d} ) $ was precalculated, using cumalitive sum I answered it in $O(1)$.</div><div><br /></div><div>Code: <a href="http://ideone.com/WrIH0C">http://ideone.com/WrIH0C</a></div><div><br /></div><div><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=3857">UVa 12426 - Counting Triangles</a></h2></div><div><b>Complexity: $O(N^2\:logN)$</b></div><div><b>Category: Two Pointer, Geo, Binary Search</b></div><div><br /></div><div>We need to choose three points of a convex hull such that the area does not exceed $K$. </div><div><br /></div><div>First, we need to know, for two points $i$ and $j$, $j>i$, what is the point $m$ such that $(i,j,m)$ forms the largest triangle possible. This can be done by selecting two points with two loops and then using ternary search. But I used two pointer technique to find $m$ for every pair of $(i,j)$ in $O(N)$.</div><div><br /></div><div>After that, for ever pair of $(i,j)$ I used binary search twice. Once on points $[j+1,m]$ and once on points $[m+1,n]$. Using binary search I found the point $x$ which gives the highest area that does not exceed $K$. Then I added $x-j-1$ in first BS and $n-1-m$ in second case.<br /><br />Careful about not adding a degenerate triangle. Hanlde the edge cases properly.</div><div><br /></div><div>Code: <a href="http://ideone.com/JtJjt9">http://ideone.com/JtJjt9</a></div><div><br /></div><div><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=3858">UVa 12427 - Donkey of the Sultan</a></h2></div><div><b>Complexity: $O(1)$</b></div><div><b>Category: Game Theory, Combinatorics, Nim Game, Bogus Nim</b></div><div><b><br /></b></div><div>This game can be converted to into a game of pile. Let the distance between the left boundary and gold coin be $x$ and the distance between diamond and silver coin be $y$. Now, instead of playing on the board, we can play with two piles of height $x$ (first pile) and $y$ (second pile).</div><div><br /></div><div>Moving gold coin left is same as removing one coin from the first pile. Moving diamond coin left is same as increasing second pile by one. Moving silver coin left is same as removing one coin from the second pile. Moving gold and silver coin left is same as removing one coin from both first and second pile.</div><div><br /></div><div>Now, note that is a state of $(x,y)$ is losing and I try to increase the second pile by one, my opponent then can simply take one from second pile and put me back to same situation. That is, increasing the second pile has no effect in this game. This is a "Bogus" move. We will ignore this move.</div><div><br /></div><div>Next I drew a $5 \times 5$ table and calculated the Win or Lose state of each possible $(x,y)$. I found that only when both $x$ and $y$ are even, the first person loses. That is the second person, me, wins the game when $x$ and $y$ both are even.</div><div><br /></div><div>What are the possible values of $x$? $x$ can only be between $a_1$ and $a_2$. Find number of even numbers in this range. Let this be $r_1$. Next, $y$ can be between $c_1$ and $c_2$. Let the number of even numbers in this range be $r_2$. The distance between gold and diamond has no effect in this game. So we can put any amount of gap there.</div><div><br /></div><div>Therefore, number of possible combination is $r_1 \times r_2 \times (b_2 - b_1 + 1 )$.</div><div><br /></div><div>Code: <a href="http://ideone.com/6q2LZT">http://ideone.com/6q2LZT</a><br /><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=3859">UVa 12428 - Enemy at the Gates</a></h2></div><div><b>Complexity: $O(\sqrt{M})$</b></div><div><b>Category: Adhoc</b></div><div><br /></div><div>This was the easiest problem in the set. In this problem, we are given $N$ nodes and $M$ edges, we have to find out how many leaves can this graph have if we maximize leaves.</div><div><br /></div><div>The problem says the graph will be connected. So I created a <a href="https://en.wikipedia.org/wiki/Star_(graph_theory)">star-shaped graph</a> with it using $N-1$ edges. Let the center node be $0$ and remaining nodes be the numbers from $1$ to $N-1$. This gives me $N-1$ critical edges. Then I checked if I still have more edges left? If I do, then I need to sacrifice a critical edge.</div><div><br /></div><div>The first edge that gets sacrifice is a special case. This is because, no matter which two nodes you connect, it will reduce the number of critical edge by $2$. There is no helping it. So let's connect $1$ and $2$ and reduce our edge by $1$ and critical edge by $2$. </div><div><br /></div><div>Then if we still have more edges left, then we need to sacrifice another critical length. From now on, at each step only one critical edge will be removed. Also, this time we can add $2$ edges to the graph by removing one critical edge. Connect an edge between $(3,1)$ and $(3,2)$. Reduce $M$ by $2$ and critical edge by $1$.</div><div><br /></div><div>Next is node $4$. We can add $3$ edges by removing $1$ critical edge now. We need to keep on doing this until all edges finish.</div><div><br /></div><div>I guess it is possible to solve the problem in $O(1)$ but I didn't bother with it. Number of test case was small anyway.</div><div><br /></div><div>Code: <a href="http://ideone.com/X3bWyZ">http://ideone.com/X3bWyZ</a><br /><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&category=&problem=3860&mosmsg=Submission+received+with+ID+15990492">UVa 12429 - Finding Magic Triplets</a></h2></div><div><b>Complexity: $O(N\:logN)$</b></div><div><b>Category: Binary Indexed Tree, Number Theory</b></div><div><br /></div><div>We have to find the number of triplets such that $a + b^2 \equiv c^3$ and $a \leq b \leq c$. Here is what I did.</div><div><br /></div><div>Iterate over $c$ from $N$ to $1$ using a loop. Let the loop iterator be $i$. Now for each $i$, we do the following:</div><div><ol><li>Find $x = ( i \times i \times i ) \:\%\: k$ and store $x$ in a BIT. BIT stores all occurrences of $c$ and by iterating in reverse order we make sure that all $c$ in BIT is greater or equal to $i$.<br /><br /></li><li>Let $b = ( i \times i ) \: \% \: k$. This is the current $b$ we are working with. In BIT we have $c$ which are $c \geq b$. Now all we need is to find possible values of $a$.<br /><br /></li><li>Now, take any value of $c$ from BIT. For this $c$ and current $b$, how many ways can we choose $a$ so that $a + b \equiv c$. Notice that $a$ can be anything from $1$ to $k$, or $0$ to $k-1$. Therefore, no matter what is the value of $c$ we can chose $1$ number from $0$ to $k-1$ so that $a + b \equiv c$.<br /><br />Let $y = \frac{i}{k}$. This means we have $y$ segments, in which value of $a$ spans from $0$ to $k-1$. So we add $y$ for each possible $c$ we take from BIT. Let $total$ be number of $c$ we have inserted till now. So we add $total \times y$ to $result$.<br /><br /></li><li>But it's not over yet. What about the values of $a$ that we did not use. To be more exact, let $r = i \: \% \: k$. If $ r > 0$, then we have unused values of $a$. We need to use these.<br /><br />But we have to be careful. We have used up values of $a$ from $1$ to $y\times k$. So the remaining $a$ we have is $1$ to $r$. $0$ is not included.<br /><br />Since we don't have all values from $0$ to $k-1$, we are not able to handle all possible values $c$ any more. How many can we handle exactly?<br /><br />$a + b \equiv c$<br />$a \equiv c - b$<br /><br />But we need $a$ can only be between $1$ and $r$.<br /><br />$1 \leq c - b \leq r$<br />$\therefore 1+b \leq c \leq r + b$<br /><br />Find out how many $c$ are there with values between $1+b$ and $r+b$ from BIT. Let this value be $z$. For each of those $c$ we can use one value of $a$ between $1$ to $r$. So add $z$ to $result$.</li></ol><div>And that's it. Be careful when handling edge cases in BIT.<br /><br />Code: <a href="http://ideone.com/iMBUmU">http://ideone.com/iMBUmU</a><br /><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&category=&problem=3861&mosmsg=Submission+received+with+ID+15990757">UVa 12430 - Grand Wedding</a></h2></div></div><div><b>Complexity: $O(N \:logN)$</b></div><div><b>Category: DFS, Bicoloring, Binary Search</b></div><div><br /></div><div>Constraint was missing. I assumed $N <= 10^6$.</div><div><br /></div><div>We have to remove some edges and then assign guards such that no edge is left unguarded and no edge has guards in both ends. We have to output the maximize the value of the edges that we remove.</div><div><br /></div><div>Binary search over the value of edges which we need to remove. Now, let's say we removed all edges with value greater or equal to $x$. How do we decide that a valid assignment of guard is possible in this problem? </div><div><br /></div><div>Suppose all the nodes in the graph are colored white. If we assign guard in a node, then we color it black. Now, two white nodes cannot be adjacent cause that would mean the edge is left unguarded. Two black nodes cannot be adjacent otherwise guards will chat with each other. That is both ends of an edge cannot have same color in the graph. This is possible when the graph has no cycle of odd length. If a graph has no cycle, then it is bicolorable. So all we need to do is check if the graph is bicolorable.</div><div><br /></div><div>Now, two special cases. If we can bicolor the graph without removing any edge, then answer is $0$. If we have to remove all edges then answer is $-1$. Why? Cause in the problem it is said we can remove <b>proper subset</b> of edges. The full set is not a proper subset.</div><div><br /></div><div>Code: <a href="http://ideone.com/R6lkje">http://ideone.com/R6lkje</a></div><div><br /></div><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=3862">UVa 12431 - Happy 10/9 Day</a></h2><div><b>Complexity: $O(\:(logN)^2\:)$</b></div><div><b>Category: Divide and Conquor, Modular Exponentiation, Repeated Squaring</b></div><div><br /></div><div>We have to find the value of:</div><div><br /></div><div>$ ( d\times b^0 + d\times b^1 + d \times b^2 + ... + d \times b^{n-1} ) \: \%\: M$</div><div>$ ( d ( b^0 + b^1 + b^2 + ... + b^{n-1} ) ) \: \%\: M$.</div><div><br /></div><div>Therefore, the problem boils down to finding the value of $( b^0 + b^1 + b^2 + ... + b^{n-1} ) \: \%\: M$. We cannot use the formula of geometric progression since $M$ and $b-1$ may not be coprime. So what do we do? We can use divide and conquor. How? Let me show you an example.</div><div><br /></div><div>Suppose $f(n) = (b^0 + b^1 + b^2 + ... + b^n)$. Let us find $f(5)$.<br /></div><div>$f(5) = b^0 + b^1 + b^2 + b^3 + b^4 + b^5$</div><div>$f(5) = (b^0 + b^1 + b^2) + b^3 ( b^0 + b^1 + b^2)$</div><div>$f(5) = f(2) + b^3 \times f(2)$</div><div>$f(5) = f(2) \times (1 + b^3)$.</div><div><br /></div><div>From this we can design a D&C in following way:</div><div><br /></div><div>$\text{n is odd: }f(n) = f(\frac{n}{2}) \times ( 1 + b^{n/2 + 1 })$</div><div>$\text{n is even: }f(n) = f(n-1) + b^n$</div><div>$f(0) = 1$</div><div><br /></div><div>Just apply modular arithmatic with the whole things. Also note that $M$ can be as larg as $10^{12}$, so $(a \times b)$ will overflow if they are multiplied directly, even after making them small by modulus $M$. So use repeated squaring metho to multiply them.</div><div><br /></div><div>Code: <a href="http://ideone.com/Fn6vzW">http://ideone.com/Fn6vzW</a></div><div><br /></div><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=3863">UVa 12432 - Inked Carpets</a></h2><div><span style="background-color: red; color: white;">Complexity:</span></div><div><span style="background-color: red; color: white;">Category:</span></div><div><br /></div><div>I didn't try this problem yet. I guess this was the stopper.</div><div><br /></div><div><hr /></div><div>It was a good set. I espcially enjoyed solving $B, D$ and $F$. </div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com2tag:blogger.com,1999:blog-811371107174953309.post-55938132127499877582015-08-20T22:44:00.000+06:002015-08-22T17:47:05.996+06:00Contest Analysis: IUT 6th National ICT Fest 2014, UVa 12830 - 12839I participated in this contest last year. Our team "NSU Shinobis" managed to solve 5 and placed in top 10 ( failed to find the rank list ). The IUT 7th ICT National Fest 2015 is going to be held next month, so I thought I would try to up solve the remaining problems. I tried to compile detailed hints for the problems here.<br /><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4695">UVa 12830 - A Football Stadium</a></h2><b>Complexity: $O(N^3)$</b><br /><b>Category: Loop, DP, Kadane's Algorithm</b><br /><b><br /></b>Given $N$ points, we have to find the largest area of the rectangle which can fit such that there is no point inside the rectangle. Points on boundaries of the rectangle are allowed.<br /><br />Now imagine a rectangle which has no point inside it and also no point on its boundary. Now focus on the top boundary of that rectangle. Since there is no point on that boundary, we can extend it up until it hits a point or limit. Same can be said for the bottom, left and right boundary. That is, the largest rectangle will have its boundary aligned with the limits or points.<br /><br />Now let us fix the top and lower boundary of the rectangle. How many possible values can they take? Each boundary can take values of $y$-coordinate from one of the $N$ points or the limits $0$ and $W$. Using two nested loops, we can fix them.<br /><br />Next we need to fix the left and right boundary. Lets us sort the points according to $x$ first and then according to $y$. Next we iterate over the points. Initially set the left boundary of rectangle aligned to $y$-axis, that is aligned with the left boundary of Sand Kingdom.<br /><br />Now, for each point, if the $y$ coordinate of the point is above or below or on the boundary we fixed then we ignore them, cause they don't cause any problem. But if they fall inside, then we need to process this point. We put the right boundary aligned to the $x$ coordinate of this point and calculate this rectangle.<br /><br />But it's not over. We shift the left boundary over to current right boundary and continue process remaining points.<br /><br />This is a classical problem of Kadane's Algorithm. With two nested loops to fix upper and lower boundary, and another loop inside iterating over $N$ points makes the complexity $O(N^3)$.<br /><br />$O(N^2)$ solution is possible, by constructing a histogram for each row ( $O(N)$ ) and the calculating area of histogram in $O(N)$. But that's not necessary here.<br /><br />Code: <a href="http://ideone.com/FjvTN2">http://ideone.com/FjvTN2</a><br /><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4696">UVa 12831 - Bob the Builder</a></h2><b>Complexity: $O(\sqrt{V}E)$</b><br /><b>Category: Max Flow, Bipartite Matching, Dinic's Algorithm, Path Cover</b><br /><br />Another classical problem, though we had a hard time solving this one. Mainly cause we didn't realize this was a path cover problem which can be solved using Max Flow.<br /><br />Take each number provided and generate all child. Consider each number a node and if it possible to generate $B$ from $A$, then add a directed edge from $A$ to $B$. When you are done, you will have DAG.<br /><br />Now split all the nodes in two. We are going to create a bipartite graph now. On the left side, we will have all the original nodes and on the right side we will have the copy we got by splitting the nodes. Now, for each edge between $A$ and $B$, add edge in the bipartite graph from $A$ on the left side to $B$ on the right side.<br /><br />Find the matching ( or maximum flow ) of this graph by running Dinic's Algorithm. The answer will be $\text{Total Nodes} - \text{Matching}$.<br /><br />Code: <a href="http://ideone.com/U1VPNB">http://ideone.com/U1VPNB</a><br /><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4697">UVa 12832 - Chicken Lover</a></h2><b>Complexity: $O(M)$</b><br /><b>Category: Expected Value, Probability, DP, Combination</b><br /><br />If a shop can make $N$ different items, and in a single day prepares $K$ items from those $N$ items, then how many different sets of menus ($total$) can they make? $total = C^N_K$. Now, if they decide to make chicken that day for sure, how many sets of menus ( $chicken$ ) can they make now? $chicken = C^{N-1}_{K-1}$. So what is the probability $P$ that if I visit a shop I will get to eat chicken? $P = \frac{chicken}{total}$. And what is the probability that I will not eat chicken? $Q = 1- P$.<br /><br />So now I know the probability of eating chicken for each shop. How do we find the expected value? We will find it using dynamic programming.<br /><br />The states of dp only be the shop number. $dp(x)$ will give me the expected number of chicken that I can eat if I visit all shops starting from $x$. Our result will be $dp(1)$.<br /><br />At each state, I have $P_i$ probability that I will eat chicken and $Q$ probability that I will not. So result for each state will be:<br /><br />$dp ( pos ) = P \times ( 1 + dp ( pos + 1 ) + Q \times dp ( pos + 1 )$<br />$dp ( pos ) = P + P \times dp ( pos + 1 ) + Q \times dp ( pos + 1 ) $<br />$dp ( pos ) = P + dp ( pos + 1 ) \times ( P + Q )$<br />$dp ( pos ) = P + dp ( pos + 1 )$.<br /><br />In order to print the result in $\frac{A}{B}$ form, we need to avoid $double$ and use integer arithmetic in all places. I implemented my own fraction class for that purpose.<br /><br />Code: <a href="http://ideone.com/mGKwuL">http://ideone.com/mGKwuL</a><br /><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4698">UVa 12833 - Daily Potato</a></h2><b>Complexity: $O(26 \times N)$</b><br /><b>Category: String, Manacher's Algorithm</b><br /><br />For each query, we have to count the number of palindromes which are substring of $S$, starts and ends with given $C$ and has exactly $X$ occurrences of $C$ in it. Since it deals with palindrome, perhaps it has something to do with Manacher's Algorithm?<br /><br />With Manacher's Algorithm, we can find out the maximum length of palindrome in $O(N)$. But what's more, we can actually generate an array $M$ which gives us, for each center in the extended string $E$, ( $aba$ when extended becomes $\text{^#a#b#a#\$} $ ) the maximum length of palindrome with that particular center. How can we use this knowledge to solve our problem?<br /><br />Let us consider the character $'a'$ only. We can easily extend it for other characters by repeating the whole process $26$ times.<br /><br />Suppose we are working with center $x$. It has a palindrom from it with length $y$. Therfore, in extended string of manacher, the palindrome starts from $x-y$ and ends in $x+y$. Now, how many times does $'a'$ occurs in this palindrome? Using Cumulative Sum it is possible to answer in $O(1)$. Let that occurance be $f = \text{# of times 'a' occurs in palindrome with center x}$. Let us mark this in another array $freq$. That is we mark it like $\text{freq[f]++}$, meaning we have $freq[f]$ palindromes where $'a'$ occurs $f$ times. But wait, what if the palindrome does not start and end with $'a'$? Simple, we keep on throwing the leading and trailing character until it starts and ends with $'a'$ and it will still have $f$ occurances of $'a'$ in it.<br /><br />So we repeat this for all possible center. Now, if the query is find number of palindromes that starts and ends with $'a'$ and $'a'$ occurs exactly $X$ times, how do we solve it?<br /><br />First of all, our result will contain $res = freq[X]$ in it. What's more, our result will also contain $res \text{+=} freq[X+2] + freq[X+4] + freq[X+6] + ...$. Why is that? Take any palindrome that contains more than $X$ occurances of $'a'$. Since they start and end with $'a'$, we can just throw them out of that palindrome and reduce the occurance of $'a'$ in it by $2$. After that, we keep on trimming down head and tail of that palindrome until we reach $'a'$ again. That is, a palindrome with $Y$ occurrences of $'a'$ can be trimmed down to palindrome with $Y-2$, $Y-4$, $Y-6$, $...$ occurrences of $'a'$.<br /><br />Instead of getting the result $res = freq[X] + freq[X+2] + freq[X+4] + ... $ we can just use cumulative sum again to find it in $O(1)$ here. Just find cumulative sum of alternate terms.<br /><br />Code: <a href="http://ideone.com/brZKXL">http://ideone.com/brZKXL</a><br /><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4699">UVa 12834 - Extreme Terror</a></h2><b>Complexity: $O(N\: logN )$</b><br /><b>Category: Adhoc</b><br /><br />This was the easiest problem. Each shop gives me some money ($income$) and then for that shop I have to give Godfather some cut ($expense$). So for each shop I get $profit = income - expense$. So I calculate profit for each shop and then sort them. I can now skip $K$ shops. I will of course skip shops for which profit is negative as long as I am allowed to skip.<br /><br />Code: <a href="http://ideone.com/s3iHGz">http://ideone.com/s3iHGz</a><br /><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4700">UVa 12835 - Fitting Pipes Again</a></h2><b>Complexity: $O(N!\: N^2)$</b><br /><b>Category: Geometry, Trigonometry, Permutation, Packing Problem</b><br /><br />It's a packing problem. At first glance it seems like a tough problem but it's not. Let us first define few variables first.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-g0TvAUL2Rh0/VdXZs8GV5UI/AAAAAAAAC3k/YzCZXvguugg/s1600/octagonModified.jpg" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="320" src="http://1.bp.blogspot.com/-g0TvAUL2Rh0/VdXZs8GV5UI/AAAAAAAAC3k/YzCZXvguugg/s320/octagonModified.jpg" width="320" /></a></div><br />This is the polygon. Each polygon has two horizontal lines through it. We will call them low and top line. Each side of the polygon has length of $x$. The radius of the polygon is $r = \frac{x}{2} + y$.<br /><br />We are given height of each polygon. But first we will need to calculate the value of $x$ and $y$. We can find their values using trigonometry.<br /><br />$y^2 + y^2 = x^2$<br />$2y^2 = x^2$<br />$x = \sqrt{2y^2}$<br /><br />We also know that $h = y + x + y$. From that we can derive:<br /><br />$y + x + y = h$<br />$2y+x = h$<br />$2y + \sqrt{2y^2} = h$<br />$2y + y\sqrt{2} = h$<br />$y( 2 + sqrt{2} ) = h$<br />$y = \frac{h}{2} +\sqrt{2}$<br /><br />With the above, we can find the value of $y = \frac{h}{2} +\sqrt{2}$ and $x = \sqrt{2y^2}$ But why did we find their values?<br /><br />Okay, what happens when we put to polygon side by side? In order to solve this problem we need to be able to find the distance between the center of two polygon $A$ and $B$.<br /><br />Now if two polygons are of same height and they are placed side by side, then the difference between their centers will be $d = r_a + r_b$. What happens when two polygons of arbitrary height come beside each other?<br /><br />$3$ things can happen when two polygon $A$ and $B$ are placed side by side. $A$ is on left of $B$.<br /><br /><ol><li>Height of bottom line of $A$ is higher than height of top line of $B$. In this case, $A$ is so big that $B$ slides inside the radius of $A$.</li><li>Height of bottom line of $B$ is higher than height of top line of $A$. In this case $A$ is so small that it slides inside the radius of $B$.</li><li>Nobody can slide inside each others radius. So $d = r_a + r_b$.</li></ol><div>We need to calculate the value of $d$ for step $1$ and $2$. That can also be easily done using trigonometry. Try to draw some diagram yourself.</div><div><br /></div><div>So we used $x$ and $y$ to find the value of $d$ between two polygon. How do we use this to find the minimum width of box? </div><div><br /></div><div>Suppose we are trying to pack the polygons in some order. Which order? We don't know which order will give us minimum width so we try them all. There will be $N!$ order. </div><div><br /></div><div>So for each order, what will be the minimum of width. First lets take the first polygon, and put it inside the empty box such that it touches the left side of the box. We will calculate the center of each polygon relative to the left side of the box. The first box will have center at $r_0$.</div><div><br /></div><div>Now take the second box. First imagine there is no polygon in the box. There where will be the center of the second polygon? Ar $r_1$. Now, since there is a polygon, we will try to put it beside that one. Where will be the center now? $r_0 + d$. We will take the maximum possible center.</div><div></div><div>Repeat this for the third polygon. Empty box, beside first polygon, beside second polygon. Take the maximum again. Repeat this for all polygon.</div><div><br /></div><div>From all polygon, find the one with maximum center position. Add radius of that polygon to it's center to find the width.</div><div><br /></div><div>Take the minimum width from all permutation. </div><div><br /></div><div>Code: <a href="http://ideone.com/kphGx5">http://ideone.com/kphGx5</a></div><br /><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4701">UVa 12836 - Gain Battle Power</a></h2><div><b>Complexity: $O(N^2)$</b></div><div><b>Category: Interval DP, LIS, Knuth Optimization</b></div><div><br /></div><div>First we need to calculate the power of each deatheater. We can do that by calculating LIS from both direction and then adding them.</div><div><br /></div><div>Once we calculate the power, we run an interval dp between $1$ and $N$. The dp will have states $start$ and $end$. Inside the dp a loop between start and end will run, choosing different cut sections. We will take the one with minimum value.</div><div><br /></div><div>But this results in $O(N^3)$ solution. Using Knuth's optimization we can reduce it to $O(N^2)$. </div><div><br /></div><div>Code: <a href="http://ideone.com/XD6hZP">http://ideone.com/XD6hZP</a></div><div><br /></div><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4702">UVa 12837 - Hasmot Ali Professor</a></h2><div><b>Complexity: $O(100 \times |S|^2 )$</b></div><div><b>Category: String, Trie, Data Structure</b></div><div><br /></div><div> We will create two trie trees. The first one will contain all the queries. Take a query and concatanate them in a special way. Take the first string of the query and add a $\text{'#'}$ and then reverse the second string of the query and attach it to result. That is, if we have $abc$ and $pqr$ as query, then special string will be $\text{abc#rqp}$. Insert all special strings for each query in the first tries.</div><div><br /></div><div>Now, let us process the main string $S$. We will take each of its suffix and insert into the second trie. Now, when inserting the suffixes, each time a new node is created, we can say that we found a unique substring. </div><div><br /></div><div>Each time we find a unique substring we will process it further. Take the unique substring, and using brute force generate all possible special strings ( the one we made using query strings ) with the first and last characters of the unique string. We don't need to take more than $10$ characters from each end.</div><div><br /></div><div>For each of those special string we made from unique substring, we will query the first trie with it and find the node where it ends. We will add one to that node number in a global array.</div><div><br /></div><div>Once we finish processing all nodes in the second trie, we will simply traverse the first trie according to each query and print the result found in that node.</div><div><br /></div><div>Code: <a href="http://ideone.com/fsYMxI">http://ideone.com/fsYMxI</a></div><div><br /></div><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4703">UVa 12838 - Identity Redemption</a></h2><div><b style="background-color: red;"><span style="color: white;">Complexity: Unknown</span></b></div><div><b style="background-color: red;"><span style="color: white;">Category: Matching on General Graph</span></b></div><div><br /></div><div>I didn't manage to solve this yet. It looks like matching on general graph.</div><div><br /></div><h2><a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4704">UVa 12839 - Judge in Queue</a></h2><div><b>Complexity: $O(N\:logN)$</b></div><div><b>Category: Data Structure, Heap</b></div><div><br /></div><div>We want to minimize the waiting time for each person. So first we will sort the people in descending order. Next we will create a priority queue, which will contain all the service center along with the information about when that service center will be free. Initially all service centers are free.</div><div><br /></div><div>Now, we take each person and take out the service center that will get free at the earliest time. The person had to originally wait $x$ minutes and it took $y$ minutes for the service to get free again. So the person had to wait $x+y$ minutes. We insert the service again inside the heap and update its free time by the time it takes to serve one person.</div><div><br /></div><div>Repeat and keep track of the highest time a person had to wait.</div><div><br /></div><div>Code: <a href="http://ideone.com/KVHnTX">http://ideone.com/KVHnTX</a></div><div><br /></div><div><hr /></div><div><br /></div><div>I hope the details are clear enough for everyone to understand. Let me know if you find any mistakes.</div><div><br /></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-10354148456654042132015-08-17T17:56:00.002+06:002015-08-21T23:19:31.925+06:00Modular Exponentiation <h2>Problem</h2><div>Given three positive integers $B, P$ and $M$, find the value of $B^P \:\%\: M$.</div><div><br /></div><div>For example, $B=2$, $P=5$ and $M=7$, then $B^P \:\%\: M = 2^5 \: \% \: 7 = 32 \: \% \: 7 = 4$.</div><div><br /></div><div>This problem is known as [1]<a href="https://en.wikipedia.org/wiki/Modular_exponentiation">Modular Exponentiation</a>. In order to solve this problem, we need to have basic knowledge about [2]<a href="http://forthright48.blogspot.com/2015/07/introduction-to-modular-arithmetic.html">Modular Arithmetic</a>. </div><br /><h2>Brute Force - Linear Solution $O(P)$</h2><div>As usual, we first start off with a naive solution. We can calculate the value of $B^P \:\%\: M$ in $O(P)$ time by running a loop from $1$ to $P$ and multiplying $B$ to result.</div><div><pre class="brush:cpp;">int bigmodLinear ( int b, int p, int m ) {<br /> int res = 1 % m;<br /> b = b % m;<br /> for ( int i = 1; i <= p; i++ ) {<br /> res = ( res * b ) % m;<br /> }<br /> return res;<br />}</pre></div><div>A simple solution which will work as long as the value of $P$ is small enough. But for large values of $P$, for example, $P > 10^9$, this code will time out. Also note that in line $2$, I wrote $res = 1 \: \% \: m$. Some people tend to write $res = 1$. This will produce wrong result when the value of $P$ is $0$ and value of $M$ is $1$.<br /><br /><h2>Divide and Conquer Approach - $O(log_2(P) )$</h2></div><div>According to [3]<a href="https://en.wikipedia.org/wiki/Divide_and_conquer_algorithms">Wiki</a>,</div><blockquote class="tr_bq">A divide and conquer algorithm works by recursively breaking down a problem into two or more sub-problems of the same (or related) type (divide), until these become simple enough to be solved directly (conquer).</blockquote>That is, instead of solving the big problem $B^P \:\%\: M$, we will solve smaller problems of similar kind and merge them to get answer to the big problem.<br /><br />So how do we apply this D&C (Divide and Conquer) idea?<br /><br />Suppose we are trying to find the value of $B^P$. Then three thing can happen.<br /><ol><li><u style="font-weight: bold;">Value of P is 0 ( Base Case )</u>: $B^0$ is $1$. This will be the base case of our recursion function. Once we reach this state, we no longer divide the problem into smaller parts.</li><li><b><u>Value of P is Even</u></b>: Since $P$ is even, we can say that $B^P = B^{\frac{P}{2}} \times B^{\frac{P}{2}}$. For example, $2^{32} = 2^{16} \times 2^{16}$, $3^6 = 3^3 \times 3^3$. Therefore, instead of calculating the value of $x = B^P$, if we find the value of $y = B^{\frac{P}{2}}$, then we can get the value of $x$ as $x = y \times y$. </li><li><b><u>Value of P is Odd</u></b>: In this case we can simply say that $B^P = B \times B^{P-1}$.</li></ol><div>Using these three states, we are can formulate a D&C code.</div><div><pre class="brush:cpp;">int bigmodRecur ( int b, int p, int m ) {<br /> if ( p == 0 ) return 1%m; ///Base Case<br /><br /> if ( p % 2 == 0 ) { ///p is even<br /> int y = bigmodRecur ( b, p / 2, m );<br /> return ( y * y ) % m; ///b^p = y * y<br /> }<br /> else {<br /> ///b^p = b * b^(p-1)<br /> return ( b * bigmodRecur ( b, p - 1, m ) ) % m;<br /> }<br />}<br /></pre></div>In line $2$ we have the base case. We return $1\:\%\:m$ in case value of $m$ is $1$. Next on line $4$ we check if $p$ is even or not. If it is even, we calculate $A^{\frac{P}{2}}$ and return after squaring it. In line $8$, we handle the case when $P$ is odd.<br /><br />At each step we either divide $P$ by $2$ or subtract $1$ from $P$ and then divide it by $2$. Therefore, there can be at most $2 \times log_2(P)$ steps in D&C approach. Hence complexity is $O(log_2(P) )$.<br /><br /><h2>A Better Solution</h2><div>A better solution for this problem exists. By using Repeated Squaring algorithm, we can find the Modular Exponentiation faster than D&C. Repeated Squaring Algorithm has the same complexity as D&C, but it is iterative, thus does not suffer from recursion overhead.</div><div><br /></div><div>We need to know about bitwise operations before we learn Repeated Squaring algorithm. Since we did not cover this topic yet, I won't write about this approach now. Hopefully, once we cover bitwise operations I will write another post</div><div><br /></div><h2>Resource</h2><ol><li>Wikipedia - <a href="https://en.wikipedia.org/wiki/Modular_exponentiation">Modular Exponentiation</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/07/introduction-to-modular-arithmetic.html">Introduction to Modular Arithmetic</a></li><li>Wikipedia - <a href="https://en.wikipedia.org/wiki/Divide_and_conquer_algorithms">Divide and conquer algorithms</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/08/introduction-to-number-systems.html">Introduction to Number Systems</a></li></ol><div><br /></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-28441969145536145662015-08-15T12:05:00.000+06:002015-09-15T19:28:03.083+06:00Leading Digits of Factorial<h1>Problem</h1>Given an integer $N$, find the first $K$ leading digits of $N!$.<br /><br />For example, for $N=10$ and $K=3$, then first $3$ leading digits of $10! = 3628800$ is $362$.<br /><br />Finding leading digits uses concepts similar to [1]<a href="http://forthright48.blogspot.com/2015/08/number-of-trailing-zeroes-of-factorial.html">Number of Trailing Zeroes of Factorial</a>.<br /><br /><h2>Brute Force Solution</h2>Finding the value of $N!$ and then printing the first $K$ digits is a simple but slow solution. Using $long\:long$ we can calculate value $N!$ up to $N \leq 20$ and using Big Integer we can calculate arbitrary $N!$ but with complexity worse than $O(N^2)$.<br /><br /><h2>Solution Using Logarithm</h2><div>In [1], we say that a logarithm of value $x$ is $y$ such that $x = 10^y$. For now let us find out leading digits of a value $x$ instead of $N!$. We will extend it to cover factorials later.</div><div><br /></div><div>So, we know that $log_{10}(x) = y$, where $y$ will be some fraction. Let us separate $y$ into its integer and decimal part and call them $p,q$. For example, if $y = 123.456$, then $p=123$ and $q=0.456$.</div><div><br /></div><div>Therefore, we can say that $log_{10}(x) = p + q$. Which means, $x = 10^y = 10^{p+q}=10^p \times 10^q$.</div><div><br /></div><div>Now expand the values of $10^p$ and $10^q$. If $A=10^p$, then $A$ will simply be a power of $10$ since $p$ is an integer. To be more exact, $A$ will be $1$ with $p$ trailing zeroes. For example, $A=10^3 = 1000$. What about $B=10^q$?</div><div><br /></div><div>Since $q$ is a fraction which is $0 \leq q < 1$, value of $B$ will be between $10^0 \leq B < 10^1$, that is, $1 \leq B < 10$.</div><div><br /></div><div>Okay, we got the value of $A$ and $B$, what now? We know that if we multiply $A$ and $B$ we will get $x$. But don't multiply them just yet. Think for a bit what will happen when we multiply a decimal number with $10$. If it is integer, it will get a trailing zero, e.g, $3 \times 10 = 30$. But if it is a fraction, its decimal point will shift to right, e.g $23.65 \times 10 = 236.5$. Actually, decimal points shifts for integer numbers too, since integer numbers are real numbers with $0$ as fraction, e.g $3 = 3.00$. So in either case multiplying $10$ shifts decimal point to the right. </div><div><br /></div><div>So what happens if we multiply, $A$, which is just $10^p$ to $B$? Since $A$ has $10$ in it $p$ times, the decimal point of $B$ will shift to right $p$ times. That is all $A$ does to $B$ is change its decimal point. It does not change the digits of $B$ in any way. Thus, $B$ contains all the leading digits of $x$.</div><div><br /></div><div>For example, $log_{10}(5420) = 3.7339993 = 3 + 0.7339993$. $\therefore B = 10^0.7339993 = 5.4200$. </div><div><br /></div><div>So, if we need first $K$ leading digits of $x$, we just need to multiply $B$ with $10^{K-1}$ and then throw away the fraction part. That is $res = \lfloor B \times 10^{K-1} \rfloor$. Why $10^{K-1}$ not just $10^K$? That's because we already have $1$ leading digit present in $10^q$ before shifting it.</div><div><br /></div><h3>Extending to Factorial</h3><div>It is easy to extend the idea above to $N!$. First we need to find out the value of $y=log_{10}(N!)$.</div><div><br /></div><div>$y= log_{10}(N!)$</div><div>$y= log_{10}(N \times (N-1) \times (N-2) \times...\times 3 \times 2 \times 1 )$</div><div>$y= log_{10}(N) + log_{10}(N-1) + log_{10}(N-2) + ... + log_{10}(2) + log_{10}(1) $</div><div><br /></div><div>So we can simply find out the value of $y$ by running a loop from $1$ to $N$ and taking its log value.</div><div><br /></div><div>After that we decompose $y$ into $p$, integer part and $q$, fraction part. The answer will be $\lfloor 10^q \times 10^{K-1}\rfloor$.</div><div><br /></div><h3>Code</h3><div><pre class="brush:cpp;">const double eps = 1e-9;<br /><br />/// Find the first K digits of N!<br />int leadingDigitFact ( int n, int k ) {<br /> double fact = 0;<br /><br /> ///Find log(N!)<br /> for ( int i = 1; i <= n; i++ ) {<br /> fact += log10 ( i );<br /> }<br /><br /> ///Find the value of q<br /> double q = fact - floor ( fact+eps );<br /><br /> double B = pow ( 10, q );<br /><br /> ///Shift decimal point k-1 times<br /> for ( int i = 0; i < k - 1; i++ ) {<br /> B *= 10;<br /> }<br /><br /> ///Don't forget to floor it<br /> return floor(B+eps);<br />}<br /></pre>The code does exactly what we discussed before. But note the $eps$ that we added when flooring value in line $12$ and $22$. This due to precision error when dealing with real numbers in C++. For example, due to precision error sometimes a value which is supposed to be $1$, becomes $0.9999999999999$. The difference between these two values is very small, but if we floor them both, the first one becomes $1$ whereas the second one becomes $0$. So in order to avoid this error, when flooring a positive value we add a small number ( epsilon = $0.000000001$ ) to the number.<br /><br /></div><h2>Summary</h2><div>We need to execute the following steps to find the first $K$ leading digits of a number $x$ ( in our problem $x=N!$ ):</div><div><ol><li>Find the log value of the number whose leading digits we are seeking. $y=log_{10}(x)$.</li><li>Decompose $y$ into two parts. Integer part $p$ and fraction part $q$.</li><li>The answer is $\lfloor 10^q \times 10^{K-1} \rfloor$.</li></ol><div><br /></div></div><h2>Resource</h2><div><ol><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/08/number-of-trailing-zeroes-of-factorial.html">Number of Trailing Zeroes of Factorial</a></li></ol></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com4tag:blogger.com,1999:blog-811371107174953309.post-34400324935571789722015-08-12T10:54:00.000+06:002015-08-12T10:54:29.351+06:00Number of Trailing Zeroes of Factorial<h2>Problem</h2><div>Given a positive integer $N$, find the number of trailing zero $N!$ has.</div><div><br /></div><div>For example, $5! = 120$ has $1$ trailing zero and $10! = 3628800$ has $2$.</div><div><br /></div><h2>Brute Force Solution</h2><div>Brute Force solution for this problem involves calculating the value of $N!$, which is overkill in this situation. Read the section on Brute Force on this article before you continue - [1] <a href="http://forthright48.blogspot.com/2015/08/number-of-digits-of-factorial.html">Number of Digits of Factorial</a>. </div><div><br /></div><div>The brute force solution for this problem faces same limitation as [1]. We cannot calculate $N!$ for $N > 20$ using long long variable and using Big Integer is too slow, $O(N^2)$, if we assume number of digits in $N!$ increase by $1$ in each step.</div><div><br /></div><h2>Solution Using Factorization</h2><div>Suppose a there is a number $x$, with $y$ trailing zeroes at its end, then can't we write that number as $x = z \times 10^y$?</div><div><br /></div><div>For example, $x = 12300$ and it has $2$ trailing zeroes. Then we can write $x = 123 \times 10^2$.</div><div><br /></div><div>Therefore, for every factor of $10$, $x$ gets a trailing zero. If we can find out number of $10$ in $N!$ then we can easily count its number of trailing zeroes.<br /><br />How do we form a $10$? We form a $10$ by multiplying $2$ and $5$. So we need to find out how many $2$ and $5$ $N!$ has. This can be found using the idea for factorizing $N!$. The idea is discussed in [2]<a href="http://forthright48.blogspot.com/2015/08/prime-factorization-of-factorial.html">Prime Factorization of Factorial</a>. We will be using $\text{factorialPrimePower}()$ function from that post.<br /><br />So all we need to do is call $x = \text{factorialPrimePower}(2)$ and $y = \text{factorialPrimePower}(5)$ to find out the frequency of those primes. For every pair of $(2,5)$ we get one $10$ factor. So how many pairs can we make if we have $x$ number of $2$'s and $y$ number of $5$'s? We can make $MIN(x,y)$ pairs.<br /><br />For example, for $10!$ we have $x = \frac{10}{2} + \frac{10}{4} + \frac{10}{8}= 5 + 2 + 1 = 8$ and $y = \frac{10}{5} = 2$. Therefore number of trailing zero is $MIN(x,y) = MIN(8,2) = 2$.<br /><br /><h2>Trailing Zeroes in Different Base</h2></div><div>We solved the problem for decimal base. Now what if we want to know how many trailing zero does $N!$ has if we convert $N!$ to base $B$?</div><div><br /></div><div>For example, how many trailing zero does $10!$ has in base $16$? In order to solve this, we need to know how number system works. Read the post [3]<a href="http://forthright48.blogspot.com/2015/08/introduction-to-number-systems.html">Introduction to Number System</a> to learn more.</div><div><br /></div><div>Previously we said that for decimal numbers, multiplying them with $10$ increase their trailing zeroes. Can something similar be said for the base $16$? Yes. Look at the following values:</div><div><br /></div><div>$(1)_{16} = 1 \times 16^0$</div><div>$(10)_{16} = 1 \times 16^1$</div><div>$(100)_{16} = 1 \times 16^2$</div><div><br /></div><div>Everytime we multiply a number with its base, all its digits shift a place to left and at the end a $0$ is added. So for a number in base $B$, if we multiply it by $B$ it gets a trailing zero at its end. </div><div><br /></div><div>So all we need to do is find out how many $B$ does $N!$ has in it. For base $16$, we need to find out how many times $16 = 2^4$ occurs in $N!$. For that we find out how many times $2$ occurs using $x = \text{factorialPrimePower}(2)$ and since we get $16$ for each $2^4$, we conclude that $N!$ has $\frac{x}{4}$ trailing zeroes.</div><div><br /></div><div>This is extended for any base. Factorize the base $B$ and find out occurances of each prime. Then figure out how many combinations we can make.</div><div><br /></div><div>For example if base is $60$, then $60 = 2^2 \times 3 \times 5$. So we find out $x = \text{factorialPrimePower}(2)$, $y = \text{factorialPrimePower}(3)$ and $z = \text{factorialPrimePower}(5)$. But then, we see that we need $2^2$ in $60$, so we can't use all the $x$ we calculated, we need to pair them. So it becomes $x = \frac{x}{2}$. Now, using $x,y,z$ we can make $60$ only $MIN(x,y,z)$ times. So this will be number of trailing zero.</div><div><br /><h2>Summary</h2><div>We find number of trailing zero using the following steps:</div><div><ol><li>Factorize the base $B$</li><li>If $B = p_1^{a_1} \times p_2^{a_2}... \times p_k^{a_k}$, then find out occurance of $x_i = \text{factorialPrimePower} ( p_i)$.</li><li>But we can't use $x_i$ directly. In order to create $B$ we will need to combine each $p_i$ into $p_i^{a_i}$. So we divide each $x_i$ by $a_i$.</li><li>Number of trailing zero is $MIN(x_1,x_2,...,x_k)$.</li></ol></div><h2>Resources</h2><br /><ol><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/08/number-of-digits-of-factorial.html">Number of Digits of Factorial</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/08/prime-factorization-of-factorial.html">Prime Factorization of Factorial</a></li><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/08/introduction-to-number-systems.html">Introduction to Number System</a> </li><li>Wikipedia - <a href="https://en.wikipedia.org/wiki/Trailing_zero">Trailing Zero</a></li></ol></div><div></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-54087663803705153762015-08-10T20:26:00.000+06:002015-08-10T20:29:22.435+06:00Prime Factorization of Factorial<h2>Problem</h2>Given a positive integer $N$, find the prime factorization of $N!$.<br /><br />For example, $5! = 5 \times 4 \times 3 \times 2 \times 1 = 120 = 2^3 \times 3 \times 5$.<br /><br /><h2>Brute Force Solution</h2>A possible solution is to calculate the value of $x = N!$ and then prime factorize $x$. But calculating $N!$ is tedious. We cannot fit $N!$ where $N > 20$ in a long long variable. We will need to use Big Integer class and that would make things slow. I will soon write a blog post on Big Integer, until then know that using Big Integer it would take more than $N^2$ steps to calculate $N!$.<br /><br />Is there a better way?<br /><br /><h2>Limits on Prime</h2><div>Before we move on to the solution, let us first decide the limit on prime. In order to factorize $x = N!$, we have to generate prime numbers. But up to which value? Should we generate all primes less than $\sqrt{x}$? </div><div><br /></div><div>Even for a small value of $N$ like $100$, $x$ can be huge with over $100$ digits in it, thus, $\sqrt{x}$ will also be huge. Generating so many primes is not feasible. Using <a href="http://forthright48.blogspot.com/2015/07/sieve-of-eratosthenes-generating-primes.html">Sieve of Eratosthenes</a> we could generate primes around $10^8$, which is nowhere near $\sqrt{100!}$.</div><div><br /></div><div>Note that $N! = N \times (N-1) \times (N-2) \times ... \times 2 \times 1$. That is, $N!$ is a product of numbers less than $N$ only. Now, can there be any prime greater than $N$ that can divide $N!$? </div><div><br /></div><div>Suppose there is a number $A$ and we factorized it. It is trivial to realize that all its prime factors will be less than or equal to $A$. So in $N!$, which is the product of numbers less than $N$, if we decompose all those numbers to their prime factors, then they will reduce to primes less than or equal to $N$.<br /><br />For example, $6! = 6 \times 5 \times 4 \times 3 \times 2 \times 1 = (2 \times 3) \times 5 \times 2^2 \times 3 \times 2 = 2^4 \times 3^2 \times 5$.<br /><br />So the prime factors of $N!$ will be less than or equal to $N$. Generating primes till $\sqrt{N!}$ is not necessary. We just need to generate all primes less than or equal to $N$.<br /><br /><h2>Prime Factorization</h2></div><div>Now that we know the limit for primes, we are ready to begin factorizing the factorial. There is more than one way to achieve this. We will see three of them and discuss which one is best.</div><div><br /></div><h3>First - Linear Loop From $1$ to $N$</h3><div>We know that $N! = N \times (N-1) \times (N-2) \times ... \times 2 \times 1$. So we could simply factorize every number from $1$ to $N$ and add to a global array that tracks the frequency of primes. Using the code for $factorize()$ from <a href="http://forthright48.blogspot.com/2015/07/prime-factorization-of-integer.html">here</a>, we could write a solution like below.</div><pre class="brush: cpp; highlight: [10,16]">vector<int> prime;<br />int primeFactor[SIZE]; ///Size should be as big as N<br /><br />void factorize( int n ) {<br /> int sqrtn = sqrt ( n );<br /> for ( int i = 0; i < prime.size() && prime[i] <= sqrtn; i++ ) {<br /> if ( n % prime[i] == 0 ) {<br /> while ( n % prime[i] == 0 ) {<br /> n /= prime[i];<br /> primeFactor[ prime[i] ]++; ///Increment global primeFactor array<br /> }<br /> sqrtn = sqrt ( n );<br /> }<br /> }<br /> if ( n != 1 ) {<br /> primeFactor[n]++;<br /> }<br />}<br /><br />void factFactorize ( int n ) {<br /> for ( int i = 2; i <= n; i++ ) {<br /> factorize( i );<br /> }<br /><br /> ///Now We can print the factorization<br /> for ( int i = 0; i < prime.size(); i++ ) {<br /> printf ( "%d^%d\n", prime[i], primeFactor[ prime[i] ] );<br /> }<br />}<br /></pre>We pass the value of $N$ to $factFactorize()$ in line $20$, and it calculates the frequency of each prime in $N!$. It starts a loop from $2$ to $N$ and factorizes each of them. In line $4$ we have the $factorize()$ function modified a bit in line $10$ and $16 $to suit our need. When those lines find a prime factor, they increase the frequency of that prime in the $primeFactor$ array.<br /><br />It is simple and straight forward, but takes $O(N \times factorize() ) $ amount of time. We can do better.<br /><br /><h3>Second - Summation of Each Prime Frequency</h3><div>Instead of factorizing from $1$ to $N$, how about we just find out how many times each prime occurs in $N!$ and list them. If $p_1$ occurs $a_1$ times, $p_2$ occurs $a_2$ times $...$ $p_x$ occurs $a_x$ times, then $N! = p_1^{a_1} \times p_2^{a_2}\times ... \times p_x^{a_x}$.</div><div><br /></div><div>That sounds nice, but how do we find the frequency of prime factors in $N!$. Let us just focus on one prime factor, for example $2$, and find out how many times it occurs in $N!$. We will extend this idea to other primes.</div><div><br /></div><div>Let $N = 12$. How many times does $2$ occur in $12!$? We know that $12! = 12 \times 10 \times 9 \times ... \times 1$. How many numbers from $1$ to $12$ has $2$ as their prime factors? $\frac{12}{2} = 6$ numbers do and they are $\{ 2, 4, 6, 8, 10, 12 \}$. So we can say that at least $2^6$ is a factor of $12!$. But is there more?</div><div><br /></div><div>Yes. Notice that $4 = 2^2$, so it has an extra $2$ in it that we did not count. That means for each number that has $2^2$ in them as a factor, we need to add $1$ to our result. How many number are there which has $2^2$ as their factor? $\frac{12}{4} = 3$ numbers, which are $\{ 4, 8, 12 \}$. So we increase our frequency to $6 + 3 = 9$ and say we have at least $2^9$ in $12!$. But is that it?</div><div><br /></div><div>No. $8 = 2^3$ and for each number with $2^3$ as factor we add $1$ to result. So our result is $9 + \frac{12}{8} = 9 + 1 = 10$. </div><div><br />Do we try with $16 = 2^4$ now? No. $12!$ cannot have any number with factor $2^4$ since $\frac{12}{16} = 0$. So we conclude that $12!$ has $2^{10}$ as its factor and no more.</div><div><br /></div><div>Now, we extend this idea to other primes. What is the frequency of prime factor $3$ in $12!$? $\frac{12}{3} + \frac{12}{9} + \frac{12}{27} = 4 + 1 + 0 = 5$. We repeat this for all primes less than equal to $12$.</div><div><br />Therefore, we can say that for a given prime $p$, $N!$ will have $p^x$ as its prime factor where $x = \frac{N}{p} + \frac{N}{p^2} + \frac{N}{p^3} + ... \text{ Until it becomes 0 }$.<br /><br /></div><div>So, using this idea our code will look as the following.</div><pre class="brush: cpp;">void factFactorize ( int n ) {<br /> for ( int i = 0; i < prime.size() && prime[i] <= n; i++ ) {<br /> int p = prime[i];<br /> int freq = 0;<br /><br /> while ( n / p ) {<br /> freq += n / p;<br /> p *= prime[i];<br /> }<br /><br /> printf ( "%d^%d\n", prime[i], freq ); ///Printing prime^freq which is factor of N!<br /> }<br />}<br /></pre>This code factorizes $N!$ as long as we can generate all primes less than or equal to $N!$. The loop in line $6$ runs until $\frac{n}{p}$ becomes 0.<br /><div><br /></div><div>This code has 3 advantages over the "First" code. </div><div><ol><li>We don't have to write $factorize()$ code.</li><li>Using this code, we can find how many times a specific prime $p$ occurs in $N!$ in $O(log_p (N))$ time. In the "First" code, we will need to run $O(N)$ loop and add occurrences of $p$ in each number.</li><li>It has a better complexity for Factorization. Assuming the loop in line $6$ runs $log_2 (N)$ times, this code has a complexity of $O(N log_2 (N))$. The code runs faster than this since we only loop over primes less than $N$ and at each prime the loop runs only $O(log_p (N))$ times. The "First" code ran with $O(N \times factorize() )$ complexity, where $factorize()$ has complexity of $O(\frac{ \sqrt{N} }{ ln ( \sqrt{N} ) } + log_2(N) \: )$. </li></ol>This idea still has a small flaw. So the next one is better than this one.<br /><br /><h3>Three - Better Code than Two</h3><div>Suppose, we want to find out how many times $1009$ occurs in $9 \times 10^{18}!$. Let us modify the "Second" code to write another function that will count the result for us.<br /><pre class="brush: cpp;">long long factorialPrimePower ( long long n, long long p ) {<br /> long long freq = 0;<br /> long long cur = p;<br /><br /> while ( n / cur ) {<br /> freq += n / cur;<br /> cur *= p;<br /> }<br /><br /> return freq;<br />}<br /></pre></div>If we pass $n = 9 \times 10^{18}$ and $p = 1009$, it will return us $8928571428571439$. But this is wrong. The line $7$ in the code overflows resulting in wrong answer. Try it yourself. Print out the value of $cur$ in each step and see when it overflows.<br /><br />We could change the condition in line $5$ into something like this to solve the situation:<br /><pre class="brush: cpp;">while ( n / cur > 0 )<br /></pre>But this remedy won't work if $cur \times p$ overflows into a positive number. If we want we could use techniques that avoids multiplying two numbers if it crosses a limit, but there is an simpler way.<br /><br />Note that $\frac{N}{p^3}$ is same as $\frac{ \frac{N}{p^2} }{p}$. So instead of saying that $res = \frac{N}{p} + \frac{N}{p^2}...$, we could rewrite it as<br /><br />$x = N: \: res = res + \frac{x}{p}$.<br />$x = \frac{N}{p} = \frac{x}{p}: \: res = res + \frac{x}{p}$.<br />$x = \frac{N}{p^2} = \frac{ \frac{N}{p} } {p} = \frac{x}{p}: \: res = res + \frac{x}{p}$.<br />$...$<br />$x = 0$<br /><br />Instead of raising the power of $p$, we divide the value of $N$ by $p$ at each step. This has the same effect.<br /><br />So our code for finding frequency of specific prime should look like following:<br /><pre class="brush: cpp">long long factorialPrimePower ( long long n, long long p ) {<br /> long long freq = 0;<br /> long long x = n;<br /><br /> while ( x ) {<br /> freq += x / p;<br /> x = x / p;<br /> }<br /><br /> return freq;<br />}</pre>There might still be inputs for which this code will overflow, but chances for that is now lower.. Now if we send in $n = 9 \times 10^{18}$ and $p = 1009$, then this time we get $8928571428571425$ as our result.<br /><br />If we apply this improvement in our $factFactorize()$ function, then it will become:<br /><pre class="brush:cpp;">void factFactorize ( int n ) {<br /> for ( int i = 0; i &lt; prime.size() &amp;&amp; prime[i] &lt;= n; i++ ) {<br /> int x = n;<br /> int freq = 0;<br /><br /> while ( x / prime[i] ) {<br /> freq += x / prime[i];<br /> x = x / prime[i];<br /> }<br /><br /> printf ( "%d^%d\n", prime[i], freq );<br /> }<br />}<br /></pre>This code has less chance to overflow so it is better.<br /><br /><h2>Resource</h2><div><ol><li>forthright48 - <a href="http://forthright48.blogspot.com/2015/07/prime-factorization-of-integer.html">Prime Factorization of an Integer</a></li></ol></div></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-24756726709752040872015-08-08T15:21:00.001+06:002015-08-10T10:54:20.502+06:00Number of Digits of Factorial<h2>Problem</h2><div>Given an integer $N$, find number of digits in $N!$.</div><div><br /></div><div>For example, for $N=3$, number of digits in $N! = 3! = 3\times 2\times 1 = 6$ is $1$. For $N=5$, number of digits in $5! = 120$ is $3$.</div><div><br /></div><h2>Brute Force Solution</h2><div>The first solution that pops into mind is to calculate $N!$ and count how many digits it has. A possible solution will look like the following:</div><pre class="brush:cpp;">int factorialDigit ( int n ) {<br /> long long fact = 1;<br /> for ( int i = 2; i <= n; i++ ) {<br /> fact *= i;<br /> }<br /><br /> int res = 0; ///Number of digit of n!<br /> while ( fact ) { /// Loop until fact becomes 0<br /> res++;<br /> fact /= 10; ///Remove last digit<br /> }<br /><br /> return res;<br />}<br /></pre>This code works, but only for $N \leq 20$. Once $N$ crosses $20$, it no longer fits in a "long long" variable.<br /><br />Since factorial of $N > 20$ overflows $long\: long$, how about we use "Big Integer" to store the value?<br /><br /><h2>Brute Force Solution with Big Integer</h2><div>If you don't know what Big Integer is, then know that it is a class for arbitrary large integer. C++ does not support this class so we will have to manually implement it if we want to use C++ to solve problems involving large integers. Or you could use a programming language that supports this class, such as Java.</div><div><br /></div><div>I will write a post on Big Integer someday, but for now you could just use the Big Int class written by <a class="g-profile" href="https://plus.google.com/114184679382905647826" target="_blank">+Anudeep Nekkanti</a> from <a href="https://github.com/anudeep2011/programming/blob/master/bigint.cpp">here</a>. </div><div><br /></div><div>Multiplication of a Big Integer and an integer takes $O(\text{number of digits of Big Integer})$. Assuming that when calculating $N!$, the number of digits of $N!$ increases by $1$ at each step, it will take $O(N^2)$ time to compute $N!$. But obviously $N!$ does not increase by $1$ digit at each step ( for e.g, multiply by $100$ increases it by $2$ digits ), so worst time complexity is worse than $O(N^2)$.</div><div><br /></div><h2>Solution Using Logarithm</h2><div>Logarithm of a number is connected to its number of digits, which might not be apparent. What is logarithm? Logarithm of a number $x$, in base $b$, is a real number $y$ such that $x = b^y$. For example: </div><div>$$log_{10}1234 = 3.0913151597 \text{ and } 10^{3.0913151597} = 1234$$</div><div>In logarithms, base of the number is important. Since we want number of digits of $N!$ in decimal, we will work with base $10$.</div><div><br /></div><h3>Number of Digit of an Integer</h3><div>First, we will apply the logarithm idea on an integer. So where is the connection of logarithm with digit number? Look at the following logarithm values:</div><div><br /></div><div>$log_{10}(x) = y$<br />$log_{10}(1) = 0$</div><div>$log_{10}(10) = 1$</div><div>$log_{10}(100 )= 2$</div><div>$log_{10}(1000) = 3$</div><div>$log_{10}(10000) = 4$</div><div><br /></div><div>As the value of $x$ increases, value of $y$ also increases. Every time we multiple $x$ by $10$, value of $y$ increases by $1$. That is, every time number of digit increases, value of $y$ increases. From this table, we can infer few things about log of other values.</div><div><br /></div><div>If $log_{10}(100)$ is 2, and $log_{10}(1000)$ is 3, then for all value of $x$ where $100 < x < 1000$ , value of $y$ will lie between $2 < y < 3$. Let us try this out.</div><div><br /></div><div>$log_{10}(100) = 2$</div><div>$log_{10}(150) = 2.17609125906$</div><div>$log_{10}(500) = 2.69897000434$</div><div>$log_{10}(999) = 2.99956548823$</div><div><br /></div><div>Now note that, for every $100 \leq x < 1000$, value of $y$ is $2 \leq y < 3$. Can you see some relation between value of $y$ and number of digits of $x$?</div><div><br /></div><div>Yes. If the value of $y$ is of form $2.XXX$, then $x$ has $2+1=3$ digits. </div><div>$$\therefore \text{number of digits of}\: x = \lfloor log_{10} (x) \rfloor + 1$$</div><div><span style="background-color: yellow;">Be careful, $\lfloor log_{10} (x) \rfloor + 1$ is not same as $\lfloor log_{10} (x) + 1 \rfloor$.</span> Try it out with $100,1000,10000$. We need to floor the log value <b>before </b>we add 1. </div><pre class="brush: cpp;">int numberDigit ( int n ) {<br /> int wrongAnswer = log10(n) + 1; ///This is wrong.<br /> int rightAnswer = ( (int) log10(n) ) + 1; ///This is right.<br /> return rightAnswer;<br />}<br /></pre>In line $3$, we type cast $log10(n)$ to integer. This has same action as $floor()$ function. Also note that we used $log10()$ function instead of $log()$ function. Unlike our calculators, in C++ $log()$ has base $2$.<br /><br /><h3>Extending To Factorial</h3><div>So how do we extend this idea to $N!$?</div><div><br /></div><div>Let $x = log_{10}(N!)$. Then our answer will be $res = \lfloor x \rfloor + 1$. So all we need to do is find value of $x$.</div><div><br /></div><div>$x = log_{10}(N!)$</div><div>$x = log_{10}(1 \times 2 \times 3 \times ... \times N)$</div><div>$\therefore x = log_{10}(1) + log_{10}(2) + log_{10}(3) + ... + log_{10}(N)$ <span style="background-color: #ffd966;">This is using the law $ log_{10}(ab) = log_{10}(a)+log_{10}(b) $</span></div><br />So in order to calculate $x = log_{10}(N!)$, we don't have to calculate value of $N!$. We can simply add log value of all numbers from $1$ to $N$. This can be achieved in $O(N)$.<br /><pre class="brush: cpp">int factorialDigit ( int n ) {<br /> double x = 0;<br /> for ( int i = 1; i <= n; i++ ) {<br /> x += log10 ( i );<br /> }<br /> int res = ( (int) x ) + 1;<br /> return res;<br />}</pre><br /><h2>Digits of $N!$ in Different Base</h2><div>Now what if we want to find how many digits $N!$ has if we convert $N!$ to some other base. </div><div><br /></div><div>For example, how many digits $3!$ has in binary number system with base $2$? We know that $(6)_{10} = (110)_2$. So $3!$ has $3$ digits in base $2$ number system.</div><div><br /></div><div>Can we use logarithms to solve this problem too? Yes. </div><div>$$\text{number of digits of x in base B} = log_B(x)$$</div><div>All we need to do is change the base of our $log$ and it will find number of digits in that base.</div><div><br /></div><div>But, how do we change base in our code? We can only use log with base $2$ and $10$ in C++. Fear not, we can use the following law to change base of logartihm from $B$ to $C$.</div><div>$$log_B(x) = \frac{log_C(x)}{log_C(B)}$$</div><div>So in C++, we will use $C = 2$ or $C=10$ to find value of $log_B(x)$.</div><pre class="brush:cpp; highlight:[4];">int factorialDigitExtended ( int n, int base ) {<br /> double x = 0;<br /> for ( int i = 1; i <= n; i++ ) {<br /> x += log10 ( i ) / log10(base); ///Base Conversion<br /> }<br /> int res = ( (int) x ) + 1;<br /> return res;<br />}<br /></pre>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com4tag:blogger.com,1999:blog-811371107174953309.post-62886576510559160932015-08-04T22:05:00.000+06:002015-08-04T22:05:10.040+06:00UVa 10407 - Simple division<h2>Problem </h2><div>Problem Link - <a href="https://uva.onlinejudge.org/index.php?option=onlinejudge&page=show_problem&problem=1348">UVa 10407 - Simple division</a></div><div><br /></div><div>Given an array of numbers, find the largest number $d$ such that, when elements of the array are divided by $d$, they leave the same remainder</div><div><br /></div><h2>Solution</h2><div>We will be using our knowledge about <a href="http://forthright48.blogspot.com/2015/07/introduction-to-modular-arithmetic.html">Congruence Relation</a> and <a href="http://forthright48.blogspot.com/2015/07/euclidean-algorithm.html">GCD</a> to solve this problem. But first, we need to make sure that we understand the problem properly.</div><div><br /></div><h3>Wrong Approach</h3><div>At first, one might think that $d$ is the $gcd$ of array elements. Surely, dividing elements of the array with $d$ will leave same remainder $0$, but that doesn't necessarily mean it is the largest possible answer. </div><div><br /></div><div>For example, suppose the array is $A = \{2,8,14,26\}$. Then what will be the $gcd$ of array A? $gcd(2,8,14,26) = 2$. Dividing elements of $A$ with $2$ leaves us $0$. So is this our answer? No.</div><div><br /></div><div>For array $A$, there exists a number greater than $2$ that leaves the same remainder. And that number is $3$. When we divide each element with $3$, it leaves us $2$ as the remainder.</div><div><br />The problem did not ask us to find the number that will leave $0$ as remainder. It asked us to find the largest number $d$ that leaves the same remainder. So for array $A$ the $gcd(A)$ is not the answer. Is $3$ our answer? No. The answer is given below in Example.</div><div><br /></div><h3>Using Congruence Relation </h3><div>Let us rephrase the problem using congruence relation. Suppose there is an array with elements, $A=\{a,b,c\}$. We need to find the largest number $d$ that leaves the same remainder for each of its element.</div><div><br /></div><div>Let us consider only $\{a,b\}$ for now. We need to find $d$ such that it leaves the same remainder when it divides $a$ and $b$. Meaning</div><div>$$a \: \equiv \: b \: \text{(mod d)} \\ a - b \: \equiv \: 0 \: \text{(mod d)}$$</div><div>What does this mean? It means, if $d$ leaves the same remainder when it divides $a$ and $b$, then it will leave no remainder when dividing $a-b$.</div><div><br />Using same logic, we can say that $b - c \: \equiv \: 0 \: \text{(mod d)}$. Therefore, if we find a number that divides $a-b$ and $b-c$, then that number will leave the same remainder when dividing $\{a,b,c\}$. </div><div><br /></div><div>This idea can be extended for any number of elements in the array. So $d$ is such a number that leaves no remainder when it divides the difference of adjacent elements in the array.</div><div><br /></div><div>But wait, there are multiple numbers that can divide the difference of adjacent elements. Which one should we take? Since we want the largest value of $d$, we will take the largest divisor that can divide all the difference of adjacent numbers. There is only one divisor that can divide all the adjacent differences, and that is $gcd( (a-b), (b-c) )$.<br /><br />Therefore, if we have an array $A=\{a_1,a_2,a_3...a_n\}$, then $d = gcd( a_2 - a_1, a_3- a_2,...a_n - a_{n_1})$.<br /><br />Careful about negative value of $gcd()$. Make sure you take the absolute value of $gcd()$.</div><div><br /></div><h2>Example</h2><div>Let us get back to the example we were looking into before. Finding $d$ for the array $A = \{2,8,14,26\}$. We know that $d$ will divide the difference of adjacent elements completely. So let us make another array that will contain the difference of adjacent elements., $B = \{8-2, 14-8, 26-14\} = \{6,6,12\}$. Therefore, $d = gcd ( 6,6,12 ) = 6$.</div><div><br /></div><h2>Summary</h2><div><ol><li>Let the array of elements be $A = \{a,b,c,d...\}$.</li><li>$res \neq gcd(a,b,c,d...)$.</li><li>$res = gcd ( a-b,b-c,c-d,...)$.</li></ol></div><div><br /></div><h2>Code</h2><pre class="brush: cpp;">#include <bits/stdc++.h><br />using namespace std;<br /><br />typedef long long vlong;<br /><br />vlong gcd ( vlong a, vlong b ) {<br /> while ( b ) {<br /> a = a % b;<br /> swap ( a, b );<br /> }<br /> return a;<br />}<br /><br />vlong arr[1010];<br /><br />int main () {<br /><br /> while ( scanf ( "%d", &arr[0] ) != EOF ) {<br /> if ( arr[0] == 0 ) break; //End of test case<br /><br /> //A new test case has started<br /> int cur = 1;<br /><br /> //Take input<br /> while ( 1 ) {<br /> scanf ( "%lld", &arr[cur] );<br /> if ( arr[cur] == 0 ) break;<br /> else cur++;<br /> }<br /><br /> vlong g = 0; //Start with 0 since gcd(0,x) = x.<br /> for ( int i = 1; i < cur; i++ ) {<br /> int dif = arr[i] - arr[i-1]; //Calculate difference<br /> g = gcd ( g, dif ); //Find gcd() of differences<br /> }<br /><br /> if ( g < 0 ) g *= -1; //In case gcd() comes out negative<br /> printf ( "%lld\n", g );<br /> }<br /><br /> return 0;<br />}<br /></pre>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-59863520795823186642015-08-03T16:00:00.000+06:002015-08-03T16:00:08.177+06:00UVa 11388 - GCD LCM<h2>Problem</h2>Problem Link - <a href="https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&category=24&page=show_problem&problem=2383">UVa 11388 - GCD LCM</a><br /><br />Given two positive integers $(G, L)$, we have to find a pair of integers $(a,b)$ such that $gcd(a, b)=G$ and $\text{lcm}(a, b)=L$. If there are multiple such pairs, we have to find the pair where $a$ is minimum. Also, both $a$ and $b$ needs to be positive.<br /><br /><h2>Solution</h2><div>We don't know the value of $(a,b)$ yet. Here is how I approached the problem.</div><div><br /></div><h3>Value of $a$</h3><div>What we know that $gcd(a,b) = G$. So, $G$ divides $a$ and $b$. Therefore, $a$ is a multiple of $G$. We also need to make sure that $a$ is as small as possible. So, what is the smallest positive number that can be divided by $G$? The answer is $G$ itself.</div><div>$$\therefore a = G$$</div><div><br /></div><h3>Existence of Solution</h3><div>Next, we know that $lcm(a,b)$ is the smallest positive number which is divisible by both $a$ and $b$. Since $a=G$ and $a \: | \: lcm(a,b)$, it follows that $a = G$ should divide $lcm(a,b) = L$. If $L$ is not divisible by $G$, then no solution exists. </div><div>$$\therefore \text{if } (G \not| L), \text{no solution exists}$$</div><div><br /></div><h3>Value of $b$</h3><div>The value of $b$ can be derived in the following way. We know that:</div><div><br /></div><div>$gcd(a,b) \times lcm(a,b) = a \times b$</div><div>$G \times L = G \times b$</div><div>$b = \frac{G\times L}{G}$</div><div>$\therefore b = L$.</div><div><br /></div><h3>Summary </h3><div><ol><li>$a = G$</li><li>If $G \not| L$, no solution exists</li><li>$b = L$</li><li>$\therefore (a,b) = (G,L)$</li></ol><div><br /></div><h2>Code</h2></div><div>Here is the code in C++</div><pre class="brush: cpp;">#include <bits/stdc++.h><br /><br />int main () {<br /> int kase;<br /> scanf ( "%d", &kase ); //Input number of case<br /><br /> while ( kase-- ) {<br /> int G, L;<br /> scanf ( "%d %d", &G, &L ); //Take input<br /><br /> int a, b; //Need to find their value<br /> <br /> a = G;<br /> <br /> if ( L % G != 0 ) {<br /> printf ( "-1\n" ); //No Solution<br /> continue;<br /> }<br /> <br /> b = L;<br /> <br /> printf ( "%d %d\n", a, b );<br /> }<br /><br /> return 0;<br />}<br /></pre>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-36924119279693953912015-08-01T21:22:00.000+06:002015-08-01T21:32:19.311+06:00Introduction to Number SystemsNumber System is a consistent way of expressing a number. It consists of three parts:<br /><ol><li>The Symbol - It can be digits, letters or any other symbol.</li><li>The Position - The position of symbols determine their weight.</li><li>The Base - The total number of different symbols supported.</li></ol><div>Number system is often named after the size of its base. For example, the "Decimal Number System", with symbols $(0,1,2,3,4,5,6,7,8,9)$. Since there are $10$ different symbols, the base is $10$.</div><div><br />The Decimal Number System is the standard number system that we use in our everyday life. We will use this system as our reference for other systems. Since we will be looking into lots of different number systems in this article, from now on I will write numbers in decimal number system as $(number)_{10}$. In general, $(x)_{y}$ means $x$ is a number with base $y$.</div><div><br /><h2>How Number System Works</h2></div><div>As mentioned before, Number System consists of three parts.<br /><br /></div><h3>Symbols in Number System</h3><div>In decimal number system, the first few numbers we encounter are $0, 1, 2, 3, 4, 5, 6, 7, 8, 9$. Then we run out of symbols. We add a new position at the front of the number and start placing number in those to continue counting: $8, 9, 10, 11, 12, 13...$</div><div><br /></div><div>The same concept is used with number system of different base. For example, let us consider a number system with 3 different symbols: $(0,1,2)$. The first few numbers in this system are: $0,1,2$. Once we run out of symbols we introduce new positions and get: $10,11,12,20,21,22,100,101,102,110$.<br /><br /><h3>Position of Symbol</h3>In a number, the position of symbols is $0$ indexed and ranked from right to left. For example, the number $54711$ will have the following index for its positions:<br /><div class="postTable"><table><tbody><tr><td>Index</td><td>$4$</td><td>$3$</td><td>$2$</td><td>$1$</td><td>$0$</td></tr><tr><td>Symbols</td><td>$5$</td><td>$4$</td><td>$7$</td><td>$1$</td><td>$1$</td></tr></tbody></table></div></div>The last digit of the number has position index $0$. The index increases as we go left.<br /><br />The position in a number represents the weight of the symbol in that position. A symbol in a higher position values more than one in a lower position. For example in the decimal number system, we know that $2 > 1$. But if we consider $1$ in a higher position then $100 > 20$. How is the weight of each position calculated exactly? The answer is in next section.<br /><br /><h3>Base of Number System</h3><div>The "Base" of number system determine the weight of each position.<br /><br />For example, $(1234)_{10} = 1000 + 200 + 30 + 4 = 1 \times 10^3 + 2 \times 10^2 + 3 \times 10^1 + 4 \times 10^0$. </div><div><br /></div><div>Here the weight of the symbol in $0_{th}$ position is $10^0$, $1_{st}$ position is $10^1$, $2_{nd}$ position is $10^2$ and $3_{rd}$ position is $10^3$. As the position increases, the weight of position increases by $10$. A symbol in $n_{th}$ position carries $10^n$ weight.</div><div><br /></div><div>But this is true only for the decimal number system. The weight is a power of $10$ for decimal system cause it has a base of $10$. A number system with base $B$ will have weight $B^n$ for $n_{th}$ position.</div><div><br /><h2>Binary Number System</h2><div>Binary Number System is a Number System with Base $2$. It has only two symbols: $(0,1)$. This number system is used by computers to store information. As programmers, we will frequently work with binary numbers, so we need to understand this system properly.<br /><br />Since we are used to the Decimal Number System, it might be awkward to work with the binary system at first. But we will get used to it pretty soon.</div><div><br /></div><div>Here are the first few numbers in Binary System.</div><div class="postTable"><table><thead></thead><tbody><tr><td>Decimal</td><td>$0$</td><td>$1$</td><td>$2$</td><td>$3$</td><td>$4$</td><td>$5$</td><td>$6$</td><td>$7$</td><td>$8$</td><td>$9$</td></tr><tr><td>Binary</td><td>$0$</td><td>$1$</td><td>$10$</td><td>$11$</td><td>$100$</td><td>$101$</td><td>$110$</td><td>$111$</td><td>$1000$</td><td>$1001$</td></tr></tbody></table></div><br /><h3>Binary to Decimal</h3>We will first look into how we can convert a number from binary to the decimal system.<br /><br />Suppose we want to find the value of $(11001)_2$ in decimal. How do we find it? Remember that we learned that a symbol in $n_{th}$ position has $Base^n$ weight. What is the base in the binary system? $2$. So this number can be written as:<br /><br />$(11001)_2 = 1\times 2^4 + 1 \times 2^3 + 0 \times 2^2 + 0\times 2^1 + 1 \times 2^0 = 16 + 8 + 0 + 0 + 1 = (25)_{10}$<br /><br /><h3>Decimal to Binary</h3></div><div>Solve $(19)_{10} = (?)_2$. We need to convert $19$ to binary.</div><div><br /></div><div>A number $x$ in binary system is of form: $x = s_n2^n +...+s_32^3+s_22^2+s_12^1+s_02^0$, where $s_i$ is the symbol in $i_{th}$ position. What happens when we find $x \: \% \: 2$? Since all position except for $0_{th}$ position is multiple of $2$, $x \: \% \: 2 = s_0$.</div><div><br /></div><div>That is, finding the value of $x \: \% \: 2$ gives us the last digit of $x$ in the binary system. </div><div><br /></div><div>Okay, we got the digit in $0_{th}$ position, how do we find the rest? In order to find the digit in next position, we divide $x$ by $2$. What happens if we divide $x$ by $2$? </div><div><br /></div><div>$\frac{x}{2} = s_n2^{n-1} +...+s_32^2+s_22^1+s_12^0+\frac{s_0}{2}$</div><div>$\therefore \lfloor \frac{x}{2} \rfloor = s_n2^{n-1} +...+s_32^2+s_22^1+s_12^0$<br /><br />That is, when we divide $x$ by $2$, the last digit of $x$ disappears and all digits of $x$ shifts one place to right. For example, if $x$ is $(11011)_2$ then $\frac{x}{2}$ is $(1101)_2$.</div><div><br /></div><div>In C++, we can easily perform this floor division by simply performing integer division. Once we divide it, we then find the value of $\lfloor \frac{x}{2} \rfloor \: \% \: 2$. By repeating this until $x$ becomes $0$, we can find all digits.</div><div><br /></div><div>Let us try to solve the problem $(19)_{10} = (?)_2$.</div><div><br /></div><div>$x=19: x \: \% \: 2 = 1$</div><div>$x = \lfloor \frac{19}{2} \rfloor = 9: x \: \% \: 2 = 1$</div><div>$x = \lfloor \frac{9}{2} \rfloor = 4: x \: \% \: 2 = 0$</div><div>$x = \lfloor \frac{4}{2} \rfloor = 2: x \: \% \: 2 = 0$</div><div>$x = \lfloor \frac{2}{2} \rfloor = 1: x \: \% \: 2 = 1$</div><div>$x = \lfloor \frac{1}{2} \rfloor = 0: \text{end}$</div><div>$\therefore (19)_{10} = (10011)_2$<br /><br /><div><h2>Hexadecimal Number System</h2><div>Hexadecimal Number System is a Number System with Base $16$. It has 16 symbols, $(0,1,2,3,4,5,6,7,8,9,A,B,C,D,E,F)$. We already saw Decimal to Binary conversion and vice versa. Looking into Decimal-Hexadecimal conversion and vice verse will help us understand the base conversions better.</div><div><br /></div><div>Here are the first few numbers in Hexadecimal System.<br /><br /></div><div class="postTable"><table><tbody><tr><td>$Decimal$</td><td>$0$</td><td>$1$</td><td>$2$</td><td>$3$</td><td>$4$</td><td>$5$</td><td>$6$</td><td>$7$</td><td>$8$</td><td>$9$</td></tr><tr><td>$Hexadecimal$</td><td>$0$</td><td>$1$</td><td>$2$</td><td>$3$</td><td>$4$</td><td>$5$</td><td>$6$</td><td>$7$</td><td>$8$</td><td>$9$</td></tr><tr><td>$Decimal$</td><td>$10$</td><td>$11$</td><td>$12$</td><td>$13$</td><td>$14$</td><td>$15$</td><td>$16$</td><td>$17$</td><td>$18$</td><td>$19$</td></tr><tr><td>$Hexadecimal$</td><td>$A$</td><td>$B$</td><td>$C$</td><td>$D$</td><td>$E$</td><td>$F$</td><td>$10$</td><td>$11$</td><td>$12$</td><td>$13$</td></tr></tbody></table></div><br />In hexadecimal system, $(A,B,C,D,E,F)$ corresponds to the decimal value $(10,11,12,13,14,15)$.<br /><br /><h3>Hexadecimal to Decimal</h3>Suppose we want to find the value of $(1A3F2B)_{16}$ in decimal. Using the same logic as Binary to Decimal conversion, we can say that:<br /><br />$(1A3F2B)_{16} = 1\times 16^5 + A \times 16^4 + 3 \times 16^3 + F\times 16^2 + 2 \times 16^1 + B \times 16^0$<br />$ (1A3F2B)_{16} = 1\times 16^5 + (10) \times 16^4 + 3 \times 16^3 + (15)\times 16^2 + 2 \times 16^1 + (11) \times 16^0$<br />$ (1A3F2B)_{16} = 1048576 + 10 \times 65536 + 3 \times 4096 + 15 \times 256 + 2 \times 16 + 11 \times 1$<br />$ (1A3F2B)_{16} = 1048576 + 655360 + 12288 + 3840 + 32 + 11$<br />$ (1A3F2B)_{16} = (1720107)_{10}$<br /><br /><h3>Decimal to Hexadecimal</h3></div><div>Solve $(123456789)_{10} = (?)_{16}$.<br /><br /></div><div>A number $x$ in hexadecimal system is of form: $x = s_n16^n +...+s_316^3+s_216^2+s_116^1+s_016^0$, where $s_i$ is the symbol in $i_{th}$ position. </div><div><br /></div><div>Just like Binary System, finding the value of $x \: \% \: 16$ gives us the last digit of $x$ in hexadecimal system. To find the rest, we divide it by $16$ and find the last digit again.</div><div><br /></div><div>$\frac{x}{16} = s_n16^{n-1} +...+s_316^2+s_216^1+s_116^0+\frac{s_0}{16}$</div><div>$\therefore \lfloor \frac{x}{16} \rfloor = s_n16^{n-1} +...+s_316^2+s_216^1+s_116^0$</div><div><br /></div><div>Let us try to solve the problem $(123456789)_{10} = (?)_{16}$.</div><div><br /></div><div>$x=123456789: x \: \% \: 16 = 5$</div><div>$x = \lfloor \frac{123456789}{16} \rfloor = 7716049: x \: \% \: 16 = 1$</div><div>$x = \lfloor \frac{7716049}{16} \rfloor = 482253: x \: \% \: 16 = 13$</div><div>$x = \lfloor \frac{482253}{16} \rfloor = 30140: x \: \% \: 16 = 12$</div><div>$x = \lfloor \frac{30140}{16} \rfloor = 1883: x \: \% \: 16 = 11$</div><div>$x = \lfloor \frac{1883}{16} \rfloor = 117: x \: \% \: 16 = 5 $<br />$x = \lfloor \frac{117}{16} \rfloor = 7: x \: \% \: 16 = 7 $<br />$x = \lfloor \frac{7}{16} \rfloor = 0: \text{end} $</div><div>$\therefore (123456789)_{10} = (75BCD15)_{16}$<br /><br />Notice how we replaced the remainder value $11,12,13$ with $B,C,D$. Since the value $(10,11,12,13,14,15)$ represents $(A,B,C,D,E,F)$ in hexadecimal, we use the proper symbols in the number.<br /><br /><h2>Number System with Base $B$</h2><div>We saw how to work with Binary and Hexadecimal number system. We will now generalize the approach.</div><div><br /></div><div>Suppose there is number $x$ in base$B$, then $x = d_nd_{n-1}..d_3d_2d_1d_0$, where $0 \leq d_i < B$. </div><div><br /></div><h3>Base to Decimal</h3><div>The number $x$ in decimal will be $d_n \times B^n + d_{n-1} \times B^{n-1} + ... + d_1 \times B^1 + d_0 \times B^0$.</div><div><br /></div><div>Here is a code that can convert a number in base $B$ into decimal.</div><pre class="brush:cpp;">int baseToDecimal ( string x, int base ) {<br /> int res = 0;<br /> int len = x.length();<br /><br /> int coef = 1; ///initially base^0<br /> for ( int i = len - 1; i >= 0; i-- ) { ///Start from reverse<br /> res += (x[i]-'0') * coef;<br /> coef *= base; //increase power of base<br /> }<br /> return res;<br />}<br /></pre>The number $x$ here is a string. That is because $x$ can have symbols that are not digits, for example, $(AB12)_{16}$. We find out the length of the number in line $3$. Then we iterate over the number from last digit to the first digit. We start from back where the coefficient of the digit is $base^0 = 1$. Every time we change position, we multiply the coefficient by $base$. This ensures we multiply $d_n$ with $base^n$.<br /><br />If we want we can convert the number to decimal from the front of the number, that is from first digit to last digit.<br /><br />Suppose we want to make the number $d_1d_2d_3$ in base $B$. We can start with the digit $d_1$ only. Then we can move it to left of its position by multiply it by $B$, $d_1B^0 \times B = d_1B^1$. Then if we add $d_2$ with it, it becomes $d_1B^1 + d_2 = d_1d_2$. We multiply it by $B$ again and then add $d_3$. It finally becomes $d_1d_2d_3$.<br /><br />For example, $(123)_{10} = (1 \times 10 + 2 ) \times 10 + 3$.<br /><br />Using this idea, we can write the following code:<br /><pre class="brush:cpp;">int baseToDecimalAlternate ( string x, int base ) {<br /> int res = 0;<br /> int len = x.length();<br /><br /> for ( int i = 0; i < len; i++ ) {<br /> res = ( res * base ) + (x[i]-'0');<br /> }<br /> return res;<br />}<br /></pre>Here we iterate the number $x$ from $0$ to $len-1$ in line $5$. Line $6$ is what I explained above.<br /><br />The alternate method is shorter and sometimes makes a problem easier to solve. Knowing both method helps us approach problems from different directions.<br /><br /><h3>Decimal To Base</h3><div>$x \: \% \: B$ will give us the last digit of $x$ in base $B$. We can get the remaining digits by repeatedly dividing $x$ by $B$ and taking its last digit until $x$ becomes 0.</div><div><br /></div><div>$x = x: d_0 = x \: \% \: B$</div><div>$x = \lfloor \frac{x}{B} \rfloor: d_1 = x \: \% \: B$</div><div>$...$</div><div><br /></div><div>Here is how we will do it in code:</div><pre class="brush:cpp;">//A list of symbol. Depending on base and number system, this list can be different.<br />char symbol[] = {'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'};<br />string decimalToBase ( int x, int base ) {<br /> string res = "";<br /><br /> while ( x ) {<br /> int r = x % base; //Find the last digit<br /> res = res + symbol[r];//Change the integer value to symbol and append to res<br /> x /= base;//Remove the last digit<br /> }<br /> if ( res == "" ) res = symbol[0];//If res is empty, that means x is 0.<br /> reverse ( res.begin(), res.end());//We found the digits in reverse order.<br /> return res;<br />}<br /></pre>In line $2$ we have a symbol list. Since we are converting to a different base, it will have symbols for each value from $0$ to $B-1$. We keep on finding the last digit of $x$ until it becomes $0$ in the loop at line $6$. The last digit is found in line $7$. We append it to the result in line $8$ and we divide $x$ by $B$ in line $9$. In line $11$ we check if $res$ is empty. If it is then $x$ was $0$ from the beginning, so we append $0$ to $res$. In line $12$ we reverse the whole string $res$ since we generated the digits in reverse order, i.e, from last digit to the first digit.<br /><br /><h2>Resource</h2><div><ol><li>tibasicdev.wikidot - <a href="http://tibasicdev.wikidot.com/binandhex">Binary, Hexadecimal and Octal number system</a></li></ol></div></div></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-84307557011549372932015-07-29T18:22:00.000+06:002015-07-29T18:22:00.890+06:00Simple Hyperbolic Diophantine Equation<h2>Problem</h2><div>Given the integers $A, B, C, D$, find a pair of $(x,y)$ such that it satisfies the equation $Axy + Bx + Cy = D$.<br /><br />For example, $2xy + 5x + 56y = -7$ has a solution $(x,y) = (-27,64)$.</div><div><br /></div><div>Hyperbolic equations are usually of form $Ax^2 + By^2 + Cxy + Dx + Ey = F$. Here we don't have coefficients of $x^2$ and $y^2$. That's why we are calling it "Simple".</div><div><br /></div><div>So how do we tackle the problem? First we will rewrite the equation to make it easier to approach.</div><div><br /></div><h2>Rewriting The Equation</h2><div>$Axy + Bx + Cy = D$</div><div>$A^2xy + ABx + ACy = AD$, <span style="background-color: #ffd966;">Multiply A with both sides</span></div><div>$A^2xy + ABx + ACy + BC = AD + BC$, <span style="background-color: #ffd966;">Add BC to both sides</span></div><div>$Ax ( Ay + B ) + C ( Ay + B ) = AD + BC$,</div><div>$\therefore (Ax + C) (Ay + B) = AD + BC$</div><div><br /></div><div>Now that we have rewritten the equation, we can see that the problem has two cases. We will handle each case separately.</div><div><br /><h2>When $AD + BC = 0$</h2></div><div>When $AD + BC = 0$, our equation becomes $(Ax + C)(Ay + B) = 0$. So we need to have either $(Ax + C) = 0$ or $(Ay+B) = 0$ to satisfy the equation.<br /><br />So we have two possible formation of solutions. When $(Ax+C) = 0$, we have $x = \frac{-C}{A}$ and $y$ any integer value. When $(Ay+B)=0$, we have $y = \frac{-B}{A}$ and $x$ any integer value.<br /><br />So, $(x,y) = (\frac{-C}{A},k)$ or $(k,\frac{-B}{A})$, where $k$ is any integer. What if it is not possible to divide $-C$ or $-B$ with $A$? Then there is no solution of that formation.<br /><br /><h2>When $AD + BC \neq 0$</h2>Let $P = AD+BC$. Since $(Ax + C)(Ay+B) = P$, $(Ax+C)$ and $(Ay+B)$ must be divisors of $P$. Let $d_1, d_2, d_3 ... d_n$ be divisors of $P$. Then for each divisor $d_i$, if $(Ax + C) = d_i$ then $(Ay + B) = \frac{P}{d_i}$, so that we get $P$ when we multiply them.<br /><br />$(Ax + C) = d_i$<br />$Ax = d_i - C$<br />$\therefore x = \frac{d_i - C}{A}$<br /><br />$(Ay+B) = \frac{P}{d_i}$<br />$Ay = \frac{P}{d_i} - B$<br />$Ay = \frac{P - Bd_i}{d_i}$<br />$\therefore y = \frac{P - Bd_i}{Ad_i}$<br /><br />Therefore for each divisor $d_i$, we have $(x,y) = ( \frac{d_i - C}{A}, \frac{P - Bd_i}{Ad_i} )$. This solution is valid only if we get integer values for $(x,y)$.<br /><br />Careful when calculating divisors of $P$. In this problem, <span style="background-color: yellow;">we have to include the negative divisors of $P$ in our calculation</span>. For example, $4$ will have the following divisors $(1,2,4,-1,-2,-4)$.<br /><br /><h2>Example</h2>Find solutions for $2xy + 2x + 2y = 4$.<br /><br />First we will find $P = AD + BC = 2 \times 4 + 2 \times 2 = 12$. $12$ has the following divisors: $(\pm 1, \pm 2, \pm 3, \pm 4, \pm 6, \pm 12 )$.<br /><br />Since $(2x + 2)(2y + 2 ) = 12$, we get the following:<br /><br />$d_1 = 1: x = \frac{1-2}{2} = \frac{-1}{2}, y = \frac{12 - 2 \times 1}{2 \times 1} = \frac{10}{2} = 5$<br />$d_2 = -1: x = \frac{-1-2}{2} = \frac{-3}{2}, y = \frac{12 - 2 \times -1}{2 \times -1} = \frac{14}{-2} = -7$<br /><span style="background-color: #ffd966;">$d_3 = 2: x = \frac{2-2}{2} = 0, y = \frac{12 - 2 \times 2}{2 \times 2} = \frac{8}{4} = 2$</span><br /><span style="background-color: #ffd966;">$d_4 = -2: x = \frac{-2-2}{2} = -2, y = \frac{12 - 2 \times -2}{2 \times -2} = \frac{16}{-4} = -4$</span><br />$d_5 = 3: x = \frac{3-2}{2} = \frac{1}{2}, y = \frac{12 - 2 \times 3}{2 \times 3} = \frac{6}{6} = 1$<br />$d_6 = -3: x = \frac{-3-2}{2} = \frac{-5}{2}, y = \frac{12 - 2 \times -3}{2 \times -3} = \frac{18}{-1} = -18$<br />$d_7 = 4: x = \frac{4-2}{2} = 1, y = \frac{12 - 2 \times 4}{2 \times 4} = \frac{4}{8} = \frac{1}{2}$<br />$d_8 = -4: x = \frac{-4-2}{2} = -3, y = \frac{12 - 2 \times -4}{2 \times -4} = \frac{20}{-8} = \frac{-5}{2}$<br /><span style="background-color: #ffd966;">$d_9 = 6: x = \frac{6-2}{2} = 2, y = \frac{12 - 2 \times 6}{2 \times 6} = \frac{0}{12} = 0$</span><br /><span style="background-color: #ffd966;">$d_{10} = -6: x = \frac{-6-2}{2} = -4, y = \frac{12 - 2 \times -6}{2 \times -6} = \frac{24}{-12} = -2$</span><br />$d_{11} = 12: x = \frac{12-2}{2} = 5, y = \frac{12 - 2 \times 12 }{2 \times 12 } = \frac{-12}{24} = \frac{-1}{2}$<br />$d_{12} = -12: x = \frac{-12-2}{2} = -7, y = \frac{12 - 2 \times -12 }{2 \times -12 } = \frac{36}{-24} = \frac{-3}{2}$<br /><br />So the solutions are $(0,2),(-2,-4),(2,0),(-4,-2)$.<br /><br /><h2>Code</h2><div>Let us try to write a function that will count the number of solutions for the given problem above. If there are infinite solutions, the function will return $-1$. </div><pre class="brush:cpp;">bool isValidSolution ( int a, int b, int c, int p, int div ) {<br /> if ( ( ( div - c )% a ) != 0 ) return false; //x = (div - c) / a<br /> if ( ( (p-b*div) % (a*div) ) != 0 ) return false;// y = (p-b*div) /(a*div)<br /> return true;<br />}<br /><br />int hyperbolicDiophantine ( int a, int b, int c, int d ) {<br /> int p = a * d + b * c;<br /><br /> if ( p == 0 ) { //ad + bc = 0<br /> if ( -c % a == 0 ) return -1; //Infinite solutions (-c/a, k)<br /> else if ( -b % a == 0 ) return -1; //Infinite solutions (k, -b/a)<br /> else return 0; //No solution<br /> }<br /> else {<br /> int res = 0;<br /><br /> //For each divisor of p<br /> int sqrtn = sqrt ( p ), div;<br /> for ( int i = 1; i <= sqrtn; i++ ) {<br /> if ( p % i == 0 ) { //i is a divisor<br /><br /> //Check if divisors i,-i,p/i,-p/i produces valid solutions<br /> if ( isValidSolution(a,b,c,p,i) )res++;<br /> if ( isValidSolution(a,b,c,p,-i) )res++;<br /> if ( p / i != i ) { //Check whether p/i is different divisor than i<br /> if ( isValidSolution(a,b,c,p,p/i) )res++;<br /> if ( isValidSolution(a,b,c,p,-p/i) )res++;<br /> }<br /> }<br /> }<br /><br /> return res;<br /> }<br />}<br /></pre>We have two functions here. The first function $\text{isValidSolution}()$ takes input $a,b,c, p, div$ and checks whether we can get a valid solution $(x,y)$ using $div$ as a divisor of $P$.<br /><br />In line $7$ we start our function to find number of solutions. In line $10$ We check our first case when $AD + BC = 0$. If it is possible to form a valid solution, we return $-1$ in line $11,12$ otherwise $0$ in line $13$.<br /><br />From line $15$ we start case when $AD + BC \neq 0$. For this case we have to find all divisors of $P$. We know that all divisors come in pairs and one part of that pair lies before $\sqrt{P}$. We we loop till $\sqrt{P}$ to find first part of the pair and process it in line $16$. The second part can be found by dividing $P$ by the first part. We process both these parts ( if second part is different than first part, line $26$ ) and their negative parts by passing them to $\text{isValidSolution}()$ functions. If they are valid we increase our result by 1.<br /><br />If we are asked to print the solutions, we could do it by modifying the functions above. As long as we understand how the solution works, there should not be any problem in modifying the functions.<br /><br /><h2>Resource</h2><ol><li>alpertron - <a href="http://www.alpertron.com.ar/METHODS.HTM">Methods to Solve $Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$</a></li></ol></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-50797495822191285372015-07-27T21:05:00.000+06:002015-07-29T13:35:51.021+06:00Linear Diophantine Equation<h2>Problem</h2><div>Given the value of integers $A, B$ and $C$ find a pair of integers $(x,y)$ such that it satisfies the equation $Ax + By = C$.</div><br /><div>For example, if $A = 2, B = 3$ and $C = 7$, then possible solution of $(x,y)$ for equation $2x + 3y = 7$ would be $(2,1)$ or $(5,-1)$.</div><br />The problem above is a type of <a href="https://en.wikipedia.org/wiki/Diophantine_equation">Diophantine problem</a>. In the Diophantine problem, only integer solutions to an equation are required. Since $Ax + By = C$ is a linear equation, this problem is a <b>Linear Diophantine Problem</b> where we have to find a solution for a <b>Linear Diophantine Equation</b>.<br /><br />For now, let us assume that $A$ and $B$ are non-zero integers.<br /><h2>Existence of Solution</h2>Before we jump in to find the solution for the equation, we need to determine whether it even has a solution. For example, is there any solution for $2x + 2y = 3$? On the left side we have $2x + 2y$ which is even no matter what integer value of $(x,y)$ is used and on the right side we have $3$ which is odd. This equation is impossible to satisfy using integer values.<br /><br />So how do we determine if the equation has a solution? Suppose $g = gcd(A,B)$. Then $Ax + By$ is a multiple of $g$. In order to have a valid solution, since left side of the equation is divisible by $g$, the right side too must be divisible by $g$. Therefore, if $g \not| \: C$, then there is no solution.<br /><br /><h2>Simplifying the Equation</h2>Since both side of equation is divisible by $g$, i.e, $g \: | \: \{ \: (Ax + By), C \: \}$, we can safely divide both side by $g$ resulting in a equivalent equation.<br /><br />Let $a = \frac{A}{g}$, $b = \frac{B}{g}$ and $c = \frac{C}{g}$. Then, $$(Ax + By = C) \: \equiv \: (ax + by = c)$$After simplification, <span style="background-color: yellow;">$gcd(a,b)$ is either $1$ or $-1$. If it is $-1$, then we need multiply $-1$ with $a,b$ and $c$ so that $gcd(a,b)$ becomes $1$</span> and the equation remains unchanged. Why did we make the $gcd(a,b)$ positive? You will find the reason below.<br /><br /><h2>Using Extended Euclidean Algorithm</h2>Recall that in a previous post "<a href="http://forthright48.blogspot.com/2015/07/extended-euclidean-algorithm.html">Extended Euclidean Algorithm</a>", we learned how to solve the Bezout's Identity $Ax + By = gcd(A, B)$. Can we apply that here in any way?<br /><br />Yes. Using $\text{ext_gcd()}$ function, we can find Bezout's coefficient for $ax + by = gcd(a,b)$. But we need to find solution for $ax + by = c$. Note that $gcd(a,b) = 1$, so when we use $\text{ext_gcd()}$ we find a solution for $ax + by = 1$. Let this solution be $(x_1,y_1)$. We can extend this solution to solve our original problem.<br /><br />Since we already have a solution where $ax_1 + by_1 = 1$, multiplying both sides with $c$ gives us $a(x_1c) + b(y_1c) = c$. So our result is $(x,y) = (x_1c, y_1c)$. This is why we had to make sure that $gcd(a,b)$ was $1$ and not $-1$. Otherwise, multiplying $c$ would have resulted $ax + by = -c$ instead.<br /><br /><h2>Summary of Solution</h2>Here is a quick summary of what I described above. We can find solution for Linear Diophantine Equation $Ax + By = C$ in 3 steps:<br /><ol><li><b><u>No Solution</u></b><br />First check if solution exists for given equation. Let $g = gcd(A,B)$. If $g \not| \: C$ then no solution exists.<br /><br /></li><li><b><u>Simplify Equation</u></b><br />Let $a = \frac{A}{g}, b = \frac{B}{g}$ and $c = \frac{C}{g}$. Then finding solution for $Ax + By = C$ is same as finding solution for $ax + by = c$. In simplified equation, make sure $GCD(a,b)$ is $1$. If not, multiply $-1$ with $a,b,c$.<br /><br /></li><li><b><u>Extended Euclidean Algorithm</u></b><br />Use $\text{ext_gcd()}$ to find solution $(x_1,y_1)$ for $ax + by = 1$. Then multiply the solution with $c$ to get solution for $ax + by = c$, where $x = x_1 \times c, y = y_1 \times c$.</li></ol>Let us try few examples.<br /><br /><h3>Example 1: $2x + 3y = 7$</h3><div><br /><b>Step $1$:</b> $g = GCD(2,3) = 1$. Since $1$ divides $7$, solution exists.</div><div><b>Step $2$:</b> Since $g$ is already $1$ there is nothing to simplify.</div><div><b>Step $3$:</b> Using $\text{ext_gcd()}$ we get $(x,y) = (-1,1)$. But this is for $ax + by = 1$. We need to multiply $7$. So our solution is $(-7,7)$.</div><div><br /></div><div>$2 \times -7 + 3 \times 7 = -14 + 21 = 7$. The solution is correct.<br /><br /><h3>Example 2: $4x + 10y = 8$</h3></div><div><br /></div><div>Step $1$: $g = GCD(4,10) = 2$. Since $2$ divides $8$, solution exists.</div><div>Step $2$: $a = \frac{4}{2}, b = \frac{10}{2}, c = \frac{8},{2}$. We will find solution of $2x + 5y = 4$.</div><div>Step $3$: Using $\text{ext_gcd()}$ we get $(x,y) = (-2,1)$. But this is for $ax + by = 1$. We need to multiply $4$. So our solution is $(-8,4)$.</div><div><br /></div><div>$ax + by = 2 \times -8 + 5 \times 4 = -16 + 20 = 4 = c$.</div><div>Also, $Ax + By = 4 \times -8 + 10 \times 4 = -32 + 40 = 8 = C$. The solution is correct. Both $ax + by = c$ and $Ax + By = C$ are satisfied.</div><br /><h2>Finding More Solutions</h2><div>We can now find a possible solution for $Ax + By = C$, but what if we want to find more? How many solutions are there? Since the solution for $Ax + By = C$ is derived from Bezout's Identity, there are infinite solutions.</div><div><br /></div><div>Suppose we found a solution $(x,y)$ for $Ax + By = C$. Then we can find more solutions using the formula: $( x + k \frac{B}{g}, y - k \frac {A}{g})$, where $k$ is any integer. </div><div><br /></div><h2>Code</h2>Let us convert our idea into code.<br /><pre class="brush:cpp;">bool linearDiophantine ( int A, int B, int C, int *x, int *y ) {<br /> int g = gcd ( A, B );<br /> if ( C % g != 0 ) return false; //No Solution<br /><br /> int a = A / g, b = B / g, c = C / g;<br /> ext_gcd( a, b, x, y ); //Solve ax + by = 1<br /><br /> if ( g < 0 ) { //Make Sure gcd(a,b) = 1<br /> a *= -1; b *= -1; c *= -1;<br /> }<br /><br /> *x *= c; *y *= c; //ax + by = c<br /> return true; //Solution Exists<br />}<br /><br />int main () {<br /> int x, y, A = 2, B = 3, C = 5;<br /> bool res = linearDiophantine ( A, B, C, &x, &y );<br /><br /> if ( res == false ) printf ( "No Solution\n" );<br /> else {<br /> printf ( "One Possible Solution (%d %d) \n", x, y );<br /><br /> int g = gcd ( A, B );<br /><br /> int k = 1; //Use different value of k to get different solutions<br /> printf ( "Another Possible Solution (%d %d)\n", x + k * ( B / g ), y - k * ( A / g ) );<br /> }<br /><br /> return 0;<br />}<br /></pre>$\text{linearDiophantine}()$ function finds a possible solution for equation $Ax + By = C$. It takes in $5$ parameters. $A,B,C$ defines the coefficients of equation and $*x, *y$ are two pointers that will carry our solution. The function will return $true$ if solution exists and $false$ if not.<br /><br />In line $2$ we calculate $gcd(A,B)$ and in line $3$ we check if $C$ is divisible by $g$ or not. If not, we return $false$.<br /><br />Next on line $5$ we define $a,b,c$ for simplified equation. On line $6$ we solve for $ax + by = 1$ using $\text{ext_gcd}$. Then we check if $g < 0$. If so, we multiply $-1$ with $a,b,c$ to make it positive. Then we multiply $c$ with $x,y$ so that our solution satisfies $ax + by = c$. A solution is found so we return true.<br /><br />In $\text{main}()$ function, we call $\text{linearDiophantine}()$ using $A=2,B=3,C=5$. In line $22$ we print a possible solution. In line $27$ we print another possible solution using formula for more solutions.<br /><br /><h2>$A$ and $B$ with Value $0$</h2><div>Till now we assumed $\{A, B\}$ have non-zero values. What happens if they have value $0$?</div><div><br /></div><h3>When Both $A = B = 0$</h3><div>When both $A$ and $B$ are zero, the value of $Ax + By$ will always be $0$. Therefore, if $C \neq 0$ then there is no solution. Otherwise, any pair of value for $(x,y)$ will act as a solution for the equation.</div><div><br /></div><h3>When $A$ or $B$ is $0$</h3><div>Suppose only $A$ is $0$. Then equation $Ax + By = C$ becomes $0x + By = C \: \equiv \: By = C$. Therefore $y = \frac {C}{B}$. If $B$ does not divide $C$ then there is no solution. Else solution will be $(x,y) = (\frac{C}{B}, k)$, where $k$ is any intger.</div><div><br /></div><div>Using same logic, when $B$ is $0$, solution will be $(x,y) = ( k, \frac{C}{A} )$.</div><div><br /></div><h2>Coding Pitfalls</h2>When we use $gcd(a,b)$ in our code, we mean the result from Euclidean Algorithm, not what we understand mathematically. $gcd(4,-2)$ is $-2$ according to Euclidean Algorithm whereas it is $2$ in common sense.<br /><br /><h2>Resources</h2><div><ol><li>Wiki - Diophantine Equation - <a href="https://en.wikipedia.org/wiki/Diophantine_equation">https://en.wikipedia.org/wiki/Diophantine_equation</a></li><li>forthright48 - Extended Euclidean Algorithm - <a href="http://forthright48.blogspot.com/2015/07/extended-euclidean-algorithm.html">http://forthright48.blogspot.com/2015/07/extended-euclidean-algorithm.html</a></li></ol></div><br />Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-36783778225374763242015-07-26T14:18:00.000+06:002015-09-23T17:10:55.176+06:00Extended Euclidean AlgorithmExtended Euclidean Algorithm is an extension of <a href="http://forthright48.blogspot.com/2015/07/euclidean-algorithm-finding-greatest.html">Euclidean Algorithm</a> which finds two things for integer $a$ and $b$:<br /><ol><li>It finds the value of $GCD(a,b)$</li><li>It finds two integers $x$ and $y$ such that, $ax + by = gcd(a,b)$. </li></ol><div>The expression $ax + by = gcd(a,b)$ is known as Bezout's identity and the pair $(x,y)$ that satisfies the identity is called Bezout coefficients. We will look into Bezout's identity at the end of this post. For now just know the name.<br /><br /><h2>How It Works</h2></div><div>In Euclidean Algorithm we worked with remainders only. In Extended Euclidean Algorithm (ext_gcd) we will use the quotient and few other extra variables. Suppose we are finding $GCD(240,46)$. Using Euclidean Algorithm the steps for finding $GCD$ would be:</div><div><br /></div><div>$GCD(240,46) = GCD(46, 240 \: \% \: 46) = GCD(46,10)$<br />$GCD(46,10) = GCD(10, 46 \: \% \: 10) = GCD(10,6)$</div><div>$GCD(10,6) = GCD(6, 10 \: \% \: 6) = GCD(6,4)$</div><div>$GCD(6,4) = GCD(4, 6 \: \% \: 4) = GCD(4,2)$</div><div>$GCD(4,2) = GCD(2, 4 \: \% \: 2) = GCD(2,0)$</div><div>$GCD(2,0) = 2$</div><div><br /></div><h3>Introducing $r$ and $q$</h3><div>We will slowly move towards ext_gcd algorithm. Let us add two new variables. $r$ represents remainder and $q$ means quotient. </div><br />Let $r_0$ be $a$ and $r_1$ be $b$. At each step, we will calculate $r_i$. Let $q_i = \lfloor \frac{r_{i-2}}{ r_{i-1}} \rfloor$. Therefore, $r_i = r_{i-2} - q_i \times r_{i-1}$.<br /><br />Then the above steps will look like the following:<br /><div class="postTable"><table><thead><tr><th>Index i</th> <th>Quotient $q_i$</th> <th>Remainder $r_i$</th> </tr></thead> <tbody><tr> <td>$0$</td> <td></td> <td>$240$</td></tr><tr> <td>$1$</td> <td></td> <td>$46$</td></tr><tr> <td>$2$</td> <td>$240 / 46 = 5$</td> <td>$240 - 5 \times 46 = 10$</td></tr><tr> <td>$3$</td> <td>$46 / 10 = 4$</td> <td>$46 - 4 \times 10 = 6$</td></tr><tr> <td>$4$</td> <td>$10 / 6 = 1$</td> <td>$10 - 1 \times 6 = 4$</td></tr><tr> <td>$5$</td> <td>$6 / 4 = 1$</td> <td>$6 - 1 \times 4 = 2$</td></tr><tr> <td>$6$</td> <td>$4 / 2 = 2$</td> <td>$4 - 2 \times 2 = 0$</td></tr></tbody> </table></div>This table is same as the calculation for Euclidean algorithm except for few extra details. Note that the line before last ( index $5$ ) contains the $GCD(a,b)$.<br /><br /><h3>Introducing $x$ and $y$</h3><div>We are trying to express $GCD(a,b) = ax + by$. So the variable $x$ and $y$ will hold the coefficients. To be exact we will write each row as terms of $a$ and $b$, i.e, $r_i = ax_i + by_i$.</div><div><br /></div><div>Initially, $(x_0, y_0) = (1,0)$ and $(x_1, y_1) = (0,1)$. But how do we calculate $(x_i,y_i)$?</div><div><br /></div><div>We know that $r_i = r_{i-2} - q_i \times r_{i-1}$. We also claimed that $r_i = ax_i + by_i$. By combining this two we get :</div><div><br /></div><div>$r_i = ( ax_{i-2} + by_{i-2} ) - q_i \times ( ax_{i-1} + by_{i-1} )$</div><div>$r_i = ax_{i-2} - q_i \times ax_{i-1} + by_{i-2} - q_i \times by_{i-1}$</div><div>$r_i = a ( x_{i-2} - q_i \times x_{i-1} ) + b ( y_{i-2} - q_i \times y_{i-1})$</div><div>$r_i = a x_i + b y_i$</div><div>$\therefore x_i = x_{i-2} - q_i \times x_{i-1} \: \text{and} \: y_i = y_{i-2} - q_i \times y_{i-1}$.</div><br />Our new table enhanced with $x$ and $y$ will now look like:<br /><div class="postTable"><table><thead><tr><th>Index i</th> <th>Quotient $q_i$</th> <th>Remainder $r_i$</th> <th>$x_i$</th> <th>$y_i$</th> </tr></thead> <tbody><tr> <td>$0$</td> <td></td> <td>$240$</td> <td>$1$</td> <td>$0$</td> </tr><tr> <td>$1$</td> <td></td> <td>$46$</td> <td>$0$</td> <td>$1$</td> </tr><tr> <td>$2$</td> <td>$240 / 46 = 5$</td> <td>$240 - 5 \times 46 = 10$</td> <td>$1 - 5 \times 0 = 1$</td> <td>$0 - 5 \times 1 = -5$</td> </tr><tr> <td>$3$</td> <td>$46 / 10 = 4$</td> <td>$46 - 4 \times 10 = 6$</td> <td>$0 - 4 \times 1 = -4 $</td> <td>$1 - 4 \times -5 = 21$</td> </tr><tr> <td>$4$</td> <td>$10 / 6 = 1$</td> <td>$10 - 1 \times 6 = 4$</td> <td>$1 − 1 × −4 = 5$</td> <td>$−5 − 1 × 21 = −26$</td> </tr><tr> <td>$5$</td> <td>$6 / 4 = 1$</td> <td>$6 - 1 \times 4 = 2$</td> <td>$−4 − 1 × 5 = −9$</td> <td>$21 − 1 × −26 = 47$</td> </tr><tr> <td>$6$</td> <td>$4 / 2 = 2$</td> <td>$4 - 2 \times 2 = 0$</td> <td>$5 − 2 × -9 = 23$</td> <td>$−26 − 2 × 47 = −120$</td> </tr></tbody> </table></div>Our answer lies on the line before last. $240 \times -9 + 46 \times 47 = 2$.<br /><br />So all we need to do now is implement these steps in code.<br /><br /><h2>Code</h2><div>Even though we will be calculating many rows in ext_gcd algorithm, in order to calculate any row we just need information from previous two rows. So we can save memory by simply storing $r_{i-2}, r_{i-1}, x_{i-2}, x_{i-1}, y_{i-2}, y_{i-1}$.</div><div><br /></div><div>In our code, $x2$ is $x_{i-2}$, $x1$ is $x_{i-1}$ and $x$ is $x_i$. Same style is followed for other variable.</div><pre class="brush: cpp;">int ext_gcd ( int A, int B, int *X, int *Y ){<br /> int x2, y2, x1, y1, x, y, r2, r1, q, r;<br /> x2 = 1; y2 = 0;<br /> x1 = 0; y1 = 1;<br /> for (r2 = A, r1 = B; r1 != 0; r2 = r1, r1 = r, x2 = x1, y2 = y1, x1 = x, y1 = y ) {<br /> q = r2 / r1;<br /> r = r2 % r1;<br /> x = x2 - (q * x1);<br /> y = y2 - (q * y1);<br /> }<br /> *X = x2; *Y = y2;<br /> return r2;<br />}<br /></pre>In line $1$ we define the function. The function $ext\_ gcd()$ takes 4 parameters. The first two parameter $A$ and $B$ represents the two integers whose gcd we wish to find. The next two parameter $*X$ and $*Y$ are two integer pointers. Our function returns us $GCD(A,B)$ but we also want to find the coefficients $x,y$ such that $ax + by = GCD(A,B)$. So we extract those values using pointers.<br /><br />In line $2$, we declare all necessary variables. We initiate some variables in line $3$ and $4$. Then a for loop starts in line $5$. We initiate few more variables in first section, $r2$ and $r1$. The loop runs till $r1$ becomes 0. In the last section of for loop, we make transitions of variables, such as $r2$ becomes $r1$ and $r1$ gets the new value $r$.<br /><br />In side the for loop we calculate values needed for current row, $q, r, x, y$.<br /><br />When the loop finishes, $r1$ contains $0$ and $r2$ is the row before it containing $GCD(A,B)$. So we send $x2,y2$ as coefficients.<br /><br /><h2>Complexity</h2><div>Same as Euclidean Algorithm.</div><h2>Uniqueness of Solution</h2><div>Using ext_gcd we found $(x,y)$ pair which satisfies $ax + by = gcd(a,b)$. But is it unique?</div><div><br /></div><div>No. Once we find a pair $(x,y)$ using ext_gcd, we can generate infinite pairs of Bezout coefficients using the formula: $$( x + k \frac { b } { \text{gcd}(a,b)}, y - k \frac { a } { \text{gcd}(a,b)} )$$Using any integer value for $k$ we can generate a different pair of Bezout coefficients $(x,y)$ which will satisfy the Bezout's identity. Here is why it works:</div><div><br /></div><div>$a ( x + \frac { kb } { \text{gcd}(a,b)} ) + b ( y - \frac { ka } { \text{gcd}(a,b)} )$</div><div>$ax + \frac { kab } { \text{gcd}(a,b)} + by - \frac { kab } { \text{gcd}(a,b)}$</div><div>$ax + by$</div><div><br /></div><div>As you can see above, the terms with $k$ in them cancel each other out. They don't change the expression $ax + by$ in any way, therefore, there are infinite pairs of Bezout coefficients.</div><div><br /></div><h2>Smallest Positive Integer of form $ax + by$</h2><div>We showed that it is possible to write $ax + by = gcd(a,b)$, now we need to find a positive number smaller than $gcd(a,b)$ of form $ax + by$. Is that even possible?</div><div><br />No. $gcd(a,b)$ divides both $a$ and $b$. Therefore, if a number is of form $ax + by$ it will be divisible by $gcd(a,b)$ since $ax + by$ is divisible by $gcd(a,b)$. And the smallest positive number which is divisible by $gcd(a,b)$ is $gcd(a,b)$ itself.</div><div><br /></div><div>So $gcd(a,b)$ is the smallest positive number of the form $ax + by$.<br /><br /></div><h2>Bezout's Identity</h2><div>This was mentioned at the beginning of the post. Almost everything related to Bezout's Identity has been explained. I will still mention them once more for the sake of completeness.</div><br /><blockquote class="tr_bq">Bézout's identity (also called Bézout's lemma) is a theorem in the elementary theory of numbers: let a and b be nonzero integers and let d be their greatest common divisor. Then there exist integers x and y such that $ax + by = d$<br /><br />In addition,<br /><ul><li>the greatest common divisor d is the smallest positive integer that can be written as $ax + by$</li><li>every integer of the form ax + by is a multiple of the greatest common divisor d. </li></ul>-<a href="https://en.wikipedia.org/wiki/B%C3%A9zout%27s_identity">Wiki</a></blockquote>Bezout's Lemma simply states that $ax + by = gcd(a,b)$ exists. We need to use Extended Euclidean Algorithm to find Bezout's Coefficients.<br /><br /><h2>Coding Pitfalls</h2><div>Careful about the equation $ax + by = \text{gcd} (a,b)$. Here $\text{gcd}()$ means the result from Euclidean Algorithm, not what we mean mathematically. For example, when $a = 4$ and $b = -2$, Extended Euclidean Algorithm finds solution for $4x - 2y = -2$. According to Euclidean algorithm $gcd(4,-2)$ returns $-2$, though in common sense it should be $2$. </div><div><br /></div><h2>Reference</h2><div><ol><li>forthright48 - Euclidean Algorithm - <a href="http://forthright48.blogspot.com/2015/07/euclidean-algorithm-finding-greatest.html">http://forthright48.blogspot.com/2015/07/euclidean-algorithm-finding-greatest.html</a></li><li>Wiki - Extended Euclidean algorithm - <a href="https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm">https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm</a></li><li>Wiki - Bezout's identity - <a href="https://en.wikipedia.org/wiki/B%C3%A9zout%27s_identity">https://en.wikipedia.org/wiki/B%C3%A9zout%27s_identity</a></li></ol></div><div><br /></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com1tag:blogger.com,1999:blog-811371107174953309.post-46385211785773107872015-07-25T15:46:00.000+06:002015-07-25T18:03:09.197+06:00Introduction to Modular Arithmetic<blockquote class="tr_bq">"In mathematics, modular arithmetic is a system of arithmetic for integers, where numbers "wrap around" upon reaching a certain value—the modulus." - <a href="https://en.wikipedia.org/wiki/Modular_arithmetic">Wiki</a></blockquote>The clock is a good example for modular arithmetic. Let's say $12$'o clock means $0$. Then clock shows the following values, $\{0,1,2,3,4,5,6,7,8,9,10,11,0,1,2,3...\}$. Every time clock hits $12$, it wraps around to $0$. Since clock wraps around $12$, we say that $12$ is the Modulus.<br /><br /><h2>General Concepts </h2><div>In order to understand modular arithmetic, we need to know the following:<br /><br /></div><h3>Modulo Operation $\%$</h3>Modular arithmetic deals with Modulo operation. It has the symbol, "$\%$". The expression $A \: \% \: B = C$ means that when we divide $A$ with $B$ we get the remainder $C$. In other words, the modulo operation finds the remainder. For example:<br /><br />$5 \: \% \: 3 = 2$<br />$15 \: \% \: 5 = 0$<br />$24 \: \% \: 30 = 24$<br />$30 \: \% \: 24 = 6$<br />$4 \: \% \: 4 = 0$<br /><br /><h3>Congruence Relation $\equiv$</h3><div>In mathematics, saying that $A \: \equiv \: B \: \text{(mod M)}$ means that both $A$ and $B$ has same remainder when they are divided by $M$. For example,</div><div><br /></div><div>$5 \: \equiv \: 10 \text{ (mod 5) }$</div><div>$2 \: \equiv \: 7 \text{ (mod 5) }$</div><div>$4 \: \equiv \: 10 \text{ (mod 3) }$</div><div><br /><h3>Range</h3></div><div>In modular arithmetic with modulus $M$, we only work with values ranging from $0$ to $M-1$. All other values can be converted to the corresponding value from the range using Modulo operation. For example for $M=3$:<br /><br /></div><div>$<br />\begin{matrix}<br />\text{Normal Arithmetic } & -3 & -2 & -1 & 0 & 1 & 2 & 3\\<br />\text{Modular Arithmetic} & 0 & 1 & 2 & 0 & 1 & 2 & 0<br />\end{matrix}<br />$<br /><br /><h2>Properties of Modular Arithmetic</h2><div>There are four properties of Modular Arithmetic that we need to learn. </div><div><br /><h3>Addition </h3></div><div>$$ \text{if,} \: a \: \equiv \: b \: \text{(mod m)} \\ \text{and,} \: c \: \equiv \: d \: \text{(mod m)} \\ \text{then,} \: a + c \: \equiv \: b + d \text{(mod m)}$$</div><div><br /></div><div>Meaning, if we have to find $(a_1+a_2+a_3...+a_n) \: \% \: m$, we can instead just find $ ( (a_1 \: \% \: m ) + (a_2 \: \% \: m ) +... (a_n \: \% \: m ) ) \: \% \: m$.</div><div><br /></div><div>For example, $ ( 13 + 24 + 44 ) \: \% \: 3 = ( 1 + 0 + 2 ) \: \% \: 3 = 3 \: \% \: 3 = 0$. </div><div><br /></div><h3>Subtraction</h3><div><div>$$ \text{if,} \: a \: \equiv \: b \: \text{(mod m)} \\ \text{and,} \: c \: \equiv \: d \: \text{(mod m)} \\ \text{then,} \: a - c \: \equiv \: b - d \text{(mod m)}$$</div><div><br /></div><div>Meaning, if we have to find $(a_1-a_2-a_3...-a_n) \: \% \: m$, we can instead just find $ ( (a_1 \: \% \: m ) - (a_2 \: \% \: m ) -... (a_n \: \% \: m ) ) \: \% \: m$.</div><div><br /></div><div>For example, $ ( -14 - 24 - 44 ) \: \% \: 3 = ( -2 - 0 - 2 ) \: \% \: 3 = -4 \: \% \: 3 = -1 \: \% \: 3 = 2$.</div></div><div><br /></div><h3>Multiplication</h3><div>$$ \text{if,} \: a \: \equiv \: b \: \text{(mod m)} \\ \text{and,} \: c \: \equiv \: d \: \text{(mod m)} \\ \text{then,} \: a \times c \: \equiv \: b \times d \: \text{(mod m)}$$</div><div><br /></div><div><div>Meaning, if we have to find $(a_1 \times a_2 \times a_3... \times a_n) \: \% \: m$, we can instead just find $ ( (a_1 \: \% \: m ) \times (a_2 \: \% \: m ) \times ... (a_n \: \% \: m ) ) \: \% \: m$.</div><div><br /></div><div>For example, $ ( 14 \times 25 \times 44 ) \: \% \: 3 = ( 2 \times 1 \times 2 ) \: \% \: 3 = 4 \: \% \: 3 = 1 \: \% \: 3 = 2$.</div></div><div><br /></div><h3>Division</h3><div>Modular Division does not work same as the above three. We have to look into a separate topic known as Modular Inverse in order to find $\frac{a}{b} \: \% \: m$. Modular Inverse will be covered in a separate post.<br /><br />For now just know that, $\frac{a}{b} \: \not\equiv \: \frac{a \: \% \: m} {b \: \% \: m} \: \text{(mod m)}$. For example, $\frac{6}{3} \: \% \: 3$ should be $2$. But using the formula above we get $\frac{0}{0} \: \% \: 3$ which is undefined.</div><div><br /></div><h2>Coding Tips</h2><div><br /></div><h3>Handling Negative Numbers</h3><div>In some programming language, when modulo operation is performed on negative numbers the result is also negative. For example in C++ we will get $-5 \: \% \: 4 = -1$. But we need to convert this into legal value. We can convert the result into legal value by adding $M$.</div><div><br /></div><div>Therefore, in C++ it is better to keep a check for negative number when doing modulo operation.</div><div><pre class="brush: cpp;">res = a % m;<br />if ( res < 0 ) res += m;<br /></pre><br /><h3>Modulo Operation is Expensive</h3><div>Modulo Operation is as expensive as division operator. It takes more time to perform division and modulo operations than to perform addition and subtraction. So it is better to avoid using modulo operation where possible.</div><div><br /></div><div>One possible situation is when we are adding lots of values. Let's say that we have two variables $a$ and $b$ both between $0$ to $m-1$. Then we can write:</div><div><pre class="brush: cpp;">res = a + b;<br />if ( res >= m ) res -= m;<br /></pre>This technique will not work when we are multiplying two values.<br /><br /><h2>Resource</h2><div><ol><li>Wiki - Modular Arithmetic - <a href="https://en.wikipedia.org/wiki/Modular_arithmetic">https://en.wikipedia.org/wiki/Modular_arithmetic</a></li><li>AoPS - Introduction to Modular Arithmetic - <a href="http://www.artofproblemsolving.com/wiki/index.php/Modular_arithmetic/Introduction">http://www.artofproblemsolving.com/wiki/index.php/Modular_arithmetic/Introduction</a></li></ol></div><br /></div></div></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0tag:blogger.com,1999:blog-811371107174953309.post-12380587721217905202015-07-20T23:42:00.000+06:002015-07-29T22:04:06.958+06:00Divisor Summatory Function<h2>Problem</h2>Given an integer $N$, find the Sum of <a href="http://forthright48.blogspot.com/2015/07/number-of-divisors-of-integer.html">Number of Divisors</a> from $1\: to \: N$. That is, find $\sum_{i=1}^{N}NOD(i)$.<br /><br />For example, find the Sum of Number of Divisors, $SNOD()$, when $N=5$. $SNOD(5)= NOD(1) + NOD(2) + NOD(3) + NOD(4) + NOD(5) = 1 + 2 + 2 + 3 + 2 = 10$.<br /><br /><h2>Divisor Summatory Function</h2><blockquote>"In number theory, the divisor summatory function is a function that is a sum over the divisor function." - <a href="https://en.wikipedia.org/wiki/Divisor_summatory_function">Wiki</a></blockquote>Divisor Summatory Function is defined as $D()$ in Wiki, but we will use $SNOD()$ in this article. Divisor Summatory Function finds Sum of Number of Divisors from $1\: to \: N$.<br /><br /><h2>Brute Force Solution - $O(N^2)$ Approach</h2>A simple solution is to run two nested loop which counts all divisors from $1 \: to \: N$. It is simple but inefficient.<br /><pre class="brush:cpp;">int SNOD( int n ) {<br /> int res = 0;<br /> for ( int i = 1; i <= n; i++ ) { //For each number from 1 - N<br /> for ( int j = 1; j <= i; j++ ) { //Find possible divisors of "i"<br /> if ( i % j == 0 ) res++;<br /> }<br /> }<br /> return res;<br />}<br /></pre><h2>Using $NOD()$ - $O( N \times \frac{\sqrt{N}}{ln(\sqrt{N})} ) $ Approach</h2>We can use the code for $NOD()$ to get the number of divisors for each integer from $1$ to $N$ instead of using a $O(N)$ loop. This will make it faster than the above solution.<br /><pre class="brush:cpp;highlight: [4]">int SNOD( int n ) {<br /> int res = 0;<br /> for ( int i = 1; i <= n; i++ ) {<br /> res += NOD(i);<br /> }<br /> return res;<br />}<br /></pre>Since $NOD()$ has same complexity as $factorize()$, the complexity for this code is $O( N \times \frac{\sqrt{N}}{ln(\sqrt{N})} ) $.<br /><br /><h2>Using Divisor List - $O(N)$ Approach</h2>Suppose we are trying to find $SNOD(10)$. Let us write down the divisors for each integer between $1$ and $10$.<br /><br />$NOD(1) = 1, \: \{1\}$<br />$NOD(2) = 2, \: \{1 , 2\}$<br />$NOD(3) = 2, \: \{1 , 3\}$<br />$NOD(4) = 3, \: \{1, 2, 4\}$<br />$NOD(5) = 2, \: \{1, 5\}$<br />$NOD(6) = 4, \: \{1, 2, 3, 6\}$<br />$NOD(7) = 2, \: \{1, 7\}$<br />$NOD(8) = 4, \: \{1, 2, 4, 8\}$<br />$NOD(9) = 3, \: \{1, 3, 9\}$<br />$NOD(10) = 4, \: \{1, 2, 5, 10\}$<br /><br />So $SNOD(10) = 1 + 2 + 2 + 3 + 2 + 4 + 2 + 4 + 3 + 4 = 27$.<br /><br />Now, look at the divisor list for each of integer between $1$ and $N$. Each divisor in the divisor lists contributes $1$ to the result. There are $27$ divisors in total.<br /><br />How many times do we find $1$ in the divisor lists? It will appear once for every integer it can divide from $1$ to $N$. $1$ can divide $\frac{N}{1}$ numbers. So it appears $\frac{10}{1} = 10$ times.<br /><br />Repeat for $2 \: to \: N$. How many times do we find $2$? $2$ can divide $\frac{10}{2} = 5$ numbers, namely $\{2,4,6,8,10\}$. So it will appear $5$ times.<br /><br />By repeating this from $1$ to $N$ we can find $SNOD(N)$.<br /><pre class="brush:cpp;highlight:[4]">int SNOD( int n ) {<br /> int res = 0;<br /> for ( int i = 1; i <= n; i++ ) {<br /> res += n / i;<br /> }<br /> return res;<br />}<br /></pre>We removed $NOD()$ from the code in line $4$ and replaced it with $\frac{n}{i}$.<br /><br /><h2>Using Divisor Pairs - $O(\sqrt{N})$ Approach</h2>Let us create another list. Instead of making divisor list, we will make "Ordered Divisor Pair" list for each value between $1$ to $N$. A divisor pair of value $X$ is a pair $(A, B)$ such that $A \times B = X$. Ordered pair means if $A \neq B$, then $(A, B)$ is a different pair than $(B, A)$.<br /><br />$NOD(1) = 1, \: (1,1)$<br />$NOD(2) = 2, \: (1,2) (2,1)$<br />$NOD(3) = 2, \: (1,3)(3,1)$<br />$NOD(4) = 3, \: (1,4)(2,2)(4,1)$<br />$NOD(5) = 2, \: (1,5)(5,1)$<br />$NOD(6) = 4, \: (1,6)(2,3)(3,2)(6,1)$<br />$NOD(7) = 2, \: (1,7)(7,1)$<br />$NOD(8) = 4, \: (1,8)(2,4)(4,2)(8,1)$<br />$NOD(9) = 3, \: (1,9)(3,3)(9,1)$<br />$NOD(10) = 4, \: (1,10)(2,5)(5,2)(10,1)$<br /><br />The "Ordered Divisor Pair" list is almost same as "Divisor List". The only difference is that each divisor $A$ is now paired with $B = \frac{X}{A}$.<br /><br />Now, $NOD(X)$ represents number of ordered divisor pairs $(A,B)$ such that $A \times B = X$. So, <span style="background-color: yellow;">$SNOD(N)$ is same as number of ordered divisor pairs such that $A \times B \leq N$</span>.<br /><br /><h3>How to find Number of Ordered Divisor Pairs $\leq N$?</h3><div>We can find this by following three steps.</div><br /><h4><b>1. Find number of divisor pairs $(A, B)$ where $A < B$</b></h4>What are the possible values for $A$? Can it ever be $> \sqrt{N}$? No. If $A > \sqrt{N}$, then with $B>A$, $B$ will be $> \sqrt{N}$ and $A \times B$ will be $ > N$. So $A$ can only be between $1$ and $\lfloor \sqrt{N} \rfloor$.<br /><br />Let $u = \lfloor \sqrt{N} \rfloor$. Then for each value of $A$ between $1$ and $u$, what are the possible values of $B$?<br /><br />For any value of $A$, there are $\frac{N}{A}$ possible values of $B$ for which $A \times B$ does not exceed $N$. For $N=10$ and $A=2$, there are $\frac{10}{2}=5$ possible values for $B$, and they are $\{1,2,3,4,5\}$. Multiplying $A=2$ with any of these values produces value $\leq N$. But we need $B > A$. So we omit $\{1,2\}$ from the list. We can only use $\{3,4,5\}$ as $B$.<br /><br />So, for any value of $A$, there will be $\lfloor \frac{N}{A} \rfloor$ possible values of $B$. The values will be from $1$ to $\lfloor \frac{N}{A} \rfloor$. Out of those, we have to omit values that are $\leq A$. So we can use $\lfloor \frac{N}{A} \rfloor - A$ values for $B$.<br /><br />For example, $N=10$. Then $u=\lfloor \sqrt{10} \rfloor = 3$. When $A=1$, number of legal values we can use is $\lfloor \frac{10}{1} \rfloor - 1 = 9$. Namely, $B$ can be one of $\{2,3,4,5,6,7,8,9,10\}$.<br /><br />When $A=2$, legal values are $\lfloor \frac{10}{2} \rfloor - 2 = 3$. $B$ can be $\{3,4,5\}$.<br />When $A=3$, legal values are $\lfloor \frac{10}{3} \rfloor - 3 = 0$.<br /><br />$A$ cannot be bigger than $u$. So number of $(A, B)$ pairs such that $A \times B \leq 10$ and $A < B$ is 12.<br /><br /><h4><b>2. Multiply the result with $2$</b></h4>Once we find the number of pairs $(A, B)$ such that $A < B$ and $A\times B \leq N$, we multiply that number with $2$. That is because we want to find the number of ordered divisor pairs. When $A \neq B$, each pair $(A, B)$ will occur in "Ordered Divisor Pair" list as $(B, A)$ again.<br /><br />Once we double the number, our result contains the number of ordered pair $(A, B)$ such that $A \neq B$ and $A \times B \leq N$.<br /><br /><h4><b>3. Add the number of pairs $(A, B)$ where $A = B$</b></h4>We still did not count pairs where $A = B$ and $A \times B \leq N$. How many pairs can there be? Since $A$ must be between $1$ and $u$, there cannot be more than $u$ such pairs. The pairs will be $(1,1),(2,2)...(u,u)$.<br /><br /><h3>Code for $O(\sqrt{N})$ Approach</h3><div><pre class="brush:cpp;">int SNOD( int n ) {<br /> int res = 0;<br /> int u = sqrt(n);<br /> for ( int i = 1; i <= u; i++ ) {<br /> res += ( n / i ) - i; //Step 1<br /> }<br /> res *= 2; //Step 2<br /> res += u; //Step 3<br /> return res;<br />}</pre></div><h2>Reference </h2><div><ol><li>forthright48 - Number of Divisors of an Integer - <a href="http://forthright48.blogspot.com/2015/07/number-of-divisors-of-integer.html">http://forthright48.blogspot.com/2015/07/number-of-divisors-of-integer.html</a></li><li>Wiki - Divisor Summatory Function - <a href="https://en.wikipedia.org/wiki/Divisor_summatory_function">https://en.wikipedia.org/wiki/Divisor_summatory_function</a></li></ol></div>Mohammad Samiul Islamhttps://plus.google.com/105318322076018763007noreply@blogger.com0