

#zadania z fl mo 15/16

#zad 1 eps maszynowy

a=1;n=0;
while(a+1>1)
 n++;
 a*=0.5;
endwhile
disp("najwieksze n : fl(2^(-n)+1)>1");
n--
a=a*2

a-eps #porownujemy z epsilonem maszynowym

#w arytmetyce pojedynczej precyzji

a=single(1);n=0;
while(a+1>1)
 n++;
 a*=0.5;
endwhile
disp("single: najwieksze n : fl(2^(-n)+1)>1");
n--
a=a*2


#wykres (a+1)-a na [0,1] a=3^k
a=3;x=linspace(0,1,500);
M=34;
e=zeros(M,1);
for k=1:M,
 a*=3; y=(a+x)-a;
 if(k>23) plot(x,y,[";y=(x+a)-a a=",num2str(a)," n=",num2str(k), ";"]);
 pause(1);
 e(k)=norm(x-y,'inf');
endif
endfor
plot(e,";||(x-3^n)+3^n-x||;");
pause(2);
#wykres (a+1)-a na [0,1] a=3^k ale w arytmetyce pojedynczej precyzji
a=3;x=linspace(0,1,500);
M=20;
e=zeros(M,1);
x=single(x);a=single(a);
for k=1:M,
 a*=3; y=(a+x)-a;
 if(k>8) 
 plot(x,y,[";y=(x+a)-a a=",num2str(a)," n=",num2str(k), ";"]);pause(1);
endif
 e(k)=norm(x-y,'inf');
endfor
plot(e,"r*;||(x-3^n)+3^n-x|| single;");
hold on;
plot(e);pause(2);
hold off




#x-sqrt(1+x.*x) na 2 sposoby
x=10.^(1:10);x=x';
f1=x-sqrt(1+x.*x);
f2=-1.0./(x+sqrt(1+x.*x));
[f1 f2]
e=f1-f2;
ew=e./f2;
[e ew]
#to samo w pojedynczej precyzji 
x=single(x);
fs1=x-sqrt(1+x.*x);
fs2=-1.0./(x+sqrt(1+x.*x));
[fs1 fs2]
es=fs1-fs2;
ews=e./fs2;
[es ews]




#wykres (x-2)^4 na 2 sposoby na [2-a,2+a] dla a=1e-2,1e-3,1e-4
f1=@(x) (x-2).*(x-2).*(x-2).*(x-2);
f2=@(x) x.*x.*x.*x - 4*2.*x.*x.*x + 6*2.*2.*x.*x -4*2.*2.*2.*x + 2.*2.*2.*2;

a=1e-1;
for k=1:4,
a*=0.2; 
x=linspace(2-a,2+a,200);
plot(x,f1(x),";(x-2)^4;",x,f2(x),[";x^4+..16 na [2-a,2+a] z a=",num2str(a),";"]);
pause(2);
endfor

#to samow pojed. precyzji
a=4;
for k=1:5,
a*=0.4; 
x=single(linspace(2-a,2+a,200));
plot(x,f1(x),";single (x-2)^4;",x,f2(x),[";x^4+..16 na [2-a,2+a] z a=",num2str(a),";"]);
pause(1.5);
endfor


#exponent z rowiniecia w szereg exp(x)=sum_k x^k/k!
function [y]=nexp(x,N=50)
 xx=y=1;
 x=x; 
 for k=1:N,
  xx=(xx.*x)/k;
  y=y+xx;
 endfor
