Wednesday, July 29, 2015

Simple Hyperbolic Diophantine Equation

Problem

Given the integers $A, B, C, D$, find a pair of $(x,y)$ such that it satisfies the equation $Axy + Bx + Cy = D$.

For example, $2xy + 5x + 56y = -7$ has a solution $(x,y) = (-27,64)$.

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

So how do we tackle the problem? First we will rewrite the equation to make it easier to approach.

Rewriting The Equation

$Axy + Bx + Cy = D$
$A^2xy + ABx + ACy = AD$, Multiply A with both sides
$A^2xy + ABx + ACy + BC = AD + BC$, Add BC to both sides
$Ax ( Ay + B ) + C ( Ay + B ) = AD + BC$,
$\therefore (Ax + C) (Ay + B) = AD + BC$

Now that we have rewritten the equation, we can see that the problem has two cases. We will handle each case separately.

When $AD + BC = 0$

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.

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.

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.

When $AD + BC \neq 0$

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.

$(Ax + C) = d_i$
$Ax = d_i - C$
$\therefore x = \frac{d_i - C}{A}$

$(Ay+B) = \frac{P}{d_i}$
$Ay = \frac{P}{d_i} - B$
$Ay = \frac{P - Bd_i}{d_i}$
$\therefore y = \frac{P - Bd_i}{Ad_i}$

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

Careful when calculating divisors of $P$. In this problem, we have to include the negative divisors of $P$ in our calculation. For example, $4$ will have the following divisors $(1,2,4,-1,-2,-4)$.

Example

Find solutions for  $2xy + 2x + 2y = 4$.

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

Since $(2x + 2)(2y + 2 ) = 12$, we get the following:

$d_1 = 1: x = \frac{1-2}{2} = \frac{-1}{2}, y = \frac{12 - 2 \times 1}{2 \times 1} = \frac{10}{2} = 5$
$d_2 = -1: x = \frac{-1-2}{2} = \frac{-3}{2}, y = \frac{12 - 2 \times -1}{2 \times -1} = \frac{14}{-2} = -7$
$d_3 = 2: x = \frac{2-2}{2} = 0, y = \frac{12 - 2 \times 2}{2 \times 2} = \frac{8}{4} = 2$
$d_4 = -2: x = \frac{-2-2}{2} = -2, y = \frac{12 - 2 \times -2}{2 \times -2} = \frac{16}{-4} = -4$
$d_5 = 3: x = \frac{3-2}{2} = \frac{1}{2}, y = \frac{12 - 2 \times 3}{2 \times 3} = \frac{6}{6} = 1$
$d_6 = -3: x = \frac{-3-2}{2} = \frac{-5}{2}, y = \frac{12 - 2 \times -3}{2 \times -3} = \frac{18}{-1} = -18$
$d_7 = 4: x = \frac{4-2}{2} = 1, y = \frac{12 - 2 \times 4}{2 \times 4} = \frac{4}{8} = \frac{1}{2}$
$d_8 = -4: x = \frac{-4-2}{2} = -3, y = \frac{12 - 2 \times -4}{2 \times -4} = \frac{20}{-8} = \frac{-5}{2}$
$d_9 = 6: x = \frac{6-2}{2} = 2, y = \frac{12 - 2 \times 6}{2 \times 6} = \frac{0}{12} = 0$
$d_{10} = -6: x = \frac{-6-2}{2} = -4, y = \frac{12 - 2 \times -6}{2 \times -6} = \frac{24}{-12} = -2$
$d_{11} = 12: x = \frac{12-2}{2} = 5, y = \frac{12 - 2 \times 12 }{2 \times 12 } = \frac{-12}{24} = \frac{-1}{2}$
$d_{12} = -12: x = \frac{-12-2}{2} = -7, y = \frac{12 - 2 \times -12 }{2 \times -12 } = \frac{36}{-24} = \frac{-3}{2}$

So the solutions are $(0,2),(-2,-4),(2,0),(-4,-2)$.

Code

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$. 
bool isValidSolution ( int a, int b, int c, int p, int div ) {
    if ( ( ( div - c )% a ) != 0 ) return false; //x = (div - c) / a
    if ( ( (p-b*div) % (a*div) ) != 0 ) return false;// y = (p-b*div) /(a*div)
    return true;
}

int hyperbolicDiophantine ( int a, int b, int c, int d ) {
    int p = a * d + b * c;

    if ( p == 0 ) { //ad + bc = 0
        if ( -c % a == 0 ) return -1; //Infinite solutions (-c/a, k)
        else if ( -b % a == 0 ) return -1; //Infinite solutions (k, -b/a)
        else return 0; //No solution
    }
    else {
        int res = 0;

        //For each divisor of p
        int sqrtn = sqrt ( p ), div;
        for ( int i = 1; i <= sqrtn; i++ ) {
            if ( p % i == 0 ) { //i is a divisor

                //Check if divisors i,-i,p/i,-p/i produces valid solutions
                if ( isValidSolution(a,b,c,p,i) )res++;
                if ( isValidSolution(a,b,c,p,-i) )res++;
                if ( p / i != i ) { //Check whether p/i is different divisor than i
                    if ( isValidSolution(a,b,c,p,p/i) )res++;
                    if ( isValidSolution(a,b,c,p,-p/i) )res++;
                }
            }
        }

        return res;
    }
}
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$.

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

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.

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.

Resource

  1. alpertron - Methods to Solve $Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$

Monday, July 27, 2015

Linear Diophantine Equation

Problem

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

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

The problem above is a type of Diophantine problem. In the Diophantine problem, only integer solutions to an equation are required. Since $Ax + By = C$ is a linear equation, this problem is a Linear Diophantine Problem where we have to find a solution for a Linear Diophantine Equation.

For now, let us assume that $A$ and $B$ are non-zero integers.

Existence of Solution

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.

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.

Simplifying the Equation

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.

Let $a = \frac{A}{g}$, $b = \frac{B}{g}$ and $c = \frac{C}{g}$. Then, $$(Ax + By = C) \: \equiv \: (ax + by = c)$$After simplification, $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$ and the equation remains unchanged. Why did we make the $gcd(a,b)$ positive? You will find the reason below.

Using Extended Euclidean Algorithm

Recall that in a previous post "Extended Euclidean Algorithm", we learned how to solve the Bezout's Identity $Ax + By = gcd(A, B)$. Can we apply that here in any way?

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.

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.

Summary of Solution

Here is a quick summary of what I described above. We can find solution for Linear Diophantine Equation $Ax + By = C$ in 3 steps:
  1. No Solution
    First check if solution exists for given equation. Let $g = gcd(A,B)$. If $g \not| \: C$ then no solution exists.

  2. Simplify Equation
    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$.

  3. Extended Euclidean Algorithm
    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$.
Let us try few examples.

Example 1: $2x + 3y = 7$


Step $1$: $g = GCD(2,3) = 1$. Since $1$ divides $7$, solution exists.
Step $2$: Since $g$ is already $1$ there is nothing to simplify.
Step $3$: 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)$.

