accuracy of intrinsic SUM function

[Schorsch MCMLX]

Dear gFortran-fans,

when comparing the results that different free F compilers yield for 
the intrinsic SUM function, I noticed that only the Salford FTN95 
compiler seems to use a decent summation algorithm. Using a sample of 
about 28k pseuso-random numbers drawn from a ND(48;4), the gFortran 
compiler V4.7.1 running on Win7x64, the gFortran compiler V4.6.1 
running on UbLi11x64, as well as the Intel iFort compiler V 12.1.2 
running on UbLi11x64 yield an error pattern that (nearly) exactly 
matches that of recurrent summation, i.e., starting with sum=value(1) 
and then iterate through all the other values by sum = sum + value(i).
The reference values were created by three independent multi-precision
tools that yielded always identical results. I always compared the 
binary representation of real*8 variables in order to avoid problems 
with possibly too few decimal digits.
The FNT95 error pattern exactly matches that of pairwise summation. 
For a description of that algorithm, see e.g. wikipedia: 
http://en.wikipedia.org/wiki/Pairwise_summation  . The gain in 
accuracy resulting from using the pairwise algo instead of the 
recurrent one is quite remarkable, as can be seen from the following 
table (max. error in ULPs). 

N=28404                    errors of different summation algorithms, 
in ULPs
Comp/OS		intrinsic	pairwise	recurrent	inline
(n_max=32)


---------------------------------------------------------------------------------------------
iFort/UbLi  min	-19	-2	-20	-3
iFort/UbLi  max	4	2	3	1


---------------------------------------------------------------------------------------------
gFort/UbLi min	-34	-2	-34	-1
gFort/UbLi max	31	2	31	1


---------------------------------------------------------------------------------------------
gFort/W7 min	-34	-2	-34	-1
gFort/W7 max	31	2	31	1


---------------------------------------------------------------------------------------------
FTN95/W7 min	-1	-2	-35	0
FTN95/W7 max	1	2	30	0

The algo referred to as inline in this table uses the pairwise algo, 
but having the formulae unrolled to a single statement for n<=32 in 
order to exploit the higher accuracy of the processor's internal 
registers.
I've not used any compiler switches in these tests, but I've checked 
that there are no switches set in FTN95 that are likely to influence 
the accuracy.
Concluding, I'd like to propose to implement the pairwise summation 
algorithm for the intrinsic SUM function in gFortran, if possible.

All the best
Schorsch