

#zad 1 eps czy 1+eps > 1 w fl?

1+eps>1
1+eps/2>1

single(1+eps)>1

#zad 2 szukamu najwiekszego k takiego ze
#1+2^{-k}>1 ale 1+2^{-k-1}=1 w fl

a=1;k=0;
while((a+1)>1)
a=a/2;
k++;
endwhile
a=a*2
k--

#to samo single

a=single(1);k=0;
while((a+1)>1)
a=a/2;
k++;
endwhile
a=a*2
k--

#zad 3 wykres y=x=(x+a)-a na [0,1] dla duzych a=5^k k=1,2..,50 itd

x=linspace(0,1);
a=mn=5;
for k=2:35,
a=a*mn;
plot(x,(x+a)-a,[";",num2str(k),";"]);
pause(0.3);
endfor

#zad 4 f(x)=x-sqrt(1+x*x)=-1/(sqrt(1+x*x)
#liczony na 2 sposoby dla x=10^k k duze

x(1)=10;
for k=2:10,
x(k)=x(k-1)*10;
endfor
a=sqrt(1+x.*x);
y1=x-a;
y2=-1./(x+a);
ba=abs(y1-y2);
bw=ba/abs(y2);

disp("blad wzgl obliczania x-sqrt(1+x*x)=-1/(sqrt(1+x*x)) na 2 sposoby");
bw


#zad 5 wykres (x-2)^4 na 2 sposoby
aa=[1e-1 1e-2 1e-3 1e-4];
for k=1:4,
 a=aa(k);
 x=linspace(2-a,2+a);
 b=(x-2);
 y1=b.*b.*b.*b;
 y2=x.*x.*x.*x - 8*x.*x.*x + 24*x.*x-32*x+16;
 op1=[";double  a=",num2str(a),"   os OX;"];
 plot(x,x-x,op1,x,y1,";(x-2)^4;",x,y2,";x^4-8x^3+24x^2-32x+16;")
 pause(3);
endfor


#to samo w pojedynczej precyzji
aa=aa*10;
for k=1:4,
 a=aa(k);
 x=linspace(2-a,2+a);
 x=single(x);
 b=(x-2);
 y1=b.*b.*b.*b;
 y2=x.*x.*x.*x - 8*x.*x.*x + 24*x.*x-32*x+16;
 op1=[";single  a=",num2str(a),"   os OX;"];
 plot(x,x-x,op1,x,y1,";(x-2)^4;",x,y2,";x^4-8x^3+24x^2-32x+16;")
 pause(3);
endfor


#zad 6 exp(x) z rozwiniecia w szereg
#przyblizamy nexp(x,N)=sum_{k=0}^N x^k/k!

function y=nexp1(x,N=1000)
 y=ones(size(x));
 a=y;
 for k=1:N,
   a=a.*x/k;
   y=y+a;
 endfor
endfunction

#testy dla x=-100,-90,..,0,10,..,100
disp("testy dla naszej f liczacej przyblizenie exp(x)z rozwiniecia")
x=-100:10:100;
n=length(x);
y1=exp(x);y2=nexp1(x);
bb=y1-y2;
bw=bb./y1
disp("dla x ujemnych duzy blad wzgledny");

#dla x ujemnych dziala zle ale mozna to poprawic bo exp(x)1/exp(-x) 

function y=nexp2(x,N=1000)
#tylko wersja skalarna
 yp=yu=zeros(size(x));
 xp=x.*(x>=0); #xp ma tylko nieujemne skaldowe x a  zera na pozostalych pozycjach 
 xu=x.*(x<0);#xu ma tylko ujemne skaldowe x a  zera na pozostalych pozycjach 
 yp=nexp1(xp,N).*(x>=0); #liczyme exp(x) ale zerujac te pozycje gdzie x byl ujemny
 yu=(1./nexp1(-xu,N)).*(x<0);#to samo zerujac pozycje gdzie x nieujemny
 y=yp+yu;
endfunction
#powtarzamy testy

#testy dla x=-100,-90,..,0,10,..,100
disp("testy dla naszej f liczacej przyblizenie exp(x) ulepszonej")
x=-100:10:100;
n=length(x);
y1=exp(x);y2=nexp2(x);
bb=y1-y2;
bw=bb./y1
disp("teraz blad wzgledny maly tez dla ujemnych x ")