$2 \times -7 + 3 \times 7 = -14 + 21 = 7$. The solution is correct.

Example 2: $4x + 10y = 8$


Step $1$: $g = GCD(4,10) = 2$. Since $2$ divides $8$, solution exists.
Step $2$: $a = \frac{4}{2}, b = \frac{10}{2}, c = \frac{8},{2}$. We will find solution of $2x + 5y = 4$.
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)$.

$ax + by = 2 \times -8 + 5 \times 4 = -16 + 20 = 4 = c$.
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.

Finding More Solutions

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.

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. 

Code

Let us convert our idea into code.
bool linearDiophantine ( int A, int B, int C, int *x, int *y ) {
    int g = gcd ( A, B );
    if ( C % g != 0 ) return false; //No Solution

    int a = A / g, b = B / g, c = C / g;
    ext_gcd( a, b, x, y ); //Solve ax + by = 1

    if ( g < 0 ) { //Make Sure gcd(a,b) = 1
        a *= -1; b *= -1; c *= -1;
    }

    *x *= c; *y *= c; //ax + by = c
    return true; //Solution Exists
}

int main () {
    int x, y, A = 2, B = 3, C = 5;
    bool res = linearDiophantine ( A, B, C, &x, &y );

    if ( res == false ) printf ( "No Solution\n" );
    else {
        printf ( "One Possible Solution (%d %d) \n", x, y );

        int g = gcd ( A, B );

        int k = 1; //Use different value of k to get different solutions
        printf ( "Another Possible Solution (%d %d)\n", x + k * ( B / g ), y - k * ( A / g ) );
    }

    return 0;
}
$\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.

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

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.

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.

$A$ and $B$ with Value $0$

Till now we assumed $\{A, B\}$ have non-zero values. What happens if they have value $0$?

When Both $A = B = 0$

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.

When $A$ or $B$ is $0$

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) = (k, \frac{C}{B})$, where $k$ is any intger.

Using same logic, when $B$ is $0$, solution will be $(x,y) = ( \frac{C}{A}, k )$.

Coding Pitfalls

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.

Resources


Sunday, July 26, 2015

Extended Euclidean Algorithm

Extended Euclidean Algorithm is an extension of Euclidean Algorithm which finds two things for integer $a$ and $b$:
  1. It finds the value of $GCD(a,b)$
  2. It finds two integers $x$ and $y$ such that, $ax + by = gcd(a,b)$. 
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.

How It Works

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:

$GCD(240,46) = GCD(46, 240 \: \% \: 46) = GCD(46,10)$
$GCD(46,10) = GCD(10, 46  \: \% \:  10) = GCD(10,6)$
$GCD(10,6) = GCD(6, 10  \: \% \:  6) = GCD(6,4)$
$GCD(6,4) = GCD(4, 6  \: \% \:  4) = GCD(4,2)$
$GCD(4,2) = GCD(2, 4  \: \% \:  2) = GCD(2,0)$
$GCD(2,0) = 2$

Introducing $r$ and $q$

We will slowly move towards ext_gcd algorithm. Let us add two new variables. $r$ represents remainder and $q$ means quotient. 

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

Then the above steps will look like the following:
Index i Quotient $q_i$ Remainder $r_i$
$0$ $240$
$1$ $46$
$2$ $240 / 46 = 5$ $240 - 5 \times 46 = 10$
$3$ $46 / 10 = 4$ $46 - 4 \times 10 = 6$
$4$ $10 / 6 = 1$ $10 - 1 \times 6 = 4$
$5$ $6 / 4 = 1$ $6 - 1 \times 4 = 2$
$6$ $4 / 2 = 2$ $4 - 2 \times 2 = 0$
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)$.

Introducing $x$ and $y$

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

Initially, $(x_0, y_0) = (1,0)$ and $(x_1, y_1) = (0,1)$. But how do we calculate $(x_i,y_i)$?

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 :

$r_i = ( ax_{i-2} + by_{i-2} ) - q_i \times  ( ax_{i-1} + by_{i-1} )$
$r_i = ax_{i-2} - q_i \times ax_{i-1} + by_{i-2} - q_i \times by_{i-1}$
$r_i = a ( x_{i-2} - q_i \times x_{i-1} ) + b ( y_{i-2} - q_i \times y_{i-1})$
$r_i = a x_i + b y_i$
$\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}$.

Our new table enhanced with $x$ and $y$ will now look like:
Index i Quotient $q_i$ Remainder $r_i$ $x_i$ $y_i$
$0$ $240$ $1$ $0$
$1$ $46$ $0$ $1$
$2$ $240 / 46 = 5$ $240 - 5 \times 46 = 10$ $1 - 5 \times 0 = 1$ $0 - 5 \times 1 = -5$
$3$ $46 / 10 = 4$ $46 - 4 \times 10 = 6$ $0 - 4 \times 1 = -4 $ $1 - 4 \times -5 = 21$
$4$ $10 / 6 = 1$ $10 - 1 \times 6 = 4$ $1 − 1 × −4 = 5$ $−5 − 1 × 21 = −26$
$5$ $6 / 4 = 1$ $6 - 1 \times 4 = 2$ $−4 − 1 × 5 = −9$ $21 − 1 × −26 = 47$
$6$ $4 / 2 = 2$ $4 - 2 \times 2 = 0$ $5 − 2 × -9 = 23$ $−26 − 2 × 47 = −120$
Our answer lies on the line before last. $240 \times -9 + 46 \times 47 = 2$.

So all we need to do now is implement these steps in code.

Code

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

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.
int ext_gcd ( int A, int B, int *X, int *Y ){
    int x2, y2, x1, y1, x, y, r2, r1, q, r;
    x2 = 1; y2 = 0;
    x1 = 0; y1 = 1;
    for (r2 = A, r1 = B; r1 != 0; r2 = r1, r1 = r, x2 = x1, y2 = y1, x1 = x, y1 = y ) {
        q = r2 / r1;
        r = r2 % r1;
        x = x2 - (q * x1);
        y = y2 - (q * y1);
    }
    *X = x2; *Y = y2;
    return r2;
}
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.

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

In side the for loop we calculate values needed for current row, $q, r, x, y$.

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.

Complexity

Same as Euclidean Algorithm.

Uniqueness of Solution

Using ext_gcd we found $(x,y)$ pair which satisfies $ax + by = gcd(a,b)$. But is it unique?

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:

$a ( x + \frac { kb } { \text{gcd}(a,b)} ) + b ( y - \frac { ka } { \text{gcd}(a,b)} )$
$ax + \frac { kab } { \text{gcd}(a,b)} + by - \frac { kab } { \text{gcd}(a,b)}$
$ax + by$

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.

Smallest Positive Integer of form $ax + by$

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?

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.

So $gcd(a,b)$ is the smallest positive number of the form $ax + by$.

Bezout's Identity

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.

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$

In addition,
  • the greatest common divisor d is the smallest positive integer that can be written as $ax + by$
  • every integer of the form ax + by is a multiple of the greatest common divisor d. 
