% Simple power
clear all;
clc;
% output_precision(16);
% pkg load statistics;

function [f] = a(n,x)
  if (n==1)
    f = x;
  elseif (n==2)
    f = x.*x;
  else
    f = (2.*x).*a(n-2,x) + (x - 2).*a(n-1,x);
  end
end

function [f] = binNb(xMin,xMax,nBins,x)
    dx = (xMax-xMin)./nBins;
    delta = x-xMin;
    f = floor(delta./dx);
end
    
function [f]	= createHisto(xMin, xMax, nBins)		
    histo = zeros(2,nBins);
    dx = (xMax-xMin)./nBins;
    histo(0+1,0+1) = xMin + dx./2.0;

    for i=0:nBins-1
        histo(0+1,i+1) = histo(0+1,0+1) + i.*dx;
        histo(1+1,i+1) = 0.0;
    end
    
    f = histo;		
end
    
function [f] = fillHisto(histo, x)		
    L = length(histo);
    dx = (histo(0+1,L-1+1)-histo(0+1,0+1))./(L-1.0);
    xMin = histo(0+1,0+1)-dx./2.0;
    xMax = histo(0+1,L-1+1)+dx./2.0;
    bin = binNb(xMin,xMax,L,x);
    if (0<=bin && bin<L) 
        histo(1+1,bin+1) =  histo(1+1,bin+1) + 1;
    end
    
    f = histo;
end
    
function [f] = normHisto(histo)		
  L = length(histo);
  sum = 0.0;
  for i=0:L-1
    sum = sum + histo(1+1,i+1);
  end

  if (sum==0)
    disp("Histogram cannot be normalized, is empty!");
    exit;
  end

  for i=0:L-1
    histo(1+1,i+1) = histo(1+1,i+1)./sum;
  end
	
  f = histo;	
end

function [f] = printHisto(histo)
  disp("");
	L = length(histo);
  for i=0:L-1	
    xVal = histo(0+1,i+1);
    yVal = histo(1+1,i+1);
    fprintf(" ( %.16e , %.16e ) \n", xVal,yVal );
  end
  
  disp("");

	disp("");
	sum = 0.0;
  for i=0:L-1
    sum = sum + histo(0+1,i+1)*histo(1+1,i+1);
  end
  fprintf("Histo average is = %.16e\n", sum);
  
end


%% PROGRAM FLOW

stat = 1E6
eps = 1E-7;
x0 = 1./3;
n = 15;
truth = 0.0000000696917193762563239137308507191523368295578192819843351;
width = 5E-13;
nHistoBins = 50;

kahanSum = 0.0;
kahanC = 0.0;

xMin = truth-width;
xMax = truth+width;
histo = createHisto(xMin,xMax,nHistoBins);

for i=0:stat-1
  x = x0 + eps.*( 2.*unifrnd(0,1) - 1);
  y = a(n,x);
  
  kahanY = y - kahanC;
  kahanT = kahanSum + kahanY;
  kahanC = (kahanT-kahanSum) - kahanY;
  kahanSum = kahanT;
  
  histo = fillHisto(histo,y);
end

single = a(n,x0);
d_st = abs(single - truth);

avg = kahanSum./stat;
d_at = abs(avg-truth);

d_rat = d_st./d_at;
fprintf("Truth = %.16e\n",truth);
fprintf("AV = %.16e, delta = %.16e\n",avg,d_at);
fprintf("Single = %.16e, delta = %.16e\n",single,d_st);
fprintf("Sing/AV delta = %.16e\n",d_rat);
fprintf("\n");

histo = normHisto(histo);
printHisto(histo);
