%% Simple fourier
clear all;
clc;

function [f] = fn(x)
  f = 1 - cos(x)./3 + cos(2.*x)./7 - cos(3.*x)/11;
end

%% PROGRAM FLOW

stat = 1E5
eps = 1E-15;
nHistoBins = 50;

%
csf = [1,-1/3,1/7,-1/11];
d = 1E-4;
x = [1,1+d,1+d+d,1+d+d+d];
for i=1:4
  y(i) = fn(x(i));
end
y = transpose(y);

Mat = zeros(4,4);
for i=1:4
  Mat(i,1) = 1;
  Mat(i,2) = cos( x(i) );
  Mat(i,3) = cos( 2.*x(i) );
  Mat(i,4) = cos( 3.*x(i) );
end

iMat = inv(Mat);
cfs_single = iMat*y;

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

for i=0:stat-1
  pertM = zeros(4,4);
  
  for k=1:4
    for l=1:4
      pertM(k,l) = Mat(k,l) + eps.*( 2.*unifrnd(0,1) -1 );
    end
  end
  
  iPertM = inv(pertM);
  cPert = iPertM*y;
  
  for j=1:4
    kahanY(j) = cPert(j) - kahanC(j);
    kahanT(j) = kahanSum(j) + kahanY(j);
    kahanC(j) = (kahanT(j)-kahanSum(j)) - kahanY(j);
    kahanSum(j) = kahanT(j);
  end
  
end


disp("Truth")
for i=1:4
  fprintf("%d\t",csf(i))
end

fprintf("\n")
fprintf("\n")
disp("AVERAGED")
for i=1:4
  av(i) = kahanSum(i)/stat;
  fprintf("%d\t",av(i));
end
fprintf("\n");
for i=1:4
  fprintf("%d\t", abs(av(i)-csf(i)) );
end

fprintf("\n")
fprintf("\n")
disp("SINGLE")
for i=1:4
  fprintf("%d\t",cfs_single(i));
end
fprintf("\n");
for i=1:4
  fprintf("%d\t",abs(cfs_single(i)-csf(i)));
end

fprintf("\n")
fprintf("\n")
disp("RATIOS")
for i=1:4
  fprintf("%d\t",  abs(cfs_single(i)-csf(i))./abs(av(i)-csf(i))  );
end
fprintf("\n");
