CMSC 435 Algorithm Analysis & Design

Divide and Conquer

Example 1. Multiplication of polynomials of the same degree

Let P(x) and Q(x) be two polynomials of degree n - 1:

P(x) = a_n-1x^n-1 + a_n-2x^n-2 + … + a₁x + a₀ and

Q(x) = b_n-1x^n-1 + b_n-2x^n-2 + … + b₁x + b₀

Let d = én/2ù

We can write the two polynomials as

P(x) = x^dP₂(x) + P₁(x)

Q(x) = x^dQ₂(x) + Q₁(x)

where

P₂(x) = a_n-1x^n-d-1 + a_n-2x^n-d-2 + … + a_d+1x + a_d and

P₁(x) = a_d-1x^d-1 + … + a₁x + a₀

and similarly for Q₁ and Q₂

a straightforward application of the distributive law of multiplication yields:

P(x)Q(x) = x^2dP₂(x)Q₂(x) + x^d(P₂(x)Q₁(x)) + P₁(x)Q₂(x)) + P₁(x)Q₁(x)

There are four products: P₂Q₂, P₂Q₁, P₁Q₂, and P₁Q₁

each of these products is an n/2 polynomial x n/2 polynomial = n²/4 products of the coefficients. The total number of multiplications is thus 4(n²/4) + Q(n).

But this is the same number of multiplications that was needed to multiply the two polynomials without dividing them into two parts each. We have divided, but we have not conquered!

However, there is a less intuitive, but more insightful way in which we can combine the split polynomials that uses only three multiplications instead of four.

P(x)Q(x) = x^2dP₂(x)Q₂(x) + x^d(P₁(x) + P₂(x))(Q₁(x) + Q₂(x))

- P₁(x)Q₁(x) - P₂(x)Q₂(x)) + P₁(x)Q₁(x)

Only three products of polynomials are needed:

P₁(x)Q₁(x), P₂(x)Q₂(x), and (P₁(x) + P₂(x))(Q₁(x) + Q₂(x))

We can then write the (pseudocode) divide and conquer algorithm as

function polyMult(P(x), Q(x), n) {

if (n = 1)

return a₀ * b₀;

else {

d = n/2 + n%2;

split(P(x), P₂(x), P₁(x)); //split P into two halves-- P₁ and P₂

split(Q(x), Q₂(x), Q₁(x)); //split Q into two halves

R(x) = polyMult(P₂(x), Q₂(x), d);

S(x) = polyMult(P₁(x) + P₂(x), Q₁(x) + Q₂(x), d);

T(x) = polyMult(P₁(x), Q₁(x), d);

return x^2d * R(x) + x^d (S(x) - R(x) - T(x)) + T(x);

}

}

Let T(n) = the number of coefficient multiplications and additions. Then:

T(n) = 3T(n/2) + cn² if n > 1 and T(n) = 1 if n = 1

from the Master Theorem we obtain: T(n) = Q(n^lg(3))

To illustrate this algorithm, consider the following problem:

P(x) = x³ + 2x² + 3x + 4

Q(x) = 5x³ - 4x² - 2x + 1

Multiplying these two polynomials in the usual way we obtain:

P(x)Q(x) = 5x⁶ + 6x⁵ + 5x⁴ + 5x³ - 20x² - 5x + 4

Now let us apply the algorithm to the same problem.

product = polyMult([1,2,3,4],[5,-4,-2,1],4)

This first call to polyMult leads to the following sequence of recursive calls:

R₁(x) = polyMult([1,2],[5,-4],2)

S₁(x) = polyMult([1+3, 2+4], [5-2, -4+1],2) = polyMult([4,6], [3, -3],2)

T₁(x) = polyMult([3,4],[-2,1],2)

Each of these calls to polyMult produces the following:

To find R₁(x), where the subscript indicates this is the value of R(x) for the original call to polyMult, we have:

R₂(x) = polyMult([1],[5],1) = 5 // from the if n == 1 block

S₂(x) = polyMult([1+2],[5-4],1) = 3

T₂(x) = polyMult([2],[-4],1) = -8

R₁ß return R₂(x)x² + (S₂(x) - R₂(x) - T₂(x))x + T₂(x) = 5x² + 6x - 8

Similarly, to find S₁(x) we have polyMult([4,6],[3,-3],2)

R₂(x) = polyMult([4],[3],1) = 12 // from the if n == 1 block

S₂(x) = polyMult([4+6],[3-3],1) = 0

T₂(x) = polyMult([6],[-3],1) = -18

S₁ß return R₂(x)x² + (S₂(x) - R₂(x) - T₂(x))x + T₂(x) = 12x² + 6x - 18

and, to find T₁(x) we have polyMult([3,4],[-2,1],2)

R₂(x) = polyMult([3],[-2],1) = -6 // from the if n == 1 block

S₂(x) = polyMult([3+4],[-2+1],1) = -7

T₂(x) = polyMult([4],[1],1) = 4

T₁ß return R₂(x)x² + (S₂(x) - R₂(x) - T₂(x))x + T₂(x) = -6x² - 5x + 4

Now given R₁, S₁, and T₁, we are ready to return the product.

PQ = R₁(x)x⁴ + (S₁(x) - R₁(x) - T₁(x))x² + T₁(x)

= (5x² +6x-8)x⁴ + (12x²+6x-18)x² -(5x² +6x-8)x² -(-6x²-5x+4)x² - 6x² - 5x + 4

= 5x⁶ +6x⁵ + (-8 + 12 - 5 +6)x⁴ + (6 - 6 +5)x³ + (-18 + 8 - 4 -6)x² - 5x + 4

= 5x⁶ + 6x⁵ + 5x⁴ + 5x³ - 20x² - 5x + 4

Here is an implementation of the previous pseudocode algorithm. The polynomials are represented by vectors of the coefficients. The coefficients of the higher power terms are in the lower indices of the array. The powers of x are implicit.

public double [ ] polyMult(double [ ] p, double [ ] q, int n) {

if (n == 1) {

double [ ] temp = {p[0] * q[0]};

return temp;

}

int d = n/2 + n%2;

//split p and q

double [ ] p2 = new double[d];

double [ ] q2 = new double [d];

for (int i = 0; i < d; i++) {

p2[i] = p[i];

q2[i] = q[i];

}

double [ ] p1 = new double[d-n%2];

double [ ] q1 = new double[d-n%2];

for (int i = 0; i < d; i++) {

p1[i] = p[i+d];

q1[i] = q[i+d];

}

double [ ] r = polyMult(p2,q2,d);

double [ ] s = polyMult(p2+p1, q2+q1,d);

double [ ] t = polyMult(p1,q1,d);

double [ ] prod = new double[2 * n];

//remember, the lowest index holds the coefficient of the highest power term

//multiplying by x^d involves shifting d places to the left which is equivalent

//to shifting the terms multiplied by a lower power d places to the right

for (int i = 0; i < d; i++) {

prod[i] = r[i];

prod[i+d] += s[i] - r[i] - t[i];

prod[i+2*d] += t[i];

}

return prod;

}

Suppose Q(x) is of degree n-1 and P(x) of degree m-1 where m < n.

Assume for convenience that n = km for some positive integer k.

Q(x) = Q_k(x)x^m(k-1) + Q_k-1(x)x^m(k-2) + … + Q₂(x) + Q₁(x)

It follows then that

P(x)Q(x) = P(x) Q_k(x)x^m(k-1) + P(x) Q_k-1(x)x^m(k-2) + … + P(x)Q₁(x)

function polyMultUnequal(P(x), Q(x), n, m) {

prodPoly(x) = 0;

for (int i = 1; i <= k; i++)

Q_i(x) = b_i(m-1)x^m-1 + b_i(m-1)x^m-2 + … + b_(i-1)m;

for (int i = 1; i <= k; i++)

prodPoly(x) += x^(i-1)m polyMult(P(x), Q_i(x), m);

return prodPoly(x);

}

The complexity of this Algorithm is Q(km^lg3) = Q((n/m) m^lg3) = Q(n m^lg(3/2))

Example 2 -- Strassen's Algorithm for matrix Multiplication

Let A and B be two NxN matrices.

Divide A and B into 4 (N/2 x N/2) submatrices

|A₁₁ A₁₂| |B₁₁ B₁₂|

A = | | B = | |

|A₂₁ A₂₂| |B₂₁ B₂₂|

The usual matrix multiplication requires 8 products of the submatrices:

A₁₁B₁₁ A₁₂B₂₁ A₁₁B₁₂ A₁₂B₂₂

A₂₁B₁₁ A₂₂B₂₁ A₂₁B₁₂ A₂₂B₂₂

Strassen found a clever, non-intuitive way to carry out matrix multiplication using only 7 matrix multiplications.

M₁ = (A₁₁ + A₂₂)(B₁₁ + B₂₂)

M₂ = (A₂₁ + A₂₂)B₁₁

M₃ = A₁₁(B₁₂ - B₂₂)

M₄ = A₂₂(B₂₁ - B₁₁)

M₅ = (A₁₁ + A₁₂)B₂₂

M₆ = (A₂₁ - A₁₁)(B₁₁ + B₁₂)

M₇ = (A₁₂ - A₂₂)(B₂₁ + B₂₂)

and

