

%epsilon maszynowy

a=eps
a+1>1 
%ok a teraz dla polowy a?
0.5*a+1>1

%szukamy najmniejszego n i a takiego ze a=2^(-n) i a+1>1
disp("eps w podwojnej precyzji");
a=1;n=0;
while((a+1)>1)
 a=a*0.5;
 n=n+1;
end
n-1
a*2

%w pojedynczej precyzji
disp("eps w pojedynczej precyzji");
a=single(1);n=0;
while((a+1)>1)
 a=a*0.5;
 n=n+1;
end
n-1
a*2

disp("wykres f(x)(x+5^k)-5^k  w podwojnej precyzji");
M=26;
a=b=5;
x=linspace(0,1,1000);
for k=1:M,
  if(k>18)
   plot(x,(x+a)-a,[";double k=" num2str(k), ";"]);
   pause(1);
  end 
  a=b*a;
end  

disp("wykres f(x)(x+3^k)-3^k  w pojed precyzji");
M=20;
a=b=single(3);
x=linspace(0,1,1000);
for k=1:M,
  if(k>5)
   plot(x,(x+a)-a,[";single k=" num2str(k), ";"]);
   pause(1);
  end 
  a=b*a;
end  


%liczymy x=sqrt(1+x*x) na 2 sposoby dla x=5^k
disp("x-sqrt(1+x^2) na 2 sposoby w podwojnej precyzji");
M=20;
x=zeros(M,1);
x=5.^(0:M);x=x';
w=sqrt(1+x.*x);
w1=x-w
w2=-1.0./(x+w)
ea=w1-w2
er=ea./w2

%powtarzamy w pojedynczej precyzji
disp("x-sqrt(1+x^2) na 2 sposoby w podwojnej precyzji");
M=20;
x=zeros(M,1);
x=5.^(0:M);x=single(x)';
w=sqrt(1+x.*x);
w1=x-w
w2=-1.0./(x+w)
ea=w1-w2
er=ea./w2


%wykres (x-2)^4=x^4-...+64 na 2 sposoby
a=0.1;
for k=1:6,
  x=linspace(2-a,2+a,1000);
  f1=(x-2);
  f1=f1.*f1;f1=f1.*f1;
  f2=x.*x.*x.*x -8.*x.*x.*x+24*x.*x-32*x+16;
  plot(x,f1,[";double a=," num2str(a)," x-2)^4;"],x,f2,";x^4...;");
  pause(1);
  a*=0.2;
end  



%wykres (x-2)^4=x^4-...+64 na 2 sposoby
a=0.5;
for k=1:5,
  x=linspace(2-a,2+a,1000);
  x=single(x);
  f1=(x-2);
  f1=f1.*f1;f1=f1.*f1;
  f2=x.*x.*x.*x -8.*x.*x.*x+24*x.*x-32*x+16;
  plot(x,f1,[";single a=," num2str(a)," (x-2)^4;"],x,f2,";x^4...;");
  pause(1);
  a*=0.5;
end  


%exponent 
function y=nexp(x,N=100)
%nasz exponent z rzwiniecia exp(x) w szereg
%przyblizamy exp(x) przez \sum_{k=0}^N x^k/k!
disp("liczymy nasz exponent i porownujemy z exp(x) - f. octave'a");
xx=1.0;
kf=1;
y=0;
 for k=1:N,
  y=y+xx/kf;
  xx=xx.*x;%xx=x^k
  kf=kf*k; %kf=k!
 end
end
 
 x=2.^[-10:10];
 x=x';
 w=nexp(x);
 w1=exp(x);
 bb=w-w1
 bw=bb./w1
%powtorz z N=1000 i w arytmetyce single


%Liczymy I(n)=\int_0^1 x^(n-1)(x+5) n=1,2,..,20
%ta calka spelnia (1/6)(1/n)<=I(n)<=(1/5)(1/n) 
%mamy wzrory rekurencyjne: I(0)=log(6/5)
% I(n+1)=\int_0^1 (x^n+5*x^(n-1))/(x+5) - 5\int_0^1 x^(n-1)/(x+5)dx =
%       = 1/n -5*I(n) 
%albo I(n)=(1/n-I(n+1)/5 startując z I(n0)=(1/5)(1/n0) dla n0>>20 np n0=30  
M=25;
jpn=1.0./(1:M+10)';
I=It=zeros(M+10,1);
%pierwszy algorytm
I(1)=log(6/5);

for k=1:M-1,
  I(k+1)=jpn(k)-5*I(k);
end
%sprawdzamy nierownosci
ln=(I(1:M)>=(1/6)*jpn(1:M));
pn=(I(1:M)<=(1/5)*jpn(1:M));
wynikwprzod=[I(1:M),(1/6)*jpn(1:M),ln,(1/5)*jpn(1:M),pn]

%2 tyl iterujemy
n0=M+10;
It(M+10)=(1/5)/(M+10);
for k=M+9:-1:1,
  It(k)=(jpn(k)-It(k+1))/5;
end
ln=(It(1:M)>=(1/6)*jpn(1:M));
pn=(It(1:M)<=(1/5)*jpn(1:M));
wynikwtyl=[It(1:M),(1/6)*jpn(1:M),ln,(1/5)*jpn(1:M),pn]
%czy pierwszy
It(1)-log(6/5)



