text: Paul Murphy asked people to run the following code using default
  compilers and no special switches, and to report the result:
  % cc -o sqr -O sqr.c -lm
  % time ./sqr
  #include <stdio.h>
  #include <math.h>
  main()
  {
  register i;
  double f=0.0;
  double sqrt();
  for(i=0;i<=1000000000;i++)
  {
  f+=sqrt( (double) i);
  }
  printf("Finished %20.14f\n",f);
  }
  The table at the end of the article summarizes some reported
  results. Studying these suggests the questions: why are the results
  different for different (IEEE754 compliant) computers and different
  compilers, what is the "correct" result, how is it found, why are some results apparently more accurate than others, and why are the results on Sun, IBM, and other RISC CPU based computers about the same?
  Murphy discovered that one method of computing a "correct" result was the utilization of the Euler-Maclaurin formula and high precision computation. It gives 21081851083600.37596259382529338 as a good approximation to the correct result.
  21081851083600.37596259382529338
  The main points to be made are:
  Next is a more detailed discussion. For a good foundation in
  floating-point computing see the ACM article by David Goldberg "What Every Computer Scientist Should Know about Floating-Point Computing" especially with the addendum by Douglas Priest that addresses IA32. Also of interest might be Joe Darcy's JavaOne talk on What Everybody Using the Java Programming Language Should Know About Floating-Point Arithmetic.
  The table presents results from some of the 2, shows a result that in decimal representation, is in
  error by -.00096259382529338. However, observe that the hex
  representation shows an error of 0 ULPs, an exact match..
  Consider now the errors encountered in the survey, based on the previous table. Consider #1 as perfect. Errors can be considered in terms of decimal representation, and ULPs.
  When considering the absolute error alone, SPARC results appear very inaccurate. However, considering the ratio error/correct, the variance between results, in fact, is discovered to be very small.
  There are two sources of error in the algorithm as provided in
  Murphy's article. The first is due to the sqrt operation, and the
  second is due to the addition in summing the square roots.
  Sqrt rounding error: the largest rounding errors for the sqrt
  operations should occur for the largest numbers. A one ULP error for
  sqrt(1,000,000,000) is on the order of 2*10^(-12). This rounding error
  is very small compared to the one that is possible due to the addition
  that follows the sqrt operations.
  The sum just before sqrt(1000000000) is 21081851051977.5977. A 1
  ULP error in the sum value results in an error of about .0029. Huge
  compared to the rounding error of the sqrt operation.
  The SPARC-quad run shows that a lot more precision will yield the
  correct result. The cost as Murphy points out, is significantly more
  compute time.
  The question is, can something be done to minimize the inaccuracy
  without pushing precision to quad?
  Yes.
  The following algorithm uses the notion of compensated summation, which
  deals with the rounding error that is associated with each add in the
  summation.
  #include <stdio.h>
  #include <math.h>
  main()
  {
  register i;
  double sum=0.0,newsum,f;
  double sumerr=0.0;
  double err;
  double sqrt();
  for(i=1;i<=1000000000;i++)
  {
  f=sqrt( (double) i);
  newsum = sum + f;
  err = (newsum - sum ) - f;
  sumerr += err;
  sum = newsum;
  }
  printf("Finished %20.14f\n",sum);
  printf("Finished %20.14f\n",sum-sumerr);
  }
  This program was compiled/linked with the same command, but
  returned results matching Correct while only taking about 20% more compute time.
  The entry "SPARC double+e" in the following table shows the result
  using compensated summation.
  We offer a bit more detail about compensated summation in the Notes
  section at the end of this article.
  There is a clustering of systems in terms of errors, identified as RISC
  machines that adhere to IEEE-754 64-bit format for double precision
  work.
  The best results were from systems based on Intel's Pentium model.
  Pentium systems have extended precision registers, 80 bits versus
  SPARC's 64-bits in double. In addition, Pentium systems allow
  intermediate results to be kept in extended format. An explicit call is
  required to force the Pentium processor to keep intermediate results in
  64-bit format.
  As the preceding section shows, the critical issue of concern here
  should be the precision used in the problem. As Murphy points out, the
  correct value was calculated using 200 digits of precision.
  While the results from the Pentium class machines were better, the
  question arguably remains: Are any of these answers bad? We contend
  that the answer to this question is no. Given the magnitude of the
  number, and the size of the error, the results are all perfectly
  reasonable.
  Finally, what should be learned from this example is that a
  precision level may work for one algorithm and data set, but fail if
  the data changes or if the algorithm is modified.
  Another interpretation of this is that the compilers of C don't have
  strong guarantees about what "double" means. In some versions of gcc
  with some compiler options, "double" can mean "80-bits" instead of
  64-bits. Some of the gcc compilers might be new enough to have some C99
  support so it is conceivable that double_t is being used instead of
  double.
  If one wants more accurate results without compensated summation,
  one should declare f as "long double" or "double_t". If one wants
  predictable results one can use Java, even on x86.
  The one pertinent observation in Murphy's analysis of his example is
  that the results on SPARC systems dating back to 1989 are all the same,
  whereas the results on x86 systems vary in ways that are not at all
  obvious.
  In the original algorithm:
  for(i=0;i<=1000000000;i++)
  {
  f+=sqrt( (double) i);
  }
  whenever the addition is done, there are two possibilities:
  With each addition, the rounding error is ignored. After a billion
  additions, the results are noticeable in this problem.
  The compensating algorithm does the following:
  newsum = sum + f;
  (newsum - sum )
  err = (newsum - sum ) - f;
  sumerr += err;
  Even the simple example of square-root summation affords the
  opportunity for appreciation of the complexity of floating-point
  computation. The complexity is inherent in the attempt to model
  arithmetic with real numbers, which in "true" arithmetic have infinite
  precision, by a finite precision approximation. Almost every arithmetic
  operation, therefore, results in an error and, in fact, it is in
  addition/subtraction that the largest errors can occur. Thus
  floating-point computation is epitomized by the need to understand and
  estimate the error to ensure the validity of what one computes.
  About the Authors
  Gregory Tarsy is the manager of the Floating Point and Numerical
  Computing group at Sun Microsystems. Prior to joining Sun he worked
  in scientific
  applications on supercomputers and was a professor of
  mathematics at the University of California and City University of
  New York.
  Neil Toda is a member of Sun's Floating Point and Numerical
  Computing
  group. Prior to joining Sun, he worked in Operating Systems
  Engineering,
  Database Engineering, and Software Performance Engineering.