-Wiki
Bezout's Lemma simply states that $ax + by = gcd(a,b)$ exists. We need to use Extended Euclidean Algorithm to find Bezout's Coefficients.

Coding Pitfalls

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

Reference


Saturday, July 25, 2015

Introduction to Modular Arithmetic

"In mathematics, modular arithmetic is a system of arithmetic for integers, where numbers "wrap around" upon reaching a certain value—the modulus." - Wiki
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.

General Concepts 

In order to understand modular arithmetic, we need to know the following:

Modulo Operation $\%$

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:

$5 \: \% \:  3 = 2$
$15  \: \% \:   5 = 0$
$24  \: \% \:  30 = 24$
$30  \: \% \:   24 = 6$
$4  \: \% \:   4 = 0$

Congruence Relation $\equiv$

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,

$5 \: \equiv \: 10 \text{ (mod 5) }$
$2 \: \equiv \: 7 \text{ (mod 5) }$
$4 \: \equiv \: 10 \text{ (mod 3) }$

Range

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$:

$
\begin{matrix}
\text{Normal Arithmetic } & -3 &  -2 &  -1 &  0 &  1 & 2 & 3\\
\text{Modular Arithmetic} & 0 &  1 & 2 & 0 & 1 & 2 &  0
\end{matrix}
$

Properties of Modular Arithmetic

There are four properties of Modular Arithmetic that we need to learn. 

Addition 

$$ \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)}$$

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

For example, $ ( 13 + 24 + 44 ) \: \% \: 3 = ( 1 + 0 + 2 ) \: \% \: 3 = 3 \: \% \: 3 = 0$. 

Subtraction

$$ \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)}$$

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

For example, $ ( -14 - 24 - 44 ) \: \% \: 3 = ( -2 - 0 - 2 ) \: \% \: 3 = -4 \: \% \: 3 = -1   \: \% \: 3 = 2$.

Multiplication

$$ \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)}$$

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

For example, $ ( 14 \times 25 \times 44 ) \: \% \: 3 = ( 2 \times 1 \times 2 ) \: \% \: 3 = 4 \: \% \: 3 = 1   \: \% \: 3 = 1$.

Division

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.

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.

Coding Tips


Handling Negative Numbers

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

Therefore, in C++ it is better to keep a check for negative number when doing modulo operation.
res = a % m;
if ( res < 0 ) res += m;

Modulo Operation is Expensive

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.

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:
res = a + b;
if ( res >= m ) res -= m;
This technique will not work when we are multiplying two values.

Resource


Monday, July 20, 2015

Divisor Summatory Function

Problem

Given an integer $N$, find the Sum of Number of Divisors from $1\: to \: N$. That is, find $\sum_{i=1}^{N}NOD(i)$.

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

Divisor Summatory Function

"In number theory, the divisor summatory function is a function that is a sum over the divisor function." - Wiki
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$.

Brute Force Solution - $O(N^2)$ Approach

A simple solution is to run two nested loop which counts all divisors from $1 \: to \: N$. It is simple but inefficient.
int SNOD( int n ) {
    int res = 0;
    for ( int i = 1; i <= n; i++ ) { //For each number from 1 - N
        for ( int j = 1; j <= i; j++ ) { //Find possible divisors of "i"
            if ( i % j == 0 ) res++;
        }
    }
    return res;
}

Using $NOD()$ - $O( N \times \frac{\sqrt{N}}{ln(\sqrt{N})} ) $ Approach

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.
int SNOD( int n ) {
    int res = 0;
    for ( int i = 1; i <= n; i++ ) {
        res += NOD(i);
    }
    return res;
}
Since $NOD()$ has same complexity as $factorize()$, the complexity for this code is $O( N \times \frac{\sqrt{N}}{ln(\sqrt{N})} ) $.

Using Divisor List - $O(N)$ Approach

Suppose we are trying to find $SNOD(10)$. Let us write down the divisors for each integer between $1$ and $10$.

$NOD(1) = 1, \: \{1\}$
$NOD(2) = 2, \: \{1 , 2\}$
$NOD(3) = 2, \: \{1 , 3\}$
$NOD(4) = 3, \: \{1, 2, 4\}$
$NOD(5) = 2, \: \{1, 5\}$
$NOD(6) = 4, \: \{1, 2, 3, 6\}$
$NOD(7) = 2, \: \{1, 7\}$
$NOD(8) = 4, \: \{1, 2, 4, 8\}$
$NOD(9) = 3, \: \{1, 3, 9\}$
$NOD(10) = 4, \: \{1, 2, 5, 10\}$

So $SNOD(10) = 1 + 2 + 2 + 3 + 2 + 4 + 2 + 4 + 3 + 4 = 27$.

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.

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.

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.

By repeating this from $1$ to $N$ we can find $SNOD(N)$.
int SNOD( int n ) {
    int res = 0;
    for ( int i = 1; i <= n; i++ ) {
        res += n / i;
    }
    return res;
}
We removed $NOD()$ from the code in line $4$ and replaced it with $\frac{n}{i}$.

Using Divisor Pairs - $O(\sqrt{N})$ Approach

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

$NOD(1) = 1, \: (1,1)$
$NOD(2) = 2, \: (1,2) (2,1)$
$NOD(3) = 2, \: (1,3)(3,1)$
$NOD(4) = 3, \: (1,4)(2,2)(4,1)$
$NOD(5) = 2, \: (1,5)(5,1)$
$NOD(6) = 4, \: (1,6)(2,3)(3,2)(6,1)$
$NOD(7) = 2, \: (1,7)(7,1)$
$NOD(8) = 4, \: (1,8)(2,4)(4,2)(8,1)$
$NOD(9) = 3, \: (1,9)(3,3)(9,1)$
$NOD(10) = 4, \: (1,10)(2,5)(5,2)(10,1)$

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

Now, $NOD(X)$ represents number of ordered divisor pairs $(A,B)$ such that $A \times B = X$. So, $SNOD(N)$ is same as number of ordered divisor pairs such that $A \times B \leq N$.

How to find Number of Ordered Divisor Pairs $\leq N$?

We can find this by following three steps.

1. Find number of divisor pairs $(A, B)$ where $A < B$

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

Let $u = \lfloor \sqrt{N} \rfloor$. Then for each value of $A$ between $1$ and $u$, what are the possible values of $B$?

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

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

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\}$.

When $A=2$, legal values are $\lfloor \frac{10}{2} \rfloor - 2 = 3$. $B$ can be $\{3,4,5\}$.
When $A=3$, legal values are $\lfloor \frac{10}{3} \rfloor - 3 = 0$.

$A$ cannot be bigger than $u$. So number of $(A, B)$ pairs such that $A \times B \leq 10$ and $A < B$ is 12.

2. Multiply the result with $2$

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.

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

3. Add the number of pairs $(A, B)$ where $A = B$

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

Code for $O(\sqrt{N})$ Approach

int SNOD( int n ) {
    int res = 0;
    int u = sqrt(n);
    for ( int i = 1; i <= u; i++ ) {
        res += ( n / i ) - i; //Step 1
    }
    res *= 2; //Step 2
    res += u; //Step 3
    return res;
}

