function Et = taylor_error2_cjp(x,n)

%TAYLOR_ERROR2_CJP computes the expression in Equation 1 of HW1
%   Et = TAYLOR_ERROR2_CJP(X,N) computes Et as a function of the scalar
%   input X and the non-negative integer input N according to the equation
%   given in Task2 of HW1:  Et is the true absolute error of an N-th order
%   Taylor series approximation of exp(x).  If N is not an integer, the
%   function will truncate it to the nearest integer smaller than N.
%
%   Inputs:
%      X: any scalar
%      N: the order of the Taylor approximation; must be a non-negative
%         integer
%
%   Outputs:
%      Et: a scalar representing the true absolute error of the Taylor
%          series approximation
%
%   Assumptions: 
%         - none.

%%
%% ME 2016 Section
%% Fall 2007
%%
%% AUTHOR:  Chris Paredis
%%

% start with some error checking
% the input n must be a non-negative integer scalar
% NOTE: it is important that the function stop executing when an error is
% detected.  The best way to accomplish this is by using the error()
% function
if length(n)>1
    error('the input must be a scalar');
end
if n<0
    error('the input must be a non-negative integer');
end
% make sure the input is an integer
% NOTE: it is equally acceptable to generate an error if n is not an
% integer
n = floor(n);

% this implementation of the function avoids most of the problems
% associated with floating point arithmetic.
% Rather than computing the difference between exp(x) and the Taylor
% series, it computes the remainder as the series from n+1 to infinity.
% When the error is small compared to exp(x), we can still compute this
% series accurately because it will consist of a sum a small number (rather
% than a difference between two large numbers).  There are a few additional
% obstacles to overcome:
% 1) Although each term in the remainder series is small, the computation
% of the term involves a division of two very large numbers: x^n and n!.
% Both of these numbers would cause and overflow in floating point
% arithmetic (Inf).  We could therefore apply the following trick:
%    x^n/n! = exp(log(x^n) - log(n!)).
% The logarithms are small numbers, and their difference is really small,
% so that floating point overflow never occurs.
% A different way to avoid the same problem is to recognize that the term
% can be rewritten as: (x/1)*(x/2)*(x/3)*...*(x/n).  Each of the numbers in
% this product is relatively small so that the computation does not result
% in an overflow.
% 2) Rather than evaluating the series from n+1 to infinity, we stop as
% soon as the additional terms are so small that they get chopped off in
% the floating point addition.


%% Approach 1: using term_n = exp(log(x^n)-log(n!))
% initialize the sum variable and the iterator
Et = 0;
i = n+1;
% we use a while loop because we don't know in advance how many terms in
% the sum need to be included.  Rather than testing for the exit criterion
% in the while loop, we test inside the loop.
while 1
    % the log of x^i is i*log(x)
    log_power_term = i*log(x);
    % using the colon operator, one can compute the log of the factorial
    % as follows
    log_fact = sum(log(1:i));

    % add the current term to the remainder sum
    Et_new = Et + exp(log_power_term-log_fact);

    % check whether the current term gets chopped off; if so, then we have
    % reached the best possible accuracy and can exit the function
    if Et_new == Et
        break;
    end

    % update the iteration variable and go to the next iteration
    Et = Et_new;
    i = i+1;
end

%% Approach 2: using term_n = (x/1)*(x/2)*...*(x/n)
Et2=0;  % use Et2 to distinguish from Et
k=n+1;
% start by calculating the (n+1)th term of the Taylor series
term = prod(x./(1:(n+1)));
% exit the while loop when the addition of term is completely chopped off
while Et2+term > Et2
    Et2 = Et2 + term;
    k=k+1;
    term = term * x/k;
end

% just out of curiosity you can take the difference between the two
% approaches and see whether you get the same result.  You will actually
% notice that there is a small difference again due to chopping.
% uncomment the next line only for debugging purposes
%diff = Et-Et2