graphics_toolkit("gnuplot");

%eps 
eps
1+eps>1
1+eps*0.5>1
eps(single(1))

%policz 2^(-k) takei ze 1+x>0 a fl(1+x/2)=0
x=1;k=0;
while(x+1>1)
 x=x*0.5;
 k++;
end
--k
x*2
%to samo w pojedynczej precyzji single
 x=single(1);k=0;
while(x+1>1)
 x=x*0.5;
 k++;
end
--k
x*2

%wykresy (x+a)-a dla a=3^k k=3,...,50

a=27;k=3;
M=50;
x=linspace(0,1,1000);
for k=3:38,
  com=[";k=",num2str(k),";"];
  if(k>26) 
   plot(x,(x+a)-a,com);pause(1);
  end 
  
  a*=3;
end  

%f(x)=x-sqrt(x.*x+1)=-1/(x+sqrt(x.*x+1)) na 2 sposoby
%dla x=10^k
N=40;
x=logspace(0,10,N)';
f1=x-sqrt(x.*x+1)
f2=-1./(x+sqrt(x.*x+1))
ind=(1:N)';
ae=[f1-f2];
re=[ind,ae./f2]


%wykres (x-2)^4 na 2 sposoby na [2-a,2+1] a=1,1e-1..,1e-4
a=1;
for k=0:4,
 x=linspace(2-a,2+a,1000);
 w1=(x-2).*(x-2); w1=w1.*w1; %(x-2)^4
 w2=x.*x.*x.*x-8*x.*x.*x+24*x.*x-32*x+16;
 plot(x,w1,x,w2,";x^4..16;"); pause(1);
 a*=0.1;
end 
%to samo w pojedynczej prezycji single



%wykres (x-2)^4 na 2 sposoby na [2-a,2+1] a=1,1e-1..,1e-4
a=1;
for k=0:4,
 x=linspace(2-a,2+a,1000);x=single(x);
 w1=(x-2).*(x-2); w1=w1.*w1; %(x-2)^4
 w2=x.*x.*x.*x-8*x.*x.*x+24*x.*x-32*x+16;
 plot(x,w1,";single;",x,w2,";x^4..16;"); pause(1);
 a*=0.1;
end 

%nasz exponent
function y=nexp(x,N=1e4)
 y=zeros(size(x));
 xx=ones(size(x)); %xx=x^k/k! k=0,...
 for k=1:N,
   y+=xx;
   xx=xx.*x/k;
 end
end  


kk=-50:2:50;

ekk=exp(kk);
ae=nexp(kk)-ekk;
re=ae./ekk;
[kk',re']

%poprawiamy nasz exponent
function y=nexp1(x,N=1e4)
 xp=(x>=0);
 xm=(x<0);
 x=xp.*x-xm.*x; % x(k)=x(k) gdy x(k)>=0 x(k)-x(k) wpp 
 y=zeros(size(x));
 xx=ones(size(x)); %xx=x^k/k! k=0,...
 for k=1:N,
   y+=xx;
   xx=xx.*x/k;
 end
 y=xp.*y+xm./y;
end  

kk=-50:2:50;

ekk=exp(kk);
ae=nexp1(kk)-ekk;
re=ae./ekk;
[kk',re']