Reference 

  1. forthright48 - Number of Divisors of an Integer - http://forthright48.blogspot.com/2015/07/number-of-divisors-of-integer.html
  2. Wiki - Divisor Summatory Function - https://en.wikipedia.org/wiki/Divisor_summatory_function

Sum of Divisors of an Integer

Problem

Given an integer N, find the sum of all its divisors.

For example, find Sum Of Divisors (SOD) of $12$. $12$ has the following divisors, $\{1, 2, 3, 4, 6, 12\}$. So $ SOD (12 ) = 1 + 2 + 3 + 4 + 6 + 12 = 28 $.

Sum of Divisor Function: $SOD(N)$ or $\sigma_{1}(N)$

The Sum of Divisors of $N$ is also called $\sigma_{1}(N)$. That's just a fancy way of saying Sum of Divisors of $N$. You can read more on it from Wiki. I will refer to Sum of Divisors of $N$ as either $\sigma_{1}(N)$ or $SOD(N)$. They both mean same.

Inefficient Solutions

Let us first look into some simple, but inefficient, solutions. These solutions are similar to the naive solutions we came up for "Number of Divisors of an Integer".

Loop Till $N$

This is the simplest solution. We loop from $1 \: to \: N$ and add all numbers that divide $N$.

Loop Till $\sqrt{N}$

If $N = A \times B$, where $A \leq B$, $A$ must be $\leq \sqrt{N}$. Using this fact we can run a loop from $1 \: to \: \sqrt{N}$ and add all numbers that divide $N$ and their complements to result.
int SOD ( int n ) {
    int sqrtn = sqrt ( n );
    int res = 0;
    for ( int i = 1; i < sqrtn; i++ ) {
        if ( n % i == 0 ) {
            res += i; //"i" is a divisor
            res += n / i; //"n/i" is also a divisor
        }
    }
    if ( n % sqrtn == 0 ) {
        if ( sqrtn * sqrtn == n ) res += sqrtn; //Perfect Square.
        else {
            res += sqrtn; //Two different divisor
            res += n / sqrtn;
        }
    }
    return res;
}
At line $10$, we handle $sqrtn$ separately cause when $N$ is a perfect square we don't get different values for $\{A, B\}$ pair. For example, for $49$ we have a divisor pair $7\times 7$. We need to add $7$ only once in our result.

Using Prime Factorization

It is possible to find $SOD(N)$ using the prime factorization of $N$. Let us see an example for it.

$Let \: N=12,\\ SOD(12)= 1 + 2  + 3 + 4 + 6 + 12,\\ SOD(12) = (2^0\times3^0) + (2 ^1\times 3^0) + (2^0 \times 3^1) + (2^2 \times 3^0) + (2^1 \times 3^1) + (2^2 \times 3^1) \\ SOD(12) = 2^0 ( 3^0 + 3^1 ) + 2^1 ( 3^0 + 3^1 ) + 2^2 ( 3^0 + 3^1 ),\\ SOD(12) = ( 2^0 + 2^1 + 2^2 ) \times ( 3^0 + 3^1)$.

This pattern emerges with any value of N. More generally, if $N= p_1^{a_1} \times p_2^{a_2} \times ... p_k^{a_k}$, then $$ SOD(N) = (p_1^0 + p_1^1 + p_1^2 ... p_1^{a_1} ) \times (p_2^0 + p_2^1 + p_2^2 ... p_2^{a_2} ) \times ... (p_k^0 + p_k^1 + p_k^2 ... p_k^{a_k} ) $$Using this formula and our code for Prime Factorization, we can write the function $SOD(N)$. This code is similar to $factorize()$. We only need to modify it in few places.
int SOD( int n ) {
    int res = 1;
    int sqrtn = sqrt ( n );
    for ( int i = 0; i < prime.size() && prime[i] <= sqrtn; i++ ) {
        if ( n % prime[i] == 0 ) {
            int tempSum = 1; //Contains value of (p^0+p^1+...p^a)
            int p = 1;
            while ( n % prime[i] == 0 ) {
                n /= prime[i];
                p *= prime[i];
                tempSum += p;
            }
            sqrtn = sqrt ( n );
            res *= tempSum;
        }
    }
    if ( n != 1 ) {
        res *= ( n + 1 ); //Need to multiply (p^0+p^1)
    }
    return res;
}
For each prime we find in factorization, we need to find the corresponding value of $(p^0+p^1+...p^a)$. For that, we use $tempSum$ in line $6$. It contains $p^0$ initially. Then for each successful division at line $8$, we increase the power of $p$ in line $10$ and add it to $tempSum$ in line $11$. Eventually, we multiply the final sum to result at line $14$.

In line $17$, we check if we are left with a sole prime remainder. In this case, we know that $n$ is a prime so we simply multiply $p^0 + p^1 = 1 + n$ to the result.

Reference

Thursday, July 16, 2015

Prime Number Theorem

In number theory, the prime number theorem (PNT) describes the asymptotic distribution of the prime numbers among the positive integers. It formalizes the intuitive idea that primes become less common as they become larger by precisely quantifying the rate at which this occurs. - Wiki
The definition above is hard to understand. To put it simply, PNT provides us with an estimation for Prime-Counting Function.

What is Prime-Counting Function?

This brings out another Wiki definition:
In mathematics, the prime-counting function is the function counting the number of prime numbers less than or equal to some real number $x$. It is denoted by $\pi (x)$ (this does not refer to the number $\pi$). - Wiki
This definition is easy to understand. Prime-Counting Function, $\pi (N)$, gives us number of primes less than or equal to $N$. For example, $\pi (10) = 4$ since there are $4$ primes, $\{ 2 , 3 ,  5 , 7 \}$, which are $\leq 10$.

PNT - Estimation for $\pi (N)$

So, as I was saying before, PNT provides us with an estimation for Prime-Counting Function. It states that: 
$$ \pi (N) \approx \frac {N}{ln(N)}$$
The accuracy of the estimation increase as $N$ becomes larger.

Application

We can use the theorem for complexity analysis.

Reference

Tuesday, July 14, 2015

Upper Bound for Number Of Divisors

Given a number N, we can find its number of divisors using prime factorization. But when the number gets too big, it gets difficult to factorize thus harder to find the number of divisors. So given a number N, can we estimate at most how many divisors N can have? That is, we want to find the upper bound for NOD.

NOD Upper Bounds

We proceed from loose to tighter bounds.

$NOD(N) \leq N$

A number which is greater than $N$ cannot divide $N$. For $N>2$, we can further say $NOD(N) < N$, but the improvement is negligible.

Can we make the limit tighter?

$NOD(N) \leq \frac{N}{2} + 1$

A divisor $D$, such that $\frac{N}{2} < D < N$  cannot divide $N$, since the result would be less than $2$ and greater than $1$, a fraction. So possible divisors are $D \leq \frac{N}{2}$ and $N$.

Better than before by half, but its possible to make it tighter.

$NOD(N) \leq 2 \times \sqrt{N}$