endfunction
x=-15:15;
ex=exp(x');
ae=ex-nexp(x');
re= ae./ex;
disp("bledy bezwzg i wzgledny obl. exp(x) jako sum_{k=0}^N x^k/k!");
[ae,re]

#poprawiamy tzn dla x>=0 uzywamy sum_{k=1}^N x^k/k! a dla 
#x<0 1.0/( sum_{k=1}^N (-x)^k/k!)
function [y]=nexp1(x,N=50)
 xu=(x<0); 
 xd=(x>=0);
 x=abs(x);
 xx=1;y=1; 
 for k=1:N,
  xx=(xx.*x)/k;
  y=y+xx;
 endfor
 y=xd.*y + xu.*(1.0./y);
endfunction
x=-15:15;
ex=exp(x');
ae1=ex-nexp1(x');
re1= ae1./ex;
disp("bledy ale dla alg poprawionego tzn dla x<0 liczymy 1/(sum_k x^k/k!)");
[ae1,re1]

#to samo w pojedynczej precyzji
function [y]=nexps(x,N=50)
 xx=y=single(1);
 x=single(x); 
 for k=1:N,
  xx=(xx.*x)/k;
  y=y+xx;
 endfor
endfunction
x=-15:15;
ex=exp(x');
aes=ex-nexps(x');
res= aes./ex;
disp("single: bledy obl. exp(x) jako sum_{k=0}^N x^k/k!");
[aes,res]

#poprawiamy tzn dla x>=0 uzywamy sum_{k=1}^N x^k/k! a dla 
#x<0 1.0/( sum_{k=1}^N (-x)^k/k!)
function [y]=nexps1(x,N=50)
 xu=(x<0); 
 xd=(x>=0);
 x=single(abs(x));
 xx=y=single(1); 
 for k=1:N,
  xx=(xx.*x)/k;
  y=y+xx;
 endfor
 y=xd.*y + xu.*(1.0./y);
endfunction
x=-15:15;
ex=exp(x');
aes1=ex-nexps1(x');
res1= aes1./ex;
disp("single: bledy ale  dla x<0 liczymy 1/(sum_k x^k/k!)");
[aes1,res1]


#liczymy calki I(n)=\int_0^1 x^{n-1}/(x+5) n=1,2,.. ze wzoru
# I(n+1)+5*I(n)=1/n i I(1)=log(6/5) dla n=1,...,21
#mamy 1/6 =n * \int_0^1 x^(n-1)/6 < n*I(n) < n * \int_0^1 x^(n-1)/5=1/5  
#czy te nierownosi zajda?
N=21;
I=zeros(N,1);
I(1)=log(6/5);
for k=2:N,
 I(k)=1/(k-1)-5*I(k-1);
endfor
ind=(1:21)';
disp("I(n)=int_0^1 x^n(x+5) z wzoru:I(n+1)=1/n-5*I(n)=1/n;I(1)=log(6/5)");
nI=ind.*I;

ok=[nI>1/6, nI<1/5];

if(sum(sum(1-ok))) 
 printf("Blad: jakas calka n*I(n) poza [1/6,1/5]\n"); 
 [i,j]=find(1-ok);
 for k=1:length(i),
   if(j(k)==1)
    printf("%d*I(%d)=%g <= 1/6=%g\n",i(k),i(k),nI(i(k)),1./6);
   else
    printf("%d*I(%d)=%g >= 1/5=%g\n",i(k),i(k),nI(i(k)),0.2);
   endif  
 end
endif



#liczymy do tylu biorac I(30)=0.5*(1/6 +0.2)
It=zeros(30,1);
# It(n)=(1/n-It(n+1))/5 czy  It(1)=log(6/5) dla n=31,..,1?
It(31)=0.2*(1/6 + 1/5);
for k=30:-1:1,
 It(k)=0.2*(1/k - It(k+1));
endfor
It=It(1:21);
disp("Liczy calki w tyl startujac z It(31)=0.2*(1/6 + 1/5); ");
disp("z tego samego wzoru do tylu: It(n)=(1/n-It(n+1))/5\n"); 
printf("Blad wz It(1)-log(6/5)=%g na poziomie bl.fl..",(It(1)-log(6/5))/log(6/5));
disp("czy pomiedzy 1/6< n*It(n) < 1/5 dla n=1,..,21?");
ind=(1:21)';
nI=ind.*It;

ok=[nI>1/6, nI<1/5];

if(sum(sum(1-ok))) 
   printf("Blad: jakas calka n*It(n) poza [1/6,1/5]\n"); 
 else 
   printf("Wszystkie calki   n*It(n) w [1/6,1/5]\n");
endif






