Why does squaring generate a different result?

The following example shows some effect of progressive underflow when compiled with Gfortran, although fails abruptly.

real    :: a,b,aa,bb
integer :: i

  a = 0.1
  b = 1.e17
  do i = 1,6
    write (*,*) i
    aa = a/b
    bb = a*b
    write ( *,*) (aa*aa)*bb, aa,bb
    write ( *,*) aa*(aa*bb), aa,bb
    b = b*10.
  end do
  end

Results with Gfortran are:

           1
   1.00000013E-20   1.00000005E-18   1.00000003E+16
   1.00000013E-20   1.00000005E-18   1.00000003E+16
           2
   1.00000007E-21   1.00000003E-19   9.99999984E+16
   9.99999968E-22   1.00000003E-19   9.99999984E+16
           3
   9.99994604E-23   1.00000005E-20   9.99999984E+17
   1.00000009E-22   1.00000005E-20   9.99999984E+17
           4
   1.00052712E-23   9.99999968E-22   9.99999998E+18
   9.99999921E-24   9.99999968E-22   9.99999998E+18
           5
   9.80908895E-25   9.99999968E-23   1.00000002E+20    ! significant error
   9.99999921E-25   9.99999968E-23   1.00000002E+20
           6
   0.00000000       1.00000000E-23   1.00000002E+21    ! total error
   9.99999958E-26   1.00000000E-23   1.00000002E+21
1 Like