Standards

Digg This!

Page 2 of 2   « previous page

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 storage to quad?

Yes.

III) Algorithm Correcting for Errors

The following algorithm uses the notion of compensated summation, which deals with the rounding error that is associated with each add in the summation.

        ------- sqrX.c -------
        #include 
        #include 
        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);
        }
        ------- sqrX.c -------

1) Correct	: 21081851083600.37596259382529338	=0x42b32c803ebb5060
!) SPARC double+e	: 21081851083600.37500000000000	=0x42b32c803ebb5060
2) SPARC-quad	: 21081851083600.37500000000000	=0x42b32c803ebb5060
3) Correct+1ULP	: 21081851083600.3789	=0x42b32c803ebb5061
4) BEST	: 21081851083600.38281250000000	=0x42b32c803ebb5062
5) SPARC double	: 21081851083600.55859375000000	=0x42b32c803ebb508f

This program was compiled/linked with the same command, but returned results matching Correct while only taking about 20% more compute time.

In the NOTES section at the end of this paper, we offer a bit more detail about compensated summation.

IV) Why are some results better than SPARC?

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.

V) IEEE-754 and Sun's Results

IEEE 754 aids understanding where errors arise in computations by enabling us to work out what error is committed in each operation just by knowing the operands and the destination format; we need not know such system-dependent details as whether the hardware uses a guard bit or whether it chops or rounds. But IEEE 754 also obscures understanding why a particular system computes the result that it does because it allows implementations to arbitrarily choose which operations are rounded to the format the programmer specifies and which are rounded to a wider format.

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.

NOTES

A) Compensated Summation

        ------- sqrX.c -------
        #include 
        #include 
        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);
        }
        ------- sqrX.c -------

In the original algorithm:

        for(i=0;i<=1000000000;i++)
        {
        f+=sqrt( (double) i);
        }

whenever the addition is done, there are two possibilities:

1. The result is exact and fits into the 64-bits allotted, or
2. the result will not exactly fit, and rounding needs to occur.

With each addition, the rounding error is ignored. After a billion additions, the results are noticeable in this problem.

The compensating algorithm is doing the following:

• Do the add, don't worry about rounding error:
  newsum = sum + f;
• Determine the part of f that contributes to the new value:
  (newsum - sum )
• f's contribution to the new sum minus f give the amount of f that is lost, i.e. the rounding error:
  err = (newsum - sum ) - f;
• Keep track of the rounding errors separately:
  sumerr += err;
• All that is left is to adjust things so that the loop can be used again.

Page 2 of 2   « previous page