%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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')


/*************************************************************************/
/* 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));
    }
}