

#eps maszynowy
a=1;
n=0;
while(a+1>1)
 a*=0.5;
 n++;
end
n
1+a>1
a=2*a;1+a>1
--n



#eps maszynowy dla pojedynczej precyzji
a=single(1);
n=0;
while(a+1>1)
 a*=0.5;
 n++;
end
n
1+a>1
a=2*a;1+a>1
--n


%wykres (1-a)+a dla a=5^n
graphics_toolkit("gnuplot"); 
a=125; 
x=linspace(0,1,300);
%x=single(x);
for k=5:25,
  com=[";k=",num2str(k),";"];
  if(k>19)%(k>8) dla pojed precyzji
   plot(x,(x+a)-a,com);pause(3);
  end 
  a*=5;
end  

%wykres (x-2)^4 na (2-a),2+a) a=10^(-k_ k=1,2,3,4

a=0.1;
for k=1:4,
 x=linspace(2-a,2+a,300);
 y1=(x-2).*(x-2).*(x-2).*(x-2);
 y2=x.*x.*x.*x-8.*x.*x.*x+24*x.*x-32*x+16;
 plot(x,y1,[";(x-2)^4a=1e-",num2str(k),");"],x,y2,";x^4-..+16;");
 pause(3);
 a*=0.1; 
end  


a=1;
for k=1:4,
 x=linspace(2-a,2+a,300);
 x=single(x);
 y1=(x-2).*(x-2).*(x-2).*(x-2);
 y2=x.*x.*x.*x-8.*x.*x.*x+24*x.*x-32*x+16;
 plot(x,y1,[";(x-2)^4a=1e-",num2str(k-1),");"],x,y2,";x^4-..+16;");
 pause(3);
 a*=0.1; 
end  


%nasz exponent
function y=nexp(x,N=1000)
%liczymy exp z rozwiniecia tzn nexp(x)=\sum_{k=0}^N x^k/k!
y=1;
xx=1;
 for k=1:N,
   xx=xx.*x/k;
   y+=xx;
 end
end    
nexp(1)-e
nexp(0)-1
(nexp(100)-exp(100))/exp(100)
(nexp(-1)-exp(-1))/exp(-1)
x=-100:10:100;x=x';
%x=single(x);
y1=nexp(x);
y=exp(x);
ae=y1-y;
re=ae./y;
[x,re]

%poprawiamy nasza funkcje 

%nasz exponent
function y=nexp1(x,N=1000)
%liczymy exp z rozwiniecia tzn nexp(x)=\sum_{k=0}^N x^k/k!
 ujem=0;
 if(x<0)
  x=-x;ujem=1;
 end 
 y=1;
 xx=1;
 for k=1:N,
   xx=xx.*x/k;
   y+=xx;
 end
 if(ujem) 
  y=1/y; 
 end 
end    
nexp1(1)-e
nexp1(0)-1
(nexp1(100)-exp(100))/exp(100)
(nexp1(-100)-exp(-100))/exp(-100)

x=-100:10:100;x=x';
x=single(x);
for k=1:length(x),  
 y1(k)=nexp1(x(k));
end 
y=exp(x);
ae=y1-y;
re=ae./y;
format long e
[x,re]


%liczymy calke I(n)=\omt_0^1 x^n/(x+5) rekurencyjnie 
%ze wzoru I(0)=log(6/5) I(n)=1/n-I(n-1) 
%oblicz I(3) sprawdz czy 1/6<=(n+1)I(n)<=1/5
%przenumerujemy J(k)=I(k-1) k=1,..,31
M=30;
J=Jr=zeros(M+1,1);
Testr=Test=-ones(length(J),2);
J(1)=log(6/5);
for k=1:M,
  J(k+1)=1/k-5*J(k);  %I(k)
end
Test(:,1)=(np1J=(1:M+1)'.*J)>=1/6; 
Test(:,2)=(1:M+1)'.*J<=1/5;
format short
disp("[k, czy k*I(k-1)>=1/6, czy k*I(k-1)<=1/5, czy I(k) w przediale]");
[(1:M+1)',Test,sum(Test')'==2]

%liczymy w tył 

Jr(M+1)=1/(6*(M+1));
for k=M:-1:1,
%  J(k+1)=1/k-5*J(k);  %I(k)
   Jr(k)=(1/k-Jr(k+1))/5;
end
disp("[k, czy k*I(k-1)>=1/6, czy k*I(k-1)<=1/5, czy I(k) w przediale]");
Testr(:,1)=(np1Jr=(1:M+1)'.*Jr)>=1/6; 
Testr(:,2)=(1:M+1)'.*Jr<=1/5;
[(1:M+1)',Testr,sum(Testr')'==2]
J(1)-Jr(1) %na poziome dokladnosci arytmetyki
(J-Jr)(1:10)


%wykresy (x-2)^2 liczonego na 2 sposoby
a=0.1;
f1=@(x) (x-2).*(x-2);
f2=@(x) x.*x-4*x+4;

for k=1:4,
 x=single(linspace(2-a,2+a,300));
 plot(x,f1(x),[";(x-2)^2 a=1e-",num2str(k),");"],x,f2(x),";x^2-4x+4;");
 pause(3);
 a*=0.1; 
end  




a=1e-4;
f1=@(x) (x-2).*(x-2);
f2=@(x) x.*x-4*x+4;

for k=1:7,
 x=linspace(2-a,2+a,300);
 plot(x,f1(x),[";(x-2)^2 a=",num2str(a),");"],x,f2(x),";x^2-4x+4;");
 pause(2);
 a*=0.2; 
end  