If we write $N = A \times B$, where $A \leq B$, then $A \leq \sqrt{N}$. Each of $\{A,B\}$ forms a divisor pair for $N$. $A$ can take any value from $1 \: to \: \sqrt{N}$, so it is possible to form only $\sqrt{N}$ pairs of divisors. Thus, there can be at most $2\times \sqrt{N}$ divisors for $N$.

$NOD(N) \approx 2 \times \sqrt[3]{N}$

I didn't know about this approximation until I read a blog post on Codeforces. From there I found out that in practice we can safely use $\sqrt[3]{N}$ as an upper bound, though I failed to understand the proof. Apparently this approximation has been tested for $N \leq 10^{18}$, which is large enough to be used in programming contests.

Using Highly Composite Numbers

Sometimes we need upper bound for small values of N. For example, in a problem you might need to find an upper bound of NOD for $N \leq 10^9$. For such small values of $N$ we could use NOD of largest Highly Composite Number (HCN), which is less than or equal to $N$, as an upper bound.

Read more about HCN here.

For programming contest, we could memorize values of HCN that comes frequently. Mainly $1344$ for $N \leq 10^9$ and $103,680$ for $N \leq 10^{18}$.

Application 

The upper bound for NOD is needed for complexity analysis.

Reference

  1. Codeforces - Upper bound for number of divisors: http://codeforces.com/blog/entry/14463
  2. forthright48 - Highly Composite Numbers: http://forthright48.blogspot.com/2015/07/highly-composite-numbers.html

Monday, July 13, 2015

Highly Composite Numbers

Definition

A Highly Composite Number (HCN) is a positive integer which has more divisors than any smaller positive integer ( Wiki ), i.e, if $NOD(N) > NOD(i)$, where $ 0 < i < N $, then $N$ is HCN.
For example, $6$ is HCN cause $NOD(6)=4$ which is bigger than $NOD\{1,2,3,4,5\} = \{1,2,2,3,2\}$.

Here are the first few HCN: $1, 2, 4, 6, 12, 24, 36, 48, 60, 120, 180, 240$ - A002182

Properties of Highly Composite Numbers

There are two properties of HCN and both of them are related to prime factorization.

Suppose we have a HCN with prime factorization $HCN=p_1^{a_1} \times p_2^{a_2} ... \times p_k^{a_k}$, then:
  1. First K Primes: The prime factorization of HCN will contain the first K consecutive primes. If it doesn't, then we can replace the $k^{th}$ prime in factorization with a smaller prime and still have same NOD. For example, $2^2 \times 5 = 20$ cannot be a HCN, since we can replace prime factor $5$ with $3$ and get $2^2 \times 3=12$ which is smaller than $20$ and has same NOD.
  2. Non-Increasing Power of Primes: Power of prime factors of HCN, i.e, $\{a_1, a_2 ... a_k\}$ will form a non-increasing sequence. Why is that? If power $a_j > a_i$, where $i < j $, then we can simply swap them to get a smaller number with same NOD. For example, $2 \times 3^2= 18$ cannot be HCN cause we can swap power of prime $2$ and $3$ to get $2^2 \times 3 = 12$ which is smaller with same NOD.
