function p1_and_2
% Plug in t^2 into the formula
% sin(x) = x/1! - x^3/3! + x^5/5!
% to prove that the series
% is the one we used in calculating
% ynew in the function IntSin2.

% Problem 1

% 1a)
% See the function antideriv below.

% 1b)
f = @(x) sin(x.^2);

antideriv_integral_test = @(x) ...
    antideriv(f, 0, x, 1e-9);
antideriv2_integral_test = @(x) ...
    antideriv2(f, 0, x, 1e-9);

% Note that it is only fair if we use the
% same tolerance 1e-9 for all tests.
% It would be like comparing apples to
% oranges if we used 1e-8 for quad
% and 1e-9 for IntSin2 as the problem suggests.
x = .02:.02:5;
tic
y_antideriv = antideriv_integral_test(x);
t_antideriv = toc
% about 3.8750s
tic
y_antideriv2 = antideriv2_integral_test(x);
t_antideriv2 = toc
% about 0.1410s

% 1c)
figure(1)
plot(x, y_antideriv2);

% Problem 2

% 2a)
% See the function IntSin2 below.

% 2b)
tic
[y_IntSin2, ntrm] = IntSin2(x);
t_IntSin2 = toc
% about 0.0150s

% 2c)
nterm_0pt4 = ntrm(20)
% 3
nterm_4pt0 = ntrm(200)
% 27

% 2d)
fprintf('x    antideriv2   IntSin2\n');
for i=10:10:length(x)
fprintf('%.02f %1.10f %1.10f\n', x(i), y_antideriv2(i), y_IntSin2(i));
end

% 2e)

F_7_antideriv2 = antideriv2_integral_test(7)
[F_7_IntSin2, ntrm, terms] = IntSin2(7, NaN, 1, 1);
F_7_IntSin2
figure(1)
semilogy(abs(terms));
title('terms in computing F(7) on a logarithmic scale');
firstbig = find(abs(terms)>1e+16, 1)
lastbig = find(abs(terms)>1e+16, 1, 'last')

% We notice that the terms being added together
% grow enormously before the Taylor series finally
% converges. The terms finally converge after n = 70.
% For n = 13 to 38, the terms overflow the
% floating point data type.
% Omitting these terms won't help either, since the
% correct value of the function relies on their presence
% when computed by a series.

% Obviously, we would need to
% take the argument x^2 modulo 2*pi in order to
% avoid this overflow.
% In fact, it is evident we may replace
% the value of the argument with a value between
% 0 and pi/4 by using trigonometry.

% This next figure illustrates what happens when
% the factorial value in the sine series
% overflows the arithmetic size.
% Print out figure 1 and figure 2.
figure(2)
x = .02:.05:6.3;
y = antideriv2_integral_test(x);
y2 = IntSin2(x, 1e-9);
plot(x, y, 'b', x, y2, 'g');
title('failure of series convergence in finite representation');

function y = antideriv(f, a, x, tol)
if nargin < 4
    tol = 1e-8;
end
x = sort(x);
y = zeros(size(x));
for i=1:length(x)
    y(i) = quad(f, a, x(i), tol);
end

function y = antideriv2(f, a, x, tol)
if nargin < 4
    tol = 1e-8;
end
x = sort(x);
y = zeros(size(x));
y(1) = quad(f, a, x(1), tol);
for i=2:length(x)
    y(i) = y(i-1) + quad(f, x(i-1), x(i), tol);
end


function [y, ntrm, extract_series] = IntSin2(x, tol, extract_v, ignore)
% does less computations when possible
% for example, x = 0 takes 0 terms to
% evaluate.
if nargin < 2 || isnan(tol)
    tol = 1e-8;
end
y = zeros(size(x));
ntrm = y;
laststep = 1;
ynew = x.^3/3;
max = 50; % maximum steps--otherwise our answer is
% inaccurate anyway, since (2*50-1)! is too big.
if nargin > 3
    max = 500;
end
for i=2:max
    if nargin > 2
        extract_series(i-1) = ynew(extract_v+laststep-1);
    end
    y(laststep:end) = y(laststep:end) + ynew;
    k = find(abs(ynew)>tol, 1);
    % we will only need to keep computing for the
    % part of the vector x beyond k.
    if k
        ntrm(laststep:laststep+k-2) = i-2;
        laststep = laststep + k - 1;
        xp = x(laststep:end);
        ynew = (-1)^(i+1)/factorial(2*i-1)/(4*i-1)*xp.^(4*i-1);
    else
        ntrm(laststep:end) = i-2;
        break;
    end
end