|M₁ + M₄ - M₅ + M₇ M₃ + M₅ |

AB = | |

|M₂ + M₄ M₁ + M₃ - M₂ + M₆ |

Let T(n) = the number of multiplications for the product of two n by n matrices.

Then

1 if n = 1

T(n) =

7T(n/2) + 18n² if n > 1

From the Master Theorem, T(n) = Q(n^lg(7)) = Q(n^2.81)

Example 3. Tiling With Triominoes

A tromino is an object composed of three 1 x 1 squares. Trominoes can appear in four different orientations:

These orientations are labeled UL -- upper left; UR -- upper right; LL -- lower left; and LR -- lower right.

The problem we are to solve is to tile with trominoes an N x N board with one deficiency, where N is 2^m for m > 1. A deficiency is where one of the unit squares in the board has been "removed".

To show that such a tiling is feasible, we first need to prove that the number of non-deficient unit squares is divisible by three. (Since a tromino covers a set of three squares.)

Exercise: Prove by induction that the number of non-deficient squares is divisible by three.

Hint: The number of non-deficient square = N² - 1 = 2^2m -1 = 4^m - 1

This can be written as 4(4^m-1) - 4 + 3 = 4(4^m-1 - 1) + 3

This last term is divisible if (4^m-1 - 1) is divisible. Now develop the inductive

proof.

We will employ a divide and conquer strategy to solve this problem.

Divide the board into four quadrants. The deficiency will be located in one of the quadrants.

Our approach is to divide the board into four quadrants, locate the quadrant with the deficiency, place a tromino with the appropriate orientation at the corner of the quadrant with the deficiency, mark the three squares covered by the tromino in the three remaining quadrants as deficient, and recursively tile the four quadrants.

In the following algorithm we will locate the board to be tiled and the deficient square by the coordinates of their lower left corner.

public void tile(int n, Pt pb, Pt, pd) {

if (n == 2) {

String str; //identify the orientation of the tromino

if (pb.equals(pd) )

str = "UR";

else if (pb.getY( ) == pd.getY( ))

str = "UL";

else if (pb.getX( ) == pd.getX( ))

str = "LR";

else

str = "LL";

place(str, pb.getX( ) + 1, pb.getY( ) + 1);

}

else {

//partition into 4 quadrants and locate the quadrant with the deficiency

Pt p = new Pt(pb.getX( ) + n/2, pb.getY( ) + n/2);

String str; //determine the orientation of the tromino

if (pd.getX( ) <= p.getX( ) && pd.getY( ) >= p.getY( )) {

//deficiency is in upper left quadrant

str = "LR";

tile(n/2, new Pt(pb.getX( ), p.getY( )), new Pt(pd)); //upper left

tile(n/2, new Pt(p), new Pt(p)); //upper right

tile(n/2, new Pt(pb), new Pt(p.getX( ) - 1, p.getY( ) -1)); //lower left

tile(n/2, new Pt(p.getX( ), pb.getY( )), new Pt(p.getX( ), p.getY( ) - 1);

}

else if (pd.getX( ) >= p.getX( ) && pd.getY( ) >= p.getY( )) {

//deficiency is in upper right quadrant

str = "LL";

//upper left

tile(n/2, new Pt(pb.getX( ), p.getY( )), new Pt(p.getX( )-1, p.getY( ));

tile(n/2, new Pt(p), new Pt(pd)); //upper right

tile(n/2, new Pt(pb), new Pt(p.getX( ) - 1, p.getY( ) -1)); //lower left

tile(n/2, new Pt(p.getX( ), pb.getY( )), new Pt(p.getX( ), p.getY( ) - 1);

}

else if (pd.getX( ) <= p.getX( ) && pd.getY( ) <= p.getY( ))

//deficiency is in lower left quadrant

str = "UR";

//upper left

tile(n/2, new Pt(pb.getX( ), p.getY( )), new Pt(p.getX( )-1, p.getY( ));

tile(n/2, new Pt(p), new Pt(p)); //upper right

tile(n/2, new Pt(pb), new Pt(pd)); //lower left

tile(n/2, new Pt(p.getX( ), pb.getY( )), new Pt(p.getX( ), p.getY( ) - 1);

}

else

//deficiency is in lower right quadrant

str = "UL";

//upper left

tile(n/2, new Pt(pb.getX( ), p.getY( )), new Pt(p.getX( )-1, p.getY( ));

tile(n/2, new Pt(p), new Pt(p)); //upper right

tile(n/2, new Pt(pb), new Pt(p.getX( ) - 1, p.getY( ) -1)); //lower left

tile(n/2, new Pt(p.getX( ), pb.getY( )), new Pt(pd);

}

place(str, p.getX( ), p.getY( ));

}

}