Therefore we conclude that Highly Composite Numbers are product of Primorials (Primorials are same as Factorials, but instead of natural number we multiply primes. $p_3\# = 2 \times 3 \times 5$ and $p_5\# = 2 \times 3 \times 5 \times 7 \times 11$ )

For example, $180$ is HCN and it can be decomposed into $180 = p_3\# \times p_2\# = (2 \times 3 \times 5) \times ( 2 \times 3 )$.

Generating Highly Composite Numbers

Problem: Given an integer N, find the largest HCN which is smaller than or equal to N.

This problem can be solved using backtrack. We just need to ensure that we satisfy the two properties of HCN.

Let's try to solve it for when $N=10^9$. We know that $p_{10} \# = 6469693230 $, so we will only have primes $ \{ 2 , 3 , 5 , 7 , 11 ,13 , 17 , 19 , 23 \} $ in HCN factorization.

Now, we will assign power to each of the prime in non increasing order to generate numbers and corresponding NOD. Out of all the numbers generated, we will keep the one with highest NOD and in case of tie, the one with smallest value.
//prime[] is a list of prime. 
int prime[] = {2, 3, 5, 7, 11, 13, 17, 19, 23 };

int resNum, resDiv, n;
void recur ( int pos, int limit, long long num, int div ) {
    if ( div > resDiv ) { //Get the number with highest NOD
        resNum = num;
        resDiv = div;
    }
    else if ( div == resDiv && num < resNum ) { //In case of tie, take smaller number
        resNum = num;
    }
    
    if ( pos == 9 ) return; //End of prime list
    
    long long p = prime[pos];

    for ( int i = 1; i <= limit; i++ ) {
        if ( num * p > n ) break;
        recur ( pos + 1, i, num * p, div * ( i + 1 ) );
        p *= prime[pos];
    }
}
Line $2$ contains a prime list. It contains first $9$ primes so it will work correctly for $N\leq 10^9$. Line $4$ contains global variables to store result. $resNum$ to hold required value and $resDiv$ to hold its NOD. Line $19$ checks if multiplying $prime[pos] ^ i$ with $res$ is becoming bigger than $N$ or not.

In line $20$, we call $recur()$ with $pos+1$ as we want to work with the next prime in list, $i$ as the new limit since we don't want the next prime to have power greater than the current one. The next two parameters are adjusted accordingly to maintain value and NOD.

In order to solve the above problem for $N=10^9$, we have to initiate and call $recur()$ from $main()$ in following manner.
int main () {
    n = 1000000000;
    resNum = 0;
    resDiv = 0;
    recur ( 0, 30, 1, 1 );
    printf ( "%d %d\n", resNum, resDiv );
}
In line $5$, we call $recur(0,40,1,1)$. $pos$ is assigned $0$ since we want to start with the first prime in prime list. $limit$ parameter is set to $30$ since $2^{30} > N$.

Reference

Saturday, July 11, 2015

Number of Divisors of an Integer

Problem

Given an integer $N$, find its number of divisors.

For example, $12$ has the following divisors $1,2,3,4,6$ and $12$. So its number of divisors is $6$.

Number Of Divisors Function: $NOD(N)$ or  $\sigma_{0}(N)$

Number of Divisors of N is also called $\sigma_{0}(N)$. That's just a fancy way of asking for the number of divisors of N. You can read more on it from Wiki. I will refer to Number of Divisors of $N$ as either $\sigma_0(N)$ or $NOD(N)$. They both mean same.

Brute Force Method O(N)

Anyways, what is the easiest way to find Number of Divisors (NOD)? We can try to divide $N$ with all numbers from $1-N$ and keep count of how many of those numbers divide $N$. This will definitely work, but obviously we need to do better.

Optimized to O($\sqrt{N}$)

In "Primality Test - Naive Methods", we established that if we can write $N=A\times B$, then one of $\{A,B\}$ must be $<=\sqrt{N}$.

So using that same idea we can find NOD by simply looping till $\sqrt{N}$ and check if that particular number divides $N$. If it doesn't then it is not a divisor and if it does then we found a $\{A,B\}$ pair of divisors. If $A \neq B$, then we add $2$ to result, otherwise we add $1$.

Examples

Suppose $N=24$. We need to loop till $\lfloor {\sqrt{24}} \rfloor= 4$. We get the following pairs of divisors, $\{1,24\}, \{2,12\}, \{3,8\}, \{4,6\}$. So our answer is $8$.

Let's try again with $N=36$.  We need to loop till $\lfloor {\sqrt{36}} \rfloor= 6$. We get the following pairs of divisors, $\{1,36\}, \{2,18\}, \{3,12\}, \{4,9\}, \{6,6\}$. So our answer is $9$. Notice that the last pair does not satisfy the condition $A \neq B$. So we add one for that pair.

This brings out an interesting observation. NOD of $N$ is always even except for when $N$ is a perfect square. Because whenever $N$ is perfect square, $\sqrt{N}$ would be its divisor and it will form a pair with itself.

Code Using Loop till $\sqrt{N}$

int NOD ( int n ) {
    int res = 0;
    int sqrtn = sqrt ( n );
    
    for ( int i = 1; i < sqrtn; i++ ) {
        if ( n % i == 0 ) res += 2; //Found a divisor pair {A,B}
    }
    
    //Need to check if sqrtn divides n
    if ( n % sqrtn == 0 ) {
        if ( sqrtn * sqrtn == n ) res++; //If n is perfect square
        else res += 2; //Not perfect square
    }

    return res;
}
We loop from $1$ to $\sqrt{N}-1$ at line 5. Then at line 10, we handle $\sqrt{N}$ seperately.

Using Prime Factorization

It is possible to find NOD from prime factorization of N.

Suppose $N=12$. Its prime factorization is $2^2 \times 3$. Is it possible to divide $12$ with $5$? No. Cause the prime factorization of $12$ does not contain $5$. We can divide $N$ with primes that appear in factorization of $N$.

Next, can we divide $12$ with $2^3$? No. The power of $2$ in prime factorization of $12$ is $2^2$. So we cannot divide $12$ with $2^3$, since $2^2$ is not divisible by $2^3$.

So, if we are trying to divide $N$, that has prime factorization $N=p_1^{a_1} \times p_2^{a_2} \times ...p_x^{a_x} $, then the divisors of $N$, $D$, will have prime factorization of form $D=p_1^{b_1} \times p_2^{b_2} \times ...p_x^{b_x}$, where $b_i <= a_i$.

For example, divisor of $12=2^2 \times 3$ will have prime factorization, $D=2^x \times 3^y$, where $x <= 2, y <= 1$. When $x = 1$ and $y=1$, $D=2^1 \times 3^1=6$ which is a divisor of $12$. If we use a different combination of $[x,y]$, then we get a different divisor. So how many different combination can we use here? We know that, $0<=x<=2$ and $0<=y<=1$. So for $x$ we can use $\{0,1,2\}$ and for y we can use $\{0,1\}$. So we can use $3 \times 2=6$ combinations of $\{x,y\}$. And that is our NOD.

So, $if \: N=p_1^{a_1} \times p_2^{a_2} \times ...p_x^{a_x}, \: then \: D=p_1^{b_1} \times p_2^{b_2} \times ...p_x^{b_x}$. We get new divisors for using different combination for $\{b_1,b_2...b_x\}$. How many different combinations can we make? $(a_1+1) \times (a_2+1) \times...(a_x+1)$. Therefore, $\sigma_0(N) = (a_1+1) \times (a_2+1) \times...(a_x+1)$.

So basically, in order to find NOD, all we need to do is factorize $N$, then take each of its prime factors power, increase them by $1$ and finally multiply all of them. Easy right?

Examples

$N=12=2^2 \times 3$. $NOD(N)= (2+1) \times (1+1) = 6$.
$N=15=3\times 5$. $NOD(N)= (1+1) \times (1+1) = 4$.
$N=252 = 2^2 \times 3^2 \times 7$. $NOD(N)=(2+1) \times (2+1) \times (1+1) = 18$.

Try out this yourself using other small values of N. 

Code Using Prime Factorization

The code for $NOD$ is not very different from the one we used for $factorize()$ in Prime Factorization of an Integer. We just need to modify it slightly to suit our needs. 
int NOD ( int n ) {
    int sqrtn = sqrt ( n );
    int res = 1;
    for ( int i = 0; i < prime.size() && prime[i] <= sqrtn; i++ ) {
        if ( n % prime[i] == 0 ) {
            int p = 0; /*Counter for power of prime*/
            while ( n % prime[i] == 0 ) {
                n /= prime[i];
                p++;
            }
            sqrtn = sqrt ( n );
            p++;/*Increase it by one at end*/
            res *= p; /*Multiply with answer*/
        }
    }
    if ( n != 1 ) {
        res *= 2; /*Remaining prime has power p^1. So multiply with 2*/
    }
    return res;
}
I have highlighted the lines that are different from $factorize()$.

Reference

Thursday, July 9, 2015

Prime Factorization of an Integer

Problem

Given an integer N, find its prime factorization.

For example, $12 = 2 \times 2 \times  3 = 2^2 \times 3$.

It is somewhat difficult to factorize an integer. If $N$ is small enough, then it is possible to factorize it using prime generate by Sieve of Eratosthenes.

Trial Division Method

Suppose we are trying to factorize $980$. How would we attempt to factorize it with pen and paper? We will try to extract the smallest possible prime factors from it and make it smaller. That is we will try to divide the number with prime numbers and keep track which are capable of dividing it.

The first prime we have is $2$. Can we divide $980$ with $2$? Yes. $980=2 \times 490$. Now we need to factorize $490$ which is smaller and thus easier.

Again, let's try to factorize $490$ with $2$. It is possible, so $980=2\times 2\times 245= 2^2 \times 245$. It is not possible to extract another $2$ from $245$ so we move on to next prime.

$3$ fails to divide $245$, so we move on to $5$. By repeating this process, we find out that $980=2^2 \times 5 \times 7^2$.

This process is kind of laborious, but we have computers to do our bidding. So for this trial division to work, all we need is a list of primes so that we can try dividing N one by one from that list. How will we get that list? Using Sieve of Eratosthenes.

But before we start generating prime, we need to answer another question. Sieve of Eratosthenes can generate prime up to a certain number. So up to which value should we generate prime numbers to factorize N?

Limit for Sieve of Eratosthenes

Well, if we are trying to factorize $N$, there is no point in generating primes that are larger than $N$. So we can generate up to $N$ and factorize using the list with primes less than or equal to $N$. But we can do even better by observing a special property of prime factors of $N$.

$N$ can have many prime factors. Some of the prime numbers will be $< \sqrt{N}$ and some $>=\sqrt{N}$. Exactly how many primes factors can $N$ have that are greater than $\sqrt{N}$?

Notice that if $N$ can never have two prime factors $>\sqrt{N}$ since, $if \: A > \sqrt{N}, then \:  A \times A > N$.

So if we generate primes up to $\sqrt{N}$ and use that list to factorize $N$ then at the end, we will either have 1 or a number greater than $\sqrt{N}$ which will be a prime for sure.

For example, $N=42$ and $\sqrt{N}=6.4 \approx 6$. So let's try to factorize using primes less than or equal to $6$ only, that is using only $[2,3,5]$. What do we get? $N = 42 = 2 \times 3 \times 7$. Since $7$ can no longer be divided with any prime from our list $[2,3,5]$ we stop. $7$ is our last missing prime factor.

So, in order to factorize $N$, we need to generate primes up to $\sqrt{N}$ only. 

Code for Sieve of Eratosthenes

vector<int>factors;
void factorize( int n ) {
    int sqrtn = sqrt ( n );
    for ( int i = 0; i < prime.size() && prime[i] <= sqrtn; i++ ) {
        if ( n % prime[i] == 0 ) {
            while ( n % prime[i] == 0 ) {
                n /= prime[i];
                factors.push_back(prime[i]);
            }
            sqrtn = sqrt ( n );
        }
    }
    if ( n != 1 ) {
        factors.push_back(n);
    }
}
Let me explain few things about the code. In line 4, inside the for loop I wrote the following condition, "i < prime.size() && prime[i] <= sqrtn". The first condition ensures that we don't go array over-bound. In order to understand it, imagine what will happen if we try to factorize a prime number $>\sqrt{N}$. Without that condition, you might get RTE.

The second condition ensures we are always trying to factorize using prime less than $\sqrt{N}$.

In line 5, we check if it is possible to divide current $N$ with prime[i]. If so we start a loop in line 6 and keep on extracting until we cannot anymore.

In line 10, we are optimizing our code. Since, after extraction of prime factors, $N$ has become smaller, we now need to only factorize with a smaller amount of prime numbers. Due to this single line, the code converges really quickly.

In line 13, we are checking if the residual number after factorization is 1 or not. If not, it is the only prime factor which is greater than $\sqrt{N}$. So we add it to our prime factor list.

Time Complexity

There are two loops in the code for $factorize()$. The first one is at line $4$. This loop runs once for every prime $\leq \sqrt{N}$. From "Prime Number Theorem", $\pi (N) \approx \frac{N}{ln(N)}$. So there are approximately $\frac{ \sqrt{N} }{ ln ( \sqrt{N} ) }$ primes $\leq \sqrt{N}$.

The second loop at line $6$ runs once for every prime factor $N$ has. So it cannot run more than $log_2(N)$ times.

So the complexity is $O(\frac{ \sqrt{N} }{ ln ( \sqrt{N} ) } + log_2(N) \: )$. In practice, if $N$ is not a prime then due to line $10$, $N$ becomes small quickly. So it is much faster than what we calculated.

Optional Optimization

There is an optional optimization we can perform only when it is possible to generate sieve up to N. Just insert a single statement in line 5.
for ( int i = 0; i < prime.size() && prime[i] <= sqrtn; i++ ) {
        if ( sieve[n] == 0 ) break; /*Checks if n is prime or not*/
        if ( n % prime[i] == 0 ) {
The new line added simply checks if the number we are trying to factorize is a prime or not. If it is prime, then obviously we cannot factorize it anymore. This optimization is highly useful when we already have sieve generated up to $N$.

Reference

  1. forthright48 - Prime Number Theorem - http://forthright48.blogspot.com/2015/07/prime-number-theorem.html

Wednesday, July 8, 2015

Sieve of Eratosthenes - Generating Primes

Problem

Given an integer N, generate all primes less than or equal to N.

Sieve of Eratosthenes - Explanation

Sieve of Eratosthenes is an algorithm that generates all prime up to N. Read this article written by +Jane Alam Jan on Generating Primes - LightOJ Tutorial. The pdf contains almost all the optimizations of Sieve of Eratosthenes.

Code

vector<int> prime; /*Stores generated primes*/
char sieve[SIZE]; /*0 means prime*/
void primeSieve ( int n ) {
    sieve[0] = sieve[1] = 1; /*0 and 1 are not prime*/

    prime.push_back(2); /*Only Even Prime*/
    for ( int i = 4; i <= n; i += 2 ) sieve[i] = 1; /*Remove multiples of 2*/

    int sqrtn = sqrt ( n );
    for ( int i = 3; i <= sqrtn; i += 2 ) {
        if ( sieve[i] == 0 ) {
            for ( int j = i * i; j <= n; j += 2 * i ) sieve[j] = 1;
        }
    }

    for ( int i = 3; i <= n; i += 2 ) if ( sieve[i] == 0 ) prime.push_back(i);
}
Some lines from the code above can be omitted depending on situation. But as a whole, the above code gives us two products, a vector<int> prime which contains all generated primes and a char[] sieve that indicates whether a particular value is prime or not. We will need sieve array for optimizations in other algorithms.

Complexity

The algorithm has runtime complexity of $O(N log (logN ) )$

Primality Test - Naive Methods

Problem

Given a number N, determine if it is a prime.

What is a Prime Number?

A prime number is a positive number greater than 1 which has no positive divisor except for 1 and itself. For example, 2, 3, 5, 11 are prime numbers. Whereas 4, 6 are non-prime or composite numbers.

O(N) Solution

From the definition, we can easily construct a function that can determine if a given number N is prime or not. Let us call it isPrime() function.
bool isPrime ( int n ) {
    if ( n <= 1 ) return false; /*n needs to be greater than 1*/
    for ( int i = 2; i < n; i++ ) {
        /*If it is possible to divide n with a number other than 1 and n, then it is not prime*/
        if ( n % i == 0 ) return false;
    }
    return true; /*Otherwise, this is prime*/
}
The code simply iterates over all values between 2 and N-1 and checks if it can divide N. It has a complexity of O(N).
We can observe that, when $i$ is greater than $n/2$, it will never divide $n$. Hence, we can optimize it a bit.
bool isPrime ( int n ) {
    if ( n <= 1 ) return false;
    for ( int i = 2; i <= n / 2; i++ ) {
        if ( n % i == 0 ) return false;
    }
    return true;
}
This, however, does not change the complexity.

O($\sqrt{N}$) Solution

We can further optimize the Primality test from the following observation. Let us consider A and B such that $N=A \times B$. Now notice that, when $A = \sqrt{N}$, B is also $\sqrt{N}$. Otherwise, $A \times B$ would not be equal to $N$. What happens when $A > \sqrt{N}$? In order for $A \times B$ to be equal to $N$, $B$ must be $< \sqrt{N}$.

So using the following observation, we can claim that if $N$ has a divisor which is $>= \sqrt{N}$, then it also has a divisor which is $<= \sqrt{N}$.

For example, $12$.
$12 = 1 \times 12$
$12 = 2 \times 6$
$12 = 3 \times 4$
$12 = \sqrt{12} \times \sqrt{12}$

Every divisor $>= \sqrt{N}$ has a complementary divisor which is $<= \sqrt{N}$.

This means that if we fail to find any divisor of $N$ below $\sqrt{N}$ then it is safe to assume we will not find any divisor above $\sqrt{N}$.
bool isPrime ( int n ) {
    if ( n <= 1 ) return false;
    int sqrtn = sqrt(n);
    for ( int i = 2; i <= sqrtn; i++ ) {
        if ( n % i == 0 ) return false;
    }
    return true;
}
Notice that I defined a new variable "sqrtn" instead of simply writing "i <= sqrt(n)" in line 4. That's because the function sqrt(n) has a complexity of O(√N) and writing it inside for loop condition would mean that the function would get called unnecessarily every time the loop iterates.

Lowest Common Multiple of Two Number

Problem

Given two number A and B, find their lowest common multiple (LCM).

We are trying to find the LCM of A and B. What is LCM? It is the smallest positive number which is divisible by both A and B.

How do we find it?

It is based on the formula that, $LCM(A,B) \times GCD(A,B) = A \times B$. How did we get this formula? I will discuss it another day. It's not that hard to figure out though. Anyways, from that formula we can derive $LCM(A,B) = \frac{A \times B}{GCD(A,B)}$.

For example, $LCM(6,15) = \frac {(6 \times 15 )}{GCD(6,15)}=  \frac{90}{3} = 30$, $LCM(3,4)= \frac{3 \times 4 }{GCD(3,4)} = \frac{12}{1} = 12$.

Code and Pitfalls

Let us try to convert the above equation for LCM to code. If we convert exactly like the equation, code would be something like the following:
int lcm ( int a, int b ) {
    return ( a * b ) / gcd ( a, b );
}
This will work, but for some cases this code will overflow. For example, if we try to find $LCM( 2^{20}, 2^{15})$, it is obvious that the LCM is $2^{20}$. But notice what happens when we follow the algorithm. We first try to multiply $2^{20}$ with $2^{15}$, which results in $2^{35}$ which is too big to fit in an "int" variable. It overflows, and returns us wrong answer. But the LCM itself should fit in the "int" variable easily.

A better way to write the above code is using the following observation. $LCM(A,B) = \frac{A \times B}{GCD(A,B)} = \frac{A}{GCD(A,B)} \times B$. Since $GCD(A,B)$ divides both A and B, the fraction $ \frac{A}{GCD(A,B)}$ will be an integer. This alternate equation avoids overflow since instead of multiplying A and B, it first reduces A to a smaller number and multiplies the resultant number with B. So, if $LCM(A,B)$ fits in "int" variable, then we will get it without the risk of intermediate products overflowing.
int lcm ( int a, int b ) {
    return ( a / gcd ( a, b ) ) * b;
}

Related Problems

Wednesday, July 1, 2015

Euclidean Algorithm - Greatest Common Divisor

Problem

Given two number A and B, find the greatest number that divides both A and B.

What we are trying to find here is the Greatest Common Divisor(GCD) of A and B. What is GCD? Like the name suggests, it the largest number ( let that be G ), that can divide both A and B. For example, $GCD(2,8) = 2$, $GCD(3,4) = 1$, $GCD(12,15) = 3$.

How do we find it?

There are many ways. One way is to list down all the divisors of A and B and then find the largest common divisors from those two lists. This is a naive method and takes too much time.

Another approach is to use Euclidean Algorithm, that works on the principle $GCD(a,b) = GCD(b,a\%b) $. Since $GCD(b,a\%b)$ is a smaller state, it is easier to find than the original. And of course, we can apply the principle on the smaller states repeatedly until the state becomes trivial. The two trivial states for GCD are $GCD(a,a) = a$ and $GCD(a,0) = a$.

Euclidean Algorithm

Recursive Version

The recursive version of GCD is simple and small. Just 4 lines of code are enough to find GCD.
int gcd ( int a, int b ) {
    if ( b == 0 ) return a;
    return gcd ( b, a % b );
}

Iterative Version

It's possible to find GCD without using recursion. This removes the recursion overhead and makes the code faster.
int gcd ( int a, int b ) {
    while ( b ) {
        a = a % b;
        swap ( a, b );
    }
    return a;
}

Built-In Version

There is also a builtin function in C++ for finding gcd. You can simply write __gcd(a,b) to find GCD(a,b).

Proof 

The Euclidean Algorithm works on the principle $GCD(a,b) = GCD(b,a\%b)$. If we can prove this, then there will be no doubt about the algorithm.

Let $g = GCD(a,b)$ and $a = k \times b + r$, where k is a non-negative integer and r is the remainder. Since $g$ divides $a$, $g$ also divides $k \times b + r$. Since $g$ divides $b$, $g$ also divides $k \times b$. Therefore, $g$ must divide $r$ otherwise $k \times b + r$ won't be divisible.

So we proved that $g$ divides $b$ and $r$.

Now lets say we have $g' = gcd (b,r)$. Since $g'$ divides both $b$ and $r$, it will divide $k \times b + r$. Therefore, $g'$ will divide $a$.

Now, can $g$ and $g'$ be two different numbers?

We will prove this using contradiction. Let's say that $g > g'$. We know that $g$ divides both $b$ and $r$. So how can $gcd(b,r)$ be $g'$ when we have a number greater than $g'$ that divides both $b$ and $r$? So $g$ cannot be greater than $g'$.

Using the same logic, we find there is a contradiction when $g < g'$. Therefore, the only possibility left is $g = g'$.

$\therefore GCD(a,b) = GCD(b, r ) = GCD( b, a \% b )$.

Complexity

It's kind of tricky. For now, just look at this answer from StackOverflow. The complexity of Eculidean Algorithm according to the answer is $O(log_{10}A + log_{10}B)$.

Apparently, there is a whole chapter in Knuth's book TAOCP. I will write a separate post on the complexity of this algorithm when I understand it. For now, just know that it works fast.

Properties

  1. Every common divisor of $a$ and $b$ is a divisor of $gcd(a, b)$.
  2. The $gcd$ is a commutative function: $gcd(a, b) = gcd(b, a)$.
  3. The $gcd$ is an associative function: $gcd(a, gcd(b, c)) = gcd(gcd(a, b), c)$.
  4. The $gcd$ of three numbers can be computed as $gcd(a, b, c) = gcd(gcd(a, b), c)$, or in some different way by applying commutativity and associativity. This can be extended to any number of numbers.
More properties can be found on Wiki page

Coding Pitfalls

Notice that the algorithm work correctly for non-negative inputs only. Try to find $GCD(4,-2)$ with the above algorithm. The correct answer should be 2. GCD will always be positive. But it returns -2. Even though the algorithm works correctly only for a non-negative number, it can easily be extended to work with negative numbers. You need to either send the absolute value of the inputs to the algorithm or use the absolute value of the return value.

Next, notice that $GCD(0,0)=0$. It should be infinity. Also, if you try to work with the return value of $GCD(0,0)$, then you might get RTE due to division by zero.

Try to be careful in the two scenarios above.

Resources

  1. Wikipedia - Greatest Common Divisor

Related Problems