Puzzle computinig arrays

[Tim Prince]

Angelo Graziosi wrote:
> The following happens on Cygwin with GFortran 4.3.1 and on GNU/Linux 
> Kubuntu 8.04 with GFortran 4.3.1 and 4.2.3.
>
> In the following test case (reduced to the simplest form), the array 
> elements of BB() and b() are computed exactly with the same expression 
> (copy/pasted), but the results seem a little different (only on 
> 14-15th digit)
>
> $ cat test_case.f90
> !
> ! gfortran test_case.f90 -o test_case
> !
> program test_case
>   implicit none
>   integer :: k
>   integer, parameter :: DP = kind(1.D0),&
>        N = 3
>   real(DP), parameter :: A(0:N) = &
>        (/1.0000000000D0,3.8860009363D0,7.4167881843D0,&
>          9.8599837415D0/)
>   real(DP), parameter :: BB(0:N-1) = &
>        (/((A(k+1)*(A(1)+k+1)),k = 0,N-2),(A(N)*(A(1)+N))/)
>   real(DP) :: b(0:N-1)
>   do k = 0,N-2
>      b(k) = A(k+1)*(A(1)+k+1)
>   enddo
>   b(N-1) = A(N)*(A(1)+N)
>   do k = 0,N-1
>      !print *, k,b(k),BB(k),b(k)-BB(k)
>      print *, k,b(k)-BB(k)
>   enddo
> end program test_case
>
> $ gfortran test_case.f90 -o test_case
>
> $ ./test_case.exe
>            0   0.0000000000000000
>            1 -7.10542735760100186E-015
>            2   0.0000000000000000
>
> Since BB() and b() are computed with the same expression, I would 
> expect the same results. But BB() is a PARAMETER (computed at compile 
> time???) and b() is computed at run time...
>
> So, what could be the difference? What explanation may you give?
>
Don't the compile time calculations get extra precision from mpfr and 
friends?  The run-time calculation would make more sense with 
parentheses (k+1) to specify a single round-off.