%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% round_off.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcolo dell'eps di macchina

meps = 1.0;

while ( 1.0 + meps > 1.0 )
     meps = meps/2.0;
end

meps = 2.0*meps;

sprintf('eps calcolato:    %.15E\neps precalcolato: %.15E\n',meps,eps)
sprintf('meps-eps = %.15E\n',meps-eps)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Xmin.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcola il piu' piccolo numero floating point della forma xmin=2^-p

x=1;
p=0;
while x>0
  xmin = x;
  p=p+1;
  x=x/2;
end
p = p - 1;
xmin, p


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Xmax.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcola il piu' grande numero floating point della forma xmax=2^r

x=1;
r=0;
while x~=inf
  xmax = x;
  r=r+1;
  x=x*2;
end
r = r - 1;
xmax, r


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cancellazione.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Cancellazione numerica
%
%Calcola (1-x)^6 con le due formule:
%y1 = (1-x)^6
%y2 = 1-6x +15x^2 -20x^3 +15x^4 -6x^5 +x^6
%e confronta y1 e y2 in un intorno di zero
%
k = 0;
for delta = [0.1, 0.01, 0.005, 0.0025]
  h = delta/100;
  x = 1-delta:h:1+delta;
  y1 = (1-x).^6;
  y2 = 1 -6*x +15*x.^2 -20*x.^3 +15*x.^4 -6*x.^5 + x.^6;
  k = k+1;
  subplot(2,2,k)
     plot(x,y1)
     hold
     plot(x,y2,'g')
     axis([1-delta 1+delta -max(abs(y2)) max(abs(y2)) ])
end

disp('OSSERVARE che nell''intorno di 1')
disp(' (1-x)^6 si comporta meglio della forma espansa')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calc_norms.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% script per la visualizzazione delle sfere unitarie nelle norme
% 1,2 e infinito ed il calcolo celle norme di matrice indotte
echo on

R1 = 1.0;
R2 = 2.0;

N = 200;

%Dtheta = 2.0*pi/(N-1);
%
%% ciclo for
%for i=1:N
%    theta(i) = (i-1)*Dtheta;
%end
%
%% istruzione colon
%theta = 0:Dtheta:2.0*pi;

theta = linspace(0,2.0*pi,N);

%rho = R1 + (R2-R1)/(2.0*pi)*theta;
rho = linspace(R1,R2,N);

x = rho.*cos(theta);
y = rho.*sin(theta);

figure
plot(x,y),grid('on'),axis('square'),axis([-2 2 -2 2])
title('Spirale di partenza, coordinate cartesiane')

disp('Strike any key to continue')
pause
figure
polar(theta,rho),grid('on')
title('Spirale di partenza, coordinate polari')

disp('Ora calcolo le norme 1, 2 ed infinito dei vettori 0 -> punto spirale')
disp('Strike any key to continue')
pause

%v_norm_1(i) = norma_1 del vettore [x(i);y(i)]
v_norm_1   = abs(x)+abs(y);
%v_norm_2(i) = norma_2 del vettore [x(i);y(i)]
v_norm_2   = sqrt(x.*x+y.*y);
%v_norm_inf(i)=norma_infinito del vettore [x(i);y(i)]
v_norm_inf = max( [abs(x);abs(y)] );

disp('Ora calcolo i corrispondenti versori')
disp('Strike any key to continue')
pause

x_1 = x./v_norm_1;
y_1 = y./v_norm_1;

x_2 = x./v_norm_2;
y_2 = y./v_norm_2;

x_inf = x./v_norm_inf;
y_inf = y./v_norm_inf;

disp('Ora rappresento le sfere unitarie nelle norme 1, 2 ed infinito')
disp('Strike any key to continue')
pause
figure
plot(x_1,y_1),grid('on'),axis('square'),axis([-2 2 -2 2])
title('Sfera unitaria in norma 1')

disp('Strike any key to continue')
pause
figure
plot(x_2,y_2),grid('on'),axis('square'),axis([-2 2 -2 2])
title('Sfera unitaria in norma 2')

disp('Strike any key to continue')
pause
figure
plot(x_inf,y_inf),grid('on'),axis('square'),axis([-2 2 -2 2])
title('Sfera unitaria in norma infinito')

disp('Ora calcolo le norme di matrice indotte')
disp('dalle precedenti norme di vettore')
disp('Strike any key to continue')
pause

disp('Genero una matrice 2x2 random')
pause
A = rand(2)

disp('Norma della matrice indotta dalla norma 1 di vettore')
disp('Strike any key to continue')
pause
A_1 = A*[x_1;y_1];

disp('Calcolo il valore approssimato della norma')
disp('il massimo e` fatto su un insieme finito di versori,')
disp('non su tutti i versori')
pause
A_norm_2_approx = max(abs(A_1(1,:))+abs(A_1(2,:)))

disp('Calcolo il valore "esatto" della norma')
pause
A_norm_2 = norm(A,1)

disp('Norma della matrice indotta dalla norma 2 di vettore')
disp('Strike any key to continue')
pause
A_2 = A*[x_2;y_2];
A_2 = A_2.*A_2;

disp('Valore approssimato')
pause
A_norm_2_approx = max(sqrt(A_2(1,:)+A_2(2,:)))

disp('Valore "esatto"')
pause
A_norm_2 = norm(A,2)

disp('Norma della matrice indotta dalla norma infinito di vettore')
disp('Strike any key to continue')
pause

disp('Valore approssimato');
pause
A_norm_inf_approx = max(max(A*[x_inf;y_inf]))

disp('Valore "esatto"')
pause
A_norm_inf = norm(A,inf)


/*************************************************************************/
/* Programma per il calcolo delle radici dell'equazione di secondo grado */
/*************************************************************************/

#include<stdio.h>
#include<math.h>

int main(void)
{
  int i;
  double a,b,c,x1,x2;
  float fla,flb,flc,flx1,flx2;
  a = 1.0;
  c = 2.0;

  fla = a;
  flc = c;
  for(i=1;i<10;i++)
    {
      b = 1.234*pow(10.0,2*i);
      flb = b;

      x1 = (-b-b/fabs(b)*sqrt(b*b-4.0*a*c))/(2.0*a);
      x2 = c/(a*x1);

      if(b/fabs(b)<0.0)
	{
	  flx1 = (-flb+sqrt(flb*flb-4.0*fla*flc))/(2.0*fla);
	  flx2 = (-flb-sqrt(flb*flb-4.0*fla*flc))/(2.0*fla);
	}
      else
	{
	  flx2 = (-flb+sqrt(flb*flb-4.0*fla*flc))/(2.0*fla);
	  flx1 = (-flb-sqrt(flb*flb-4.0*fla*flc))/(2.0*fla);
	}
      
      printf("b=%.6E\n",b);
      printf("  x1=%.6E\tflx1=%.6E\tabs err=%.6E\trel err=%.6E\n",
	     x1,flx1,fabs(x1-flx1),fabs(x1-flx1)/fabs(x1));
      printf("  x2=%.6E\tflx2=%.6E\tabs err=%.6E\trel err=%.6E\n",
	     x2,flx2,fabs(x2-flx2),fabs(x2-flx2)/fabs(x2));
    }
}