function [ y ] = etetr (x)
%y=etetr(x) 
% e <tetrate> x, which equals to exp(0), iterated x times
% x may be real, in this case smooth approximation will be used

%polynomial for smooth approximation
P=[-0.032486264307116,...
   0.161609516997898,...
  -0.272374056966456,...
   0.365751133821221,...
  -0.310443586372316,...
   1.087943256826769,...
   0.000000000000000];


xi=floor(x);
xf=x-xi;

y=polyval(P,xf);

%now iterating exp and log
x_neg=xi<0;
if any(x_neg)
  y(x_neg)=log(y(x_neg));
end

%iterating exp
do
  x_plus=xi>0.5; %those x that are 1 and more
  if ~any(x_plus)
    return
  end
  y(x_plus)=exp(y(x_plus));
  xi=xi-1;
until false

endfunction
