Spring 2001 CS100M Project 4 Solutions

Project 4: Part 1. Rightmost, Longest, Arithmetic Subsequence

Before we present the solution, here are some preliminaries.

This problem uses the same ideas in Project 3.1: maintain the answer so far (rightmost, longest, arithmetic subsequence so far) and maintain information about the previous input value (longest arithmetic subsequence ending on the previous input value). This is essentially the loop invariant. Note that a "loop invariant" is not Matlab code. It is the comment you write before a loop to specify what properties remain unchanged through each iteration of the loop.

As stated in the hint, you can represent an arithmetic sequence simply as the initial value, the gap, and the length.

Whenever you write code, you need think about what to do about boundary conditions, aka degenerate cases. This means you need to get as precise a specification of the problem as possible, and then to intepret the specification as generously as possible. That is, if the specification does not mention a restriction, then you may not assume that restriction.

In this case, the question is what to do about short input sequences. Since they were not mentioned in the specification, we must assume they are legal. It turns out the specification already covers this: we will want the longest arithmetic subsequence, where an arithmetic sequence is one in which all gaps between successive elements are the same. Recall from the 2/08 lecture that all can be thought of as there are no counter-examples. Thus,

An empty sequence is an arithmetic sequence.
A sequence with just one element is an arithmetic sequence.

Now we're ready to start writing code.

We start out by defining the stopping value and the prompt users see when inputting values. We can easily change the stopping value by just changing the definition of stopping.

Next, we define the variables to hold the information we maintain during the loop.

In the loop, we look at each number the user enters and calculate the gap between this number and the previous number. If the gap is consistent with the current arithmetic sequence (the sequence that ended on the previous input value), then we extend the current sequence. If this current sequence now turns out to be equal to or longer than the rightmost, longest, arithmetic subsequence encountered so far, then this current sequence is promoted to the rank of rla (right most longest arithmetic subsequence). If the gap is not consistent with the current arithmetic sequence, then we have a new (different gap) current arithmetic sequence.


stopping = -1; % stopping value
prompt = sprintf('next number (%d to stop)? ', stopping); % input prompt
num = input(prompt); % current input value
prev = num; % previous input value
% rightmost, longest arithmetic subsequence so far -- up to but excluding $num$
    longest.lo = num; % start
    longest.gap = 0;  % gap
    longest.len = 0;  % length
% "current" subsequence: longest arithmetic subsequence ending on $prev$
    cur = longest; 
% process $stopping$-terminated input sequence
% inv: maintain variable definitions above
while num ~= stopping
    % update current subsequence
        if cur.gap == num - prev
            cur.len = cur.len + 1;
        else
            cur.lo = prev;
            cur.gap = num - prev;
            cur.len = 2;
        end
    % update rightmost, longest arithmetic subsequence so far
        if cur.len >= longest.len
            longest = cur;
        end
    % read next input value
        prev = num;
        num = input(prompt);
end
% display rightmost, longest arithmetic subsequence
    fprintf('rightmost, longest arithmetic subsequence: ['); 
    fprintf(' %d', (0:longest.len-1)*longest.gap+longest.lo);
    fprintf(' ]\n');

The rightmost, longest arithmetic subsequence for user input

    17 0 -2 -4 -6 -9 -12 15 11 10 9 2 1 5 9 13 17 -1

is 1 5 9 13 17

Project 4: Part 2. Polynomials

Observe that the problem asks that our functions return results in canonical form:

To make the polynomial representation unique, we will impose the condition that the representation be in canonical form...

It is good general programming practice to make programs deal inputs that are as general as possible (here: polynomials in non-canonical form), but to provide results in form as as uniform as possible (here: polynomials in canonical form).

makepoly(p)

Except for polynomial 0, all polynomials have at least one non-zero coefficient. All we do here is find the right-most non-zero coefficient and copy everything up to and include it, effectively "cutting off" the coefficients to the right of it (we know that all these are zero). If the polynomial is in fact the constant 0, we can not apply this idea, as it would yield an empty array, which is not a legal polynomial according to our definition.

function q = makepoly(p)
    % q = makepoly(p): return canonical form $q$ of polynomial $p$

    if all(p == 0)
      q = 0;
    else
      q = p(1:max(find(p ~= 0)));
    end

Alternative solution (function body only):

    while length(p)>1 & p(end) == 0
        p(end) = [];
    end

The other functions that return a polynomial result should return it in canonical form. To avoid duplicating code, we do so by calling makepoly.

addpoly(p1,p2)

We don't have to worry whether the polynomials are in canonical form or not, we just have to append zeros to the one with lower degree, so that we can add up the corresponding coefficients as arrays.

