Chapter 4. Summing Lists of Numbers

	Java Number Cruncher: The Java Programmer's Guide to Numerical Computing By Ronald Mak
	Table of Contents

	Part II. Iterative Computations

Chapter 4. Summing Lists of Numbers

A typical assignment in Computer Programming 101 is to write a program that sums a list of numbers. If you add many numbers, say hundreds or thousands of them, and there are no overflows, then you might reasonably expect that the roundoff errors will cancel each other out, and the sum will be fairly accurate. What could possibly go wrong?

This chapter shows that even such a simple operation can go wrong if we're careless about floating-point arithmetic. Several summation problems will illustrate some of the typical pitfalls. We'll take an experimental approach by writing programs that try different summation techniques to find ones that either compensate for the roundoff errors or reduce the opportunities for the errors to occur. (After all, isn't this computer science? )

Top

	Java Number Cruncher: The Java Programmer's Guide to Numerical Computing By Ronald Mak
	Table of Contents

	Chapter 4. Summing Lists of Numbers

4.1 A Summing Mystery?athe Magnitude Problem

Let's start with a simple summing problem. The formula for the sum s _n of the first n counting numbers, 1 through n, is

For example, the sum of the first 10 counting numbers is 55.

Program 4-1 adds a slight twist to the summation. See Listing 4-1. It adds fractions instead, where, for each summation, the numerators are the first n counting numbers and the denominator is the sum s _n . Therefore, the sum should be 1. For example, the fraction sum for s ₁₀ is

The program computes the fraction sums with denominators s ₁₀ , s ₁₀₀ , s ₁₀₀₀ , and so on through s _100,000,000 .

We can see from the output that things go fairly smoothly, with the computed sums close to 1, until denominator s _100,000,000 ?aand disaster strikes. In the last line of output, how did the sum manage to be 0.5? That's not even close!

Listing 4-1 Fraction sums.

package numbercruncher.program4_1;

import numbercruncher.mathutils.AlignRight;

/**
 * PROGRAM 4-1: Fraction Sums
 *
 * For each n, compute the sum 1/d + 2/d + 3/d + ... + n/d = d/d
 * where:
 *     n = 10, 100, 1000, ..., 100,000,000
 *     d = 1 + 2 + 3 + ... + n = (n/2)(n + 1)
 *
 * See how close each sum is to 1.
 */
public class FractionSums
{
    public static void main(String args[])
    {
        AlignRight ar = new AlignRight();

        ar.print("n", 9); ar.print("Denom", 14);
        ar.print("1/Denom", 14); ar.print("n/Denom", 13);
        ar.print("Sum", 11); ar.print("% Error", 13);
        ar.underline();

        for (int n = 10; n <= 100000000; n *= 10) {
            float sum   = 0;
            float denom = (0.5f*n)*(n + 1);

            // Sum fractions.
            for (int i = 1; i <= n; ++i) sum += i/denom;

            ar.print(n, 9); ar.print(denom, 14); ar.print(1/denom, 14);
            ar.print(n/denom, 13); ar.print(sum, 11);
            ar.print(100*Math.abs(sum - 1), 13);
            ar.println();
        }
    }
}

Output:

n         Denom       1/Denom      n/Denom        Sum      % Error
---------------------------------------------------------------------------
       10          55.0   0.018181818   0.18181819 0.99999994 5.9604645E-6
      100        5050.0   1.980198E-4   0.01980198        1.0          0.0
     1000      500500.0   1.998002E-6  0.001998002 0.99999994 5.9604645E-6
    10000      5.0005E7     1.9998E-8    1.9998E-4  0.9999999 1.1920929E-5
   100000   5.0000502E9 1.9999799E-10   1.99998E-5 0.99999875 1.2516975E-4
  1000000 5.00000489E11  1.999998E-12  1.999998E-6  0.9998997  0.010031462
 10000000  5.0000004E13 1.9999998E-14 1.9999999E-7   1.002663   0.26630163
100000000  5.0000001E15 1.9999999E-16       2.0E-8        0.5         50.0

Why should the fraction sum for denominator s _100,000,000 fail so miserably when the previous sums all did so well? It isn't underflow or overflow. For each sum, the program prints out the values of the first fraction and last fraction, and their exponent values are well within the range (-45 through +38) for the float type. Each sum is computed from many addends, and so we might expect most of the high and low roundoff errors to cancel each other out.

For one possible source of the problem, recall our computation of the machine epsilon at the end of Chapter 3. Since the sum is supposed to be 1, and the float value is around 6.0 x 10 ^- ⁸ , we can see that, with denominator s _100,000,000 , a significant percentage of the addends may have no effect on the sum as it approaches 1.

Program 4-2 provides some insight. See Listing 4-2. It computes the fraction sum for denominator s _100,000,000 by dividing the series of fractions into 20 equal- sized groups, and at the end of each group it prints the current fraction and the current value of the running sum. Before adding each fraction to the running sum, it calculates the difference between the fraction's binary exponent and the running sum's binary exponent. At the end of each group, it also prints the minimum exponent difference of that group and the percentage of exponent differences within the group that exceeded 24. (Recall that, at the end of Chapter 3, we saw that, for type float, the unbiased binary exponent of the machine epsilon is -24.)

Listing 4-2 Fraction sums for n = 100,000,000.

package numbercruncher.program4_2;

import numbercruncher.mathutils.AlignRight;

/**
 * PROGRAM 4-2: Fraction Sum 100M
 *
 * Compute the sum 1/d + 2/d + 3/d + ... + n/d = d/d
 * where:
 *     n = 100,000,000
 *     d = 1 + 2 + 3 + ... + n = (n/2)(n + 1)
 *
 * See why the sum ends up being 0.5
 */
public class FractionSum100M
{
    private static final int GROUPS = 20;
    private static final int MAX    = 100000000;    // 100M

    public static void main(String args[])
    {
        AlignRight ar = new AlignRight();

        ar.print("i", 9);
        ar.print("Fraction", 15);
        ar.print("Running sum", 15);
        ar.print("Min exp diff", 15);
        ar.print("% ExpDiff>24", 15);
        ar.underline();

        float sum   = 0;
        float denom = (0.5f*MAX)*(MAX + 1);

        int gSize   = MAX/GROUPS;         // group size
        int gEnd    = gSize;              // index of group end
        int minDiff = Integer.MAX_VALUE;  // min exponent difference
        int exceeds = 0;                  // # of exponent diff > 24

        // Loop to sum the fractions.
        for (int i = 1; i <= MAX; ++i) {
            float fraction = i/denom;

            // Compute the exponent difference.
            int expSum      = Float.floatToIntBits(sum)      >> 23;
            int expFraction = Float.floatToIntBits(fraction) >> 23;
            int diff        = Math.abs(expSum - expFraction);

            if (i > 1) {
                minDiff = Math.min(minDiff, diff);
                if (diff > 24) ++exceeds;
            }

            sum += fraction;

            // Printout at the end of each group.
            if (i == gEnd) {
                ar.print(i, 9);
                ar.print(fraction, 15);
                ar.print(sum, 15);
                ar.print(minDiff, 15);
                ar.print((100*exceeds)/gSize, 15);
                ar.println();

                minDiff = Integer.MAX_VALUE;
                exceeds = 0;
                gEnd += gSize;
            }
        }
    }
}

Output:

i       Fraction    Running sum   Min exp diff   % ExpDiff>24
---------------------------------------------------------------------
  5000000         1.0E-9    0.002493547              0              0
 10000000         2.0E-9    0.009991142             21              0
 15000000         3.0E-9    0.022081949             22              0
 20000000         4.0E-9      0.0407084             23              0
 25000000         5.0E-9     0.05933485             23              0
 30000000         6.0E-9     0.09342261             23              0
 35000000   6.9999997E-9          0.125             24             15
 40000000         8.0E-9     0.16593489             24             45
 45000000         9.0E-9      0.2404407             24              0
 50000000         1.0E-8           0.25             24             87
 55000000         1.1E-8           0.25             25            100
 60000000         1.2E-8           0.25             25            100
 65000000   1.2999999E-8           0.25             25            100
 70000000   1.3999999E-8           0.25             25            100
 75000000         1.5E-8     0.26472795             24             90
 80000000         1.6E-8     0.41373956             24              0
 85000000         1.7E-8            0.5             24             42
 90000000         1.8E-8            0.5             25            100
 95000000         1.9E-8            0.5             25            100
100000000         2.0E-8            0.5             25            100

Now we begin to see what happened . If the difference in magnitude between two addends is too great, the smaller addend loses much of its contribution to the sum. Indeed, as we saw at the end of Chapter 3, smaller addends will have no effect at all on the sum if the difference between their binary exponents is greater than 24 when doing float addition.

Program 4-2 shows that, whenever the value of the running sum reaches the next power of 2, the sum's binary exponent increases by 1. If that widens the exponent difference between the sum and the upcoming fractions to more than 24, then those fractions have no effect on the sum until they become large enough to lower the exponent difference. Thus, the running sum "stalls" at 0.25 and again at 0.5 when the minimum exponent differences are 25. Figure 4-1 shows this clearly by graphing the values of the running sum versus the minimum exponent differences in each group.

Figure 4-1. A graph of the output of Program 4 §C2. The running sum "stalls" at 0.25 and at 0.5 when the minimum exponent differences exceed 24.

graphics/04fig01.gif

Would it make a difference if we summed the fractions in the reverse order, from the largest down to the smallest? Program 4-3 shows that the running sum still stalls at 0.5, and we have even greater percentages of exponent differences exceeding 24. See Listing 4-3. So we've made the problem worse ?athe values of the running sum increase while the fractions decrease, making the exponent differences even larger.

In general, when adding a sorted list of numbers all of the same sign, it is better to add them in the order from the smallest magnitude up to the largest.

Listing 4-3 Fraction sums for n = 100,000,000 in reversed order.

package numbercruncher.program4_3;

import numbercruncher.mathutils.AlignRight;

/**
 * PROGRAM 4-3: Fraction Sum 100M in the Reversed Order
 *
 * Compute the sum n/d + (n-1)/d + (n-2)/d + ... + 2/d + 1/d = d/d
 * where:
 *     n = 100,000,000
 *     d = 1 + 2 + 3 + ... + n = (n/2)(n + 1)
 *
 * See if the sum is closer to 1.
 */
public class FractionSum100MReversed
{
    private static final int GROUPS = 20;
    private static final int MAX    = 100000000;    // 100M

    public static void main(String args[])
    {
        AlignRight ar = new AlignRight();

        ar.print("i", 9);
        ar.print("Fraction", 15);
        ar.print("Running sum", 15);
        ar.print("Min exp diff", 15);
        ar.print("% ExpDiff>24", 15);
        ar.underline();

        float sum   = 0;
        float denom = (0.5f*MAX)*(MAX + 1);

        int gSize   = MAX/GROUPS;         // group size
        int gStart  = (MAX + 1) - gSize;  // index of group start
        int minDiff = Integer.MAX_VALUE;  // min exponent difference
        int exceeds = 0;                  // # of exponent diff > 24

        // Loop to sum the fractions.
        for (int i = MAX; i >= 1; - i) {
            float fraction = i/denom;

            int expSum      = Float.floatToIntBits(sum)      >> 23;
            int expFraction = Float.floatToIntBits(fraction) >> 23;
            int diff        = Math.abs(expSum - expFraction);

            if (i < MAX) {
                minDiff = Math.min(minDiff, diff);
                if (diff > 24) ++exceeds;
            }

            sum += fraction;

            // Printout at the start of each group.
            if (i == gStart) {
                ar.print(i, 9);
                ar.print(fraction, 15);
                ar.print(sum, 15);
                ar.print(minDiff, 15);
                ar.print((100*exceeds)/gSize, 15);
                ar.println();

                minDiff = Integer.MAX_VALUE;
                exceeds = 0;
                gStart -= gSize;
            }
        }
    }
}

Output:

i       Fraction    Running sum   Min exp diff   % ExpDiff>24
---------------------------------------------------------------------
 95000001         1.9E-8     0.10206167              0              0
 90000001         1.8E-8     0.18421358             22              0
 85000001         1.7E-8     0.26743877             23              0
 80000001         1.6E-8     0.41645038             24              0
 75000001         1.5E-8            0.5             24             43
 70000001   1.3999999E-8            0.5             25            100
 65000001   1.2999999E-8            0.5             26            100
 60000001         1.2E-8            0.5             26            100
 55000001         1.1E-8            0.5             26            100
 50000001         1.0E-8            0.5             26            100
 45000001         9.0E-9            0.5             26            100
 40000001         8.0E-9            0.5             26            100
 35000001   6.9999997E-9            0.5             26            100
 30000001         6.0E-9            0.5             27            100
 25000001         5.0E-9            0.5             27            100
 20000001         4.0E-9            0.5             27            100
 15000001         3.0E-9            0.5             27            100
 10000001   2.0000002E-9            0.5             28            100
  5000001   1.0000002E-9            0.5             28            100
        1  1.9999999E-16            0.5             29            100

We have a clear challenge. In order for a summation to be accurate, we need to keep the running sum and each addend sufficiently close together in magnitude, so that the difference between their binary exponents is not greater than 24.

Program 4-4 shows one way to do this. See Listing 4-4. The program partitions the fractions into 20 equal-sized groups. It computes a separate subtotal for each group; then it adds the subtotal to the running sum. Each group's fractions are closer in magnitude to the group subtotal, so the minimum difference between the binary exponents of a fraction and its group subtotal is under 24. Each group subtotal has enough "heft" to make its contribution to the overall sum. The trick is to get the right size for the groups. The final sum is very close to 1.

Listing 4-4 Fraction sums for n = 100,000,000 using group subtotals.

package numbercruncher.program4_4;

import numbercruncher.mathutils.AlignRight;
/**
 * PROGRAM 4-4: Fraction Sum 100M Using Group Subtotals
 *
 * Compute the sum 1/d + 2/d + 3/d + ... + n/d = d/d
 * where:
 *     n = 100,000,000
 *     d = 1 + 2 + 3 + ... + n = (n/2)(n + 1)
 *
 * Compute the sum by adding the group subtotals.
 * See if the final sum is closer to 1.
 */
public class FractionSum100MGrouped
{
    private static final int GROUPS = 20;
    private static final int MAX    = 100000000;    // 100M

    public static void main(String args[])
    {
        AlignRight ar = new AlignRight();

        ar.print("i", 9);
        ar.print("Group sum", 15);
        ar.print("Running sum", 15);
        ar.print("% ExpDiff>24", 15);
        ar.underline();

        float sum      = 0;           // running sum
        float subtotal = 0;           // group subtotal
        int   gSize    = MAX/GROUPS;  // group size
        int   gEnd     = gSize;       // index of group end
        int   exceeds  = 0;           // # of exponent diff > 24

        float denom = (0.5f*MAX)*(MAX + 1);

        // Sum the fractions by groups.
        for (int i = 1; i <= MAX; ++i) {
            float fraction = i/denom;

            int expSubtotal = Float.floatToIntBits(subtotal) >> 23;
            int expFraction = Float.floatToIntBits(fraction) >> 23;
            int diff        = Math.abs(expSubtotal - expFraction);

            if ((subtotal > 0) && (diff > 24)) ++exceeds;
            subtotal += fraction;

            // Subtotal and printout at the end of each group.
            if (i == gEnd) {
                sum += subtotal;

                ar.print(i, 9);
                ar.print(subtotal, 15);
                ar.print(sum, 15);
                ar.print((100*exceeds)/gSize, 15);
                ar.println();

                subtotal = 0;
                exceeds  = 0;
                gEnd += gSize;
            }
        }
        System.out.println("\n% error = " + 100*Math.abs(sum - 1));
    }
}

Output:

i      Group sum    Running sum   % ExpDiff>24
------------------------------------------------------
  5000000    0.002493547    0.002493547              0
 10000000    0.007556428   0.0100499755              0
 15000000    0.012413543     0.02246352              0
 20000000    0.017349107    0.039812624              0
 25000000    0.023679595     0.06349222              0
 30000000     0.02718998      0.0906822              0
 35000000    0.034128636     0.12481084              0
 40000000    0.036770158       0.161581              0
 45000000     0.04185812     0.20343912              0
 50000000    0.049500026     0.25293913              0
 55000000     0.05435951     0.30729866              0
 60000000    0.055825584     0.36312425              0
 65000000     0.05929653     0.42242077              0
 70000000     0.07003953      0.4924603              0
 75000000     0.07346739      0.5659277              0
 80000000     0.07455075     0.64047843              0
 85000000     0.07784402     0.71832246              0
 90000000     0.08436947     0.80269194              0
 95000000    0.098553196      0.9012451              0
100000000     0.09986824      1.0011134              0

% error = 0.11134148

Top