title: How are these double precision values accurate to 20 decimals?
body: <p>I am testing some very simple equivalence errors when precision is an issue and was hoping to perform the operations in extended double precision (so that I knew what the answer would be in ~19 digits) and then perform the same operations in double precision (where there would be roundoff error in the 16th digit), but somehow my double precision arithmetic is maintaining 19 digits of accuracy.</p>
  
  <p>When I perform the operations in extended double, then hardcode the numbers into another Fortran routine, I get the expected errors, but is there something strange going on when I assign an extended double precision variable to a double precision variable here?</p>
  
  <pre><code>program code_gen
      implicit none 
      integer, parameter :: Edp = selected_real_kind(17)
      integer, parameter :: dp = selected_real_kind(8)
      real(kind=Edp) :: alpha10, x10, y10, z10 
      real(kind=dp) :: alpha8, x8, y8, z8
  
      real(kind = dp) :: pi_dp = 3.1415926535897932384626433832795028841971693993751058209749445
  
      integer :: iter
      integer :: niters = 10
  
      print*, 'tiny(x10) = ', tiny(x10)
      print*, 'tiny(x8)  = ', tiny(x8)
      print*, 'epsilon(x10) = ', epsilon(x10)
      print*, 'epsilon(x8)  = ', epsilon(x8)
  
      do iter = 1,niters
          x10 = rand()
          y10 = rand()
          z10 = rand()
          alpha10 = x10*(y10+z10)
  
          x8 = x10 
          x8 = x8 - pi_dp
          x8 = x8 + pi_dp
          y8 = y10 
          y8 = y8 - pi_dp
          y8 = y8 + pi_dp
          z8 = z10 
          z8 = z8 - pi_dp
          z8 = z8 + pi_dp
          alpha8 = alpha10
  
          write(*, '(a, es30.20)') 'alpha8 .... ', x8*(y8+z8)
          write(*, '(a, es30.20)') 'alpha10 ... ', alpha10
  
          if( alpha8 .gt. x8*(y8+z8) ) then
              write(*, '(a)') 'ERROR(.gt.)'
          elseif( alpha8 .lt. x8*(y8+z8) ) then
              write(*, '(a)') 'ERROR(.lt.)'
          endif
      enddo
  end program code_gen
  </code></pre>
  
  <p>where <code>rand()</code> is the gfortran function found <a href="https://gcc.gnu.org/onlinedocs/gfortran/RAND.html" rel="nofollow">here</a>.</p>
  
  <p>If we are speaking about only one precision type (take, for example, double), then we can denote machine epsilon as <code>E16</code> which is approximately <code>2.22E-16</code>. If we take a simple addition of two Real numbers, <code>x+y</code>, then the resulting machine expressed number is <code>(x+y)*(1+d1)</code> where <code>abs(d1) &lt; E16</code>. Likewise, if we then multiply that number by <code>z</code>, the resulting value is really <code>(z*((x+y)*(1+d1))*(1+d2))</code> which is nearly <code>(z*(x+y)*(1+d1+d2))</code> where <code>abs(d1+d2) &lt; 2*E16</code>. If we now move to extended double precision, then the only thing that changes is that <code>E16</code> turns to <code>E20</code> and has a value of around <code>1.08E-19</code>.</p>
  
  <p>My hope was to perform the analysis in extended double precision so that I could compare two numbers which <em>should</em> be equal but show that, on occasion, roundoff error will cause comparisons to fail. By assigning <code>x8=x10</code>, I was hoping to create a double precision 'version' of the extended double precision value <code>x10</code>, where only the first ~16 digits of <code>x8</code> conform to the values of <code>x10</code>, but upon printing out the values, it shows that all 20 digits are the same and the expected double precision roundoff error is not occurring as I would expect.</p>
  
  <p>It should also be noted that before this attempt, I wrote a program which actually writes <em>another</em> program where the values of <code>x</code>, <code>y</code>, and <code>z</code> are 'hardcoded' to 20 decimal places. In this version of the program, the comparisons of <code>.gt.</code> and <code>.lt.</code> failed as expected, but I am not able to duplicate the same failures by casting an extended double precision value as a double precision variable.</p>
  
  <p>In an attempt to further 'perturb' the double precision values and add roundoff error, I have added, then substracted, <code>pi</code> from my double precision variables which should leave the remaining variables with some double precision roundoff error, but I am still not seeing that in the final result.</p>