%% Simple f4
clear all;
clc;

function [f] = fun(x)
  f = exp(x)./( (sin(x)).^3  + (cos(x)).^3 );
end

function [f] = d1(h,x0)
  x = [x0-h-h,x0-h,x0,x0+h,x0+h+h];
  y = [ fun(x(1)),fun(x(2)),fun(x(3)),fun(x(4)),fun(x(5))];
  
  f = (y(1)./12-2.*y(2)./3+2.*y(4)./3-y(5)./12)./h;
end

function [f] = d2(h,x0)
  x = [x0-h-h,x0-h,x0,x0+h,x0+h+h];
  y = [ fun(x(1)),fun(x(2)),fun(x(3)),fun(x(4)),fun(x(5))];
  
  f = (-y(1)./12+4.*y(2)./3-5.*y(3)./2+4.*y(4)./3-y(5)./12.0) ./h.^2;
end

function [f] = d3(h,x0)
  x = [x0-h-h-h,x0-h-h,x0-h,x0,x0+h,x0+h+h,x0+h+h+h];
  y = [ fun(x(1)),fun(x(2)),fun(x(3)),fun(x(4)),fun(x(5)),fun(x(6)),fun(x(7))];
  
  f = (y(1)./8-y(2)+13.*y(3)./8-13.*y(5)./8+y(6)-y(7)./8) ./h.^3;
end


function [f] = d4(h,x0)
  x = [x0-h-h-h,x0-h-h,x0-h,x0,x0+h,x0+h+h,x0+h+h+h];
  y = [ fun(x(1)),fun(x(2)),fun(x(3)),fun(x(4)),fun(x(5)),fun(x(6)),fun(x(7))];
  
  f = (-y(1)./6+2.*y(2)-13.*y(3)./2+28.*y(4)./3-13.*y(5)./2+2.*y(6)-y(7)./6) ./h.^4;
end

% Program flow
h_isDisturbed = true;
%h_isDisturbed = false;

stat = 1E5

dx = 1E-6; %plati co inde

h0_d1 = 1E-4;
h0_d2 = 1E-3;
h0_d3 = 1E-3;
h0_d4 = 1E-2;
%h0_d4 = 1E-2/2;
  
x0 = 0;

truth_d1 = 1;
truth_d2 = 4;
truth_d3 = 4;
truth_d4 = 28;

nHistoBins = 50;

kahanSum = [0,0,0,0];
kahanC = [0,0,0,0];

for i=0:stat-1
  
    x = x0;
    di = [0,0,0,0];
  
    if (h_isDisturbed)
        h = h0_d1 .*( unifrnd(0,1)+0.5 );
        di(1) = d1(h,x);

        h = h0_d2 .*( unifrnd(0,1)+0.5 );
        di(2) = d2(h,x);

        h = h0_d3 .*( unifrnd(0,1)+0.5 );
        di(3) = d3(h,x);

        h = h0_d4 .*( unifrnd(0,1)+0.5 );
        %h = h0_d4 + 1E-10.*( 2.*unifrnd(0,1)-1.0 );
        di(4) = d4(h,x);

        for j=1:4
            kahanY(j) = di(j) - kahanC(j);
            kahanT(j) = kahanSum(j) + kahanY(j);
            kahanC(j) = (kahanT(j)-kahanSum(j)) - kahanY(j);
            kahanSum(j) = kahanT(j);
        end
  
    else
        x = x0 + dx.*( 2.*unifrnd(0,1) -1);
    
        di(1) = d1(1E-4,x);
        di(2) = d2(1E-3,x);
        di(3) = d3(1E-3,x);
        di(4) = d4(1E-2,x);

    for j=1:4
        kahanY(j) = di(j) - kahanC(j);
        kahanT(j) = kahanSum(j) + kahanY(j);
        kahanC(j) = (kahanT(j)-kahanSum(j)) - kahanY(j);
        kahanSum(j) = kahanT(j);
    end
  
    end
  
end

df1_si = d1(1E-4,x0);
d1_si = abs(truth_d1 -df1_si);
df1_av = kahanSum(1)/stat;
d1_av = abs(truth_d1 -df1_av);
		
df2_si = d2(1E-3,x0);
d2_si = abs(truth_d2 -df2_si);
df2_av = kahanSum(2)/stat;
d2_av = abs(truth_d2 -df2_av);

df3_si = d3(1E-3,x0);
d3_si = abs(truth_d3 -df3_si);
df3_av = kahanSum(3)/stat;
d3_av = abs(truth_d3 -df3_av);

df4_si = d4(1E-2,x0);
d4_si = abs(truth_d4 -df4_si);
df4_av = kahanSum(4)/stat;
d4_av = abs(truth_d4 -df4_av);

%

fprintf("\n");
fprintf("First Derivative\n");
fprintf("Truth = %d\n",truth_d1);
fprintf("AV = %d, delta = %d\n",df1_av,d1_av);
fprintf("Single = %d, delta = %d\n",df1_si,d1_si);
fprintf("Sing/AV delta = %d\n",d1_si./d1_av);

fprintf("\n");
fprintf("Second Derivative\n");
fprintf("Truth = %d\n",truth_d2);
fprintf("AV = %d, delta = %d\n",df2_av,d2_av);
fprintf("Single = %d, delta = %d\n",df2_si,d2_si);
fprintf("Sing/AV delta = %d\n",d2_si./d2_av);

fprintf("\n");
fprintf("Third Derivative\n");
fprintf("Truth = %d\n",truth_d3);
fprintf("AV = %d, delta = %d\n",df3_av,d3_av);
fprintf("Single = %d, delta = %d\n",df3_si,d3_si);
fprintf("Sing/AV delta = %d\n",d3_si./d3_av);

fprintf("\n");
fprintf("Fourth Derivative\n");
fprintf("Truth = %d\n",truth_d4);
fprintf("AV = %d, delta = %d\n",df4_av,d4_av);
fprintf("Single = %d, delta = %d\n",df4_si,d4_si);
fprintf("Sing/AV delta = %d\n",d4_si./d4_av);
