  function y = DiagPade(x,q)
% x is a column n-vector
% y is a column n-vector with the property that y(k) = Rqq(x(k)) where
% Rqq is the qq-Pade aproximation to exp().

% We need a recursiopn for the coefficients. 
% If c(k) = (2q-k)!q! / (2q)!k!(q-k)! then
%          c(0) = 1
%          c(k+1) = c(k)*(q-k)/(2q-k)(k+1))  k>=0

% Initialize running sums for both the numerator and denominator...
n = length(x);
Numer = ones(n,1);
Denom = ones(n,1);
% xpower will be used to house powers of x, i.e., x.^k
xpower =  ones(n,1);
% c will be used to house the coefficients c(k)
c = 1;
for k=0:q-1
    % add in the (k+1)-st terms
    xpower = xpower.*x;
    c = -c*(q-k)/((2*q-k)*(k+1));
    Numer = Numer + abs(c)*xpower;
    Denom = Denom + c*xpower;
end
y = Numer./Denom;

% Partial Credit
% Correct result: 3 points.
% Vectorization: 1 point.
% If the computing of each term in the summation is based on the previous term, 
% 1 more point. 
% Any calling of factorial, power or self-written function that makes 
% the complexity O(q^2) will make you lose this point.