function p = addpoly(p1, p2)
% p = addpoly(p1, p2): return canonical polynomial sum $p=p1+p2$;
%                      $p1,p2$ are not necessarily canonical

    if length(p1) > length(p2)
        p = p1 + [p2 zeros(1, length(p1)-length(p2))];
    else
        p = [p1 zeros(1, length(p2) - length(p1))] + p2;
    end
    p = makepoly(p);

Here is the same idea implemented more concisely (function body only):

    p = makepoly( [ p1 zeros(1, length(p2) - length(p1)) ] ...
                + [ p2 zeros(1, length(p1) - length(p2)) ]);

Here is another way to implement the same idea, using the property that x(a:b)=0 has no effect when a>b since in that case a:b=[] is empty:

    p1(end+1:length(p2)) = 0;
    p2(end+1:length(p1)) = 0;
    p = makepoly(p1 + p2);

Here is an alternative function body, where we add up only the coefficients in positions common to both polynomials, and we tack on the leftover coefficients at the end:

    if length(p1) > length(p2)
      p = [p1(1:length(p2))+p2  p1(length(p2)+1:end)];
    else
      p = [p2(1:length(p1))+p1  p2(length(p1)+1:end)];
    end
    p = makepoly(p);

Here is the same idea implemented more concisely, using the property that p(m+1:m)=[] since m+1:m=[]:

    m = min(length(p1), length(p2));
    p = makepoly([p1(1:m)+p2(1:m) p1(m+1:end) p2(m+1:end));

subpoly(p1,p2)

Function subpoly is almost identical to addpoly above. This suggests that we would like to somehow avoid duplicating code and/or effort, which it turns out is easy to do.

function p = subpoly(p1,p2)
% p = subpoly(p1, p2): return canonical polynomial difference $p=p2-p1$;
%                      $p1,p2$ are not necessarily canonical
    p = addpoly(p2, -p1 );

Note: In the P4.2 write-up, the table of functions requests p2-p1, but the Example show p1-p2, so we will accept either implementation.

mulpoly(p1,p2)

We can reduce multiplication to a series of polynomial additions, which we perform in a loop. To explain, let p2 be represented as [a₀, a₁, a₂,...,a_n-1,a_n] and let the original value of p1 be q. Below are two versions of the algorithm in pseudocode:


    p <-- 0
    for i <-- 0 .. n-1
        p <-- p + a_i * q * xⁱ


    p <-- 0
    for i <-- 0 .. n-1
        p <-- p + a_i * p1  
        p1 <-- p1 * x

When we translate to Matlab code, to multiply by a_i is easy enough, but how do we multiply by a power of x? The answer is to concatenate an appropriate number of 0s in front of the vector of coefficients. Here is an implementation of the pseudocode above, left:

function p = mulpoly(p1, p2)
% p = mulpoly(p1, p2): return canonical polynomial product $p=p1*p2$;
%                      $p1,p2$ are not necessarily canonical
    p = 0;
    for i = 0 : length(p2)-1
        p = addpoly(p, [zeros(1, i) p2(i+1) * p1]);
    end
    p = makepoly(p);

Observe there is a little trickiness about whether to use i or i+1 or i-1. Here is an implementation of the pseudocode above, right (function body only), which somewhat avoids this trickiness:

    p = 0;
    for a = p2
        p = addpoly(p, a * p1); 
        p1 = [0 p1];
    end
    p = makepoly(p);

Here's another possible solution, using the idea of collecting like terms a_ix^i-1 * b_jx^j-1 = (a_i*b_j) x^i-1+j-1 of the same power i-1+j-1 of x:

function p = mulpoly(p1,p2)
    p = zeros(1,length(p1)+length(p2)-1);
    for i = 1:length(p1)
        for j = 1:length(p2)
            p(i+j-1) = p(i+j-1) + p1(i) * p2(j);
        end
    end
    p = makepoly(p); % in case p1 or p2 was 0

Although we will probably not require you write recursive functions in CS100M, some of you might sometimes find them easier to understand. Using the observation that

(x * p1 + a)  times  p2   =   x * (p1 times p2)  +  a * p2,

here is a recursive solution:

function p = rmulpoly(p1,p2)
% p = rmulpoly(p1,p2): recursively compute canonical polynomial product $p=p1*p2$

    p = p1 * p2(1);
    if length(p2) > 1
        p = addpoly(p, [0 rmulpoly(p1, p2(2:end))]);
    else
        p = makepoly(p);
    end

[q, r] = divpoly(p1,p2)

While the other operations are defined on all polynomials, division is not. We have to take special precautions to handle division by polynomial p2=0. But what should the result be?

If we divided numbers, x/0 would might reasonable take three values: +inf, if x were positive, -inf, if x were negative, and nan, if x were 0, since the result of the latter operation is undefined. Note that for a specific value of x, we get back a single, specific answer.

Turning our attention to polynomials, p(x)/0 again can take three different values if we plug in a value for x: +inf, if p(x) were positive, -inf, if p(x) were negative, and nan, if p(x) were 0. However, note that to get an answer, we needed to plug in a value for x. Without a specific value for x, p(x)/0 makes no sense: none of +inf, -inf, and nan is the single, well-defined answer. Therefore, we have to return nan or call the error function to halt the program.

The division algorithm is the same one we (people) use. We determine the coefficients of the quotient from right to left (from high powers down to lower powers), one at a time. The next coefficient of the quotient is the ratio of the leading (i.e. righmost) coefficients of the two polynomials (p1 will change at every step). Once the coefficient is determined, it is used to multiply p2, which is then "shifted" by concatenating the results with an array of zeros to align corresponding coefficients (this is similar to the idea employed in the multiplication function).

One tricky point is that although the result of the subtraction should eliminate the most significant (i.e. rightmost) coefficient of p1, due to the finite precision of computations in a computer, it is not at all unlikely that a very small error will occur (say, on the order of p1(end)*10^(-15)). You can get an intuition for why this happens if you think of computing (1/3)*3 using a precision of only two decimals. You would get 0.99, wouldn't you? Because we know that p1(end) should be 0, we can just delete it.

The loop is performed as long as what is left of p1 has a degree not less than that of p2. When the condition doesn't hold, the current value of p1 is the remainder of the division.

Note how the quotient array is initialized to 0, a vector containing the single element, 0. This is so that if the loop's body is not executed even once, the quotient polynomial is assigned the correct answer 0. However, note that this can yield a leading coefficient of 0, but that is taken care of by a call to makepoly. We could perhaps avoid the "extra" work of inserting then deleting a 0, but compared to the rest of the work the function does, this is not a lot of work and hence not worth worrying about.

function [q, r] = divpoly(p1, p2)
% [q, r] = divpoly(p1, p2): compute quotient $q$ and remainder $r$ of $p1/p2$
%                           $p1,p2$ are polynomials, not necessarily canonical
%                           $q,r$ are canonical polynomials

    if p2 == 0
        q = nan;
        r = nan;
        return; % exit the function
    end

    % we now know that p2 is not 0

    q = 0;

    while length(p1) >= length(p2)
        nextQ = p1(end) / p2(end);
        p1 = p1 - [zeros(1, length(p1) - length(p2)) nextQ*p2];
        % deal with roundoff errors
            p1(end) = [];
        q = [nextQ q];
    end

    q = makepoly(q);
    r = makepoly(p1);

Although we will probably not require you write recursive functions in CS100M, some of you might sometimes find them easier to understand. Since the recursive solution for multiplication was short, we might wonder if there is a short recursive solution for division. There doesn't seem to be a short, recursive solution.

Nevertheless, for completeness, we present a recursive solution. Here are the key observations:

if p1/p2	=	quotient q, remainder r	i.e. *p1 = q p2 + r**
and if *(xr+b)/p2**	=	quotient a, remainder s,	i.e. *xr + b = a * p2 + s**
then *(xp1 + b)/p2**	=	quotient *x q + a**, remainder s	i.e. *xp1 + b = (xq+a) p2 + s**

The last equality x*p1 + b = (x*q+a) * p2 + s follows from plugging in p1=q*p2+r and x*r+b=a*p2+s.

Here is a recursive solution using the ideas above:

function [q,s] = rdivpoly(p1,p2)
% [q,s] = divpoly(p1,p2): rec. compute quotient $q$ and remainder $s$ of $p1/p2$
%                         $p1,p2$ are polynomials, not necessarily canonical
%                         $q,r$ are canonical polynomials
    if length(p2)>length(p1)
        q = 0; s = makepoly(p1);
    else
        b = p1(1); p1(1) = []; % original p1 = b + x * current p1
        [q,r] = rdivpoly(p1, p2); % $p1 = q * p2 + r$, as in observations above
        xrb = [b r]; % x * r + b
        a = xrb(end) / p2(end); $a$ = quotient $xrb/p2$ in observations above
        q = [a q];
        s = addpoly(xrb, -a * p2); % $s = xrb - a*p2$ from observations above
        % deal with roundoff errors
            s(end) = []; 
    end

evalpoly(p,x0)

function v = evalpoly(p,x0)
% v = evalpoly(p,x0): evaluate polynomial $p(x)$ at $x=x0$

    x0 = x0 .^ (0:length(p)-1);
    v = sum(x0 .* p);

testpoly

We now turn our attention to testing these functions. In the testpoly script shown below we test the most relevant cases for the functions at hand. Note that we have hard-coded the expected result, to ease the task of those who run the program.

% Canonical form of 0.
    p = makepoly(0);
    disp([ '[' num2str(p) ']. Expected: [0].' ]);
% Canonical form of a non-zero constant.
    p = makepoly(17);
    disp([ '[' num2str(p) ']. Expected: [17].' ]);
% Canonical form of a polynomial *in* canonical form.
    p = makepoly([1 0 2 0 0 3]);
    disp([ '[' num2str(p) ']. Expected: [1  0  2  0  0  3].' ]);
% Canonical form of a polynomial *not* in canonical form.
    p = makepoly([1 0 2 0 0 3 0 0]);
    disp([ '[' num2str(p) ']. Expected: [1  0  2  0  0  3].' ]);
% Add two polynomials, one degenerate.
    p = addpoly([1 0 2 0 0 3], 5);
    disp([ '[' num2str(p) ']. Expected: [6  0  2  0  0  3].' ]);
% Add two non-degenerate polynomials.
    p = addpoly([1 0 2 0 0 3], [0 7 8]);
    disp([ '[' num2str(p) ']. Expected: [1  7 10  0  0  3].' ]);
% Add the same two polynomials, in reverse order.
    p = addpoly2([0 7 8], [1 0 2 0 0 3]);
    disp([ '[' num2str(p) ']. Expected: [1  7 10  0  0  3].' ]);
% Add two polynomials, one degenerate.
    p = addpoly([1 0 2 0 0 3], 5);
    disp([ '[' num2str(p) ']. Expected: [6  0  2  0  0  3].' ]);
% Add two non-degenerate polynomials.
    p = addpoly([1 0 2 0 0 3], [0 7 8]);
    disp([ '[' num2str(p) ']. Expected: [1  7 10  0  0  3].' ]);
% Add the same two polynomials, in reverse order.
    p = addpoly([0 7 8], [1 0 2 0 0 3]);
    disp([ '[' num2str(p) ']. Expected: [1  7 10  0  0  3].' ]);
% Add two polynomials where leading coefficient cancel.
    p = addpoly([0 7 8], [2 0 -8]);
    disp([ '[' num2str(p) ']. Expected: [2  7].' ]);
% Subtract a degenerate polynomial (result depends on subtraction order)
    p = subpoly([1 0 2 0 0 3], 5);
    disp([ '[' num2str(p) ']. Expected: [-4  0  2  0  0  3].' ]);
% Subtract *from* a degenerate polynomial (result depends on subtraction order)
    p = subpoly(5, [1 0 2 0 0 3]);
    disp([ '[' num2str(p) ']. Expected: [4  0 -2  0  0 -3].' ]);
% Subtract two polynomials, none degenerate (result depends on subtraction order)
    p = subpoly([1 0 2 0 0 3], [0 7 8]);
    disp([ '[' num2str(p) ']. Expected: [1 -7 -6  0  0  3].' ]);
% Subtract the same polynomials, with reversed positions (result depends on subtraction order)
    p = subpoly([0 7 8], [1 0 2 0 0 3]);
    disp([ '[' num2str(p) ']. Expected: [-1  7  6  0  0 -3].' ]);
% Test that subtract undoes addition: same polynomials used for addition tests
    p = subpoly([0 7 8], addpoly([1 0 2 0 0 3], [0 7 8]);
    disp([ '[' num2str(p) ']. Expected: [1 0 2 0 0 3].' ]);

    p = subpoly([1 0 2 0 0 3], addpoly2([0 7 8], [1 0 2 0 0 3]);
    disp([ '[' num2str(p) ']. Expected: [0 7 8].' ]);

    p = subpoly([5], addpoly([1 0 2 0 0 3], [5]);
    disp([ '[' num2str(p) ']. Expected: [1 0 2 0 0 3].' ]);

    p = subpoly([0 7 8], addpoly([1 0 2 0 0 3], [0 7 8]);
    disp([ '[' num2str(p) ']. Expected: [1 0 2 0 0 3].' ]);

    p = subpoly([1 0 2 0 0 3], addpoly([0 7 8], [1 0 2 0 0 3]);
    disp([ '[' num2str(p) ']. Expected: [0 7 8].' ]);

    p = subpoly([2 0 -8], addpoly([0 7 8], [2 0 -8]);
    disp([ '[' num2str(p) ']. Expected: [0 7 8].' ]);
% Multiply by 0 - tests whether result will be in canonical form.
    p = mulpoly([1 0 2 0 0 3], 0);
    disp([ '[' num2str(p) ']. Expected: [0].' ]);
% Multiply two polynomials, none degenerate.
    p1 = [1 0 2 0 0 3]; p2 = [0 7 8];
    p = mulpoly(p1, p2);
    disp([ '[' num2str(p) ']. Expected: [' num2str(fliplr(conv(fliplr(p1),fliplr(p2)))) '].' ]);
% Multiply the same polynomials, in reverse order.
    p1 = [0 7 8]; p2 = [1 0 2 0 0 3];
    p = mulpoly(p1, p2);
    disp([ '[' num2str(p) ']. Expected: [ num2str(fliplr(conv(fliplr(p1),fliplr(p2)))) '].' ]);
% Division by 0.
    p = divpoly([1 0 2 0 0 3], 0);
    disp([ num2str(p) '. Expected: NaN.' ]);
% Division of two non-degenerate polynomials.
    p1 = [3 2 0 7]; p2 = [1 2 1];
    [q r] = divpoly(p1,p2);
    [q2 r2] = deconv(fliplr(p1), fliplr(p2));
    q2 = fliplr(q2); r2 = fliplr(r2);
    disp([ ' q == [' num2str(q) ']. Expected: [' num2str(q2) '].' ]);
    disp([ ' r == [' num2str(r) ']. Expected: [' num2str(r2) '].' ]);
% Division of the same polynomials, in reverse order.
    p1 = [1 2 1]; p2 = [3 2 0 7];
    [q r] = divpoly(p1,p2);
    [q2 r2] = deconv(fliplr(p1), fliplr(p2));
    q2 = fliplr(q2); r2 = fliplr(r2);
    [q r] = divpoly(p1, p2);
    disp([ ' q == [' num2str(q) ']. Expected: [' num2str(q2) '].' ]);
    disp([ ' r == [' num2str(r) ']. Expected: [' num2str(r2) ']' ]);
% Evaluation of a degenerate polynomial.
    p = 5;
    val = evalpoly(p, 9);
    disp(['val == ' num2str(val) '. Expected: ' polyval(fliplr(p,9)) '.']);
% Evaluation of a non-degenerate polynomial.
    p = [3 2 0 7];
    val = evalpoly(p, 9);
    disp(['val == ' num2str(val) '. Expected: ' polyval(fliplr(p,9)) '.']);

If you want to generate more test cases automatically, you can run a loop and generate polynomials with random coefficients. Here's one possible implementation:

function bad = superPolyTest% this does lots of stuff you didn't need to do
    d = 6;
    n = 5;
    bad = 0;
    poke = 1/3;
    epsilon = 1e-8;
    for i = 1:1000
        p1 = ceil(rand(1,d) * (2*n+1)) - n;
        p2 = ceil(rand(1,d) * (2*n+1)) - n;
        % poke lots of 0s into p1, p2
             p1(rand(1,d)<poke) = 0;  p2(rand(1,d)<poke) = 0;
        p1 = makepoly(p1); p2 = makepoly(p2);
        unhappy = ''; % operations that failed

        if ~eqpoly(p1, subpoly(addpoly(p1,p2),p2))
            unhappy = [unhappy '+'];
        end

        if ~eqpoly(conv(p1,p2), mulpoly(p1,p2))
            unhappy = [unhappy '*'];
        end

        if length(p2) ~= 1 & p2 ~= 0
            [q1,r1] = deconv(fliplr(p1),fliplr(p2));
            r1(abs(r1)<epsilon) = 0;
            r1 = makepoly(fliplr(r1)); q1 = fliplr(q1);
            [q2,r2] = divpoly(p1,p2);
            if ~eqpoly(q1,q2)
                unhappy = [unhappy 'q'];
            end
            if ~eqpoly(r1,r2)
                unhappy = [unhappy 'r'];
            end
        end

        if polyval(fliplr(p1),2) ~= evalpoly(p1,2)
            unhappy = [unhappy 'e'];
        end

        if length(unhappy)
            fprintf('%s failed: p1 = ', unhappy);
            fprintf('%+d ', p1);
            fprintf(' ; p2 = ');
            fprintf('%+d ', p2);
            fprintf('\n');
            bad = bad+length(unhappy);
            if bad > 20, return; end
        end
    end

function t = eqpoly(p,q)
    epsilon = 1e-8;
    if length(p) ~= length(q)
        t = 0;
    else
        t = all(abs(p - q) < epsilon);
    end