! Questo programma calcola l'espressione a*a-b*b con due formule
! differenti e mette in luce il diverso errore algoritmico generato.
! Provare ad assegnare valori di a e b con molte cifre in comune,
! ad esempio a=1.23456 b=1.23457, e valori molto diversi ad esempio
! a=12.3456 e 0.123456

program primo_programma
  implicit none
  real :: a, b, c, s1, s2           ! floating point in semplice precisione
  real(kind(0.d0)) :: da, db, dc    ! floating point in doppia precisione

  write(*,*) "assegna il valore di a"
  read(*,*) a
  write(*,*) "assegna il valore di b"
  read(*,*) b
 
  ! calcolo a*a-b*b in doppia precisione
  da=a
  db=b
  dc=(da-db)*(da+db)
  
  ! calcolo a*a-b*b, scrivo risultato ed errore relativo
  s1=a*a
  s2=b*b
  c=s1-s2
  write(*,*) "a*a-b*b=",c, "err=",(c-dc)/dc
  
  ! calcolo (a-b)*(a+b), scrivo risultato ed errore relativo
  s1=a-b
  s2=a+b
  c=s1*s2
  write(*,*) "(a-b)*(a+b)=",c,(c-dc)/dc
end program primo_programma

