4.1 A Summing Mystery?athe Magnitude Problem
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 10 is The program computes the fraction sums with denominators s 10 , s 100 , s 1000 , 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 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. 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 |
at the end of Chapter 3. Since the sum is supposed to be 1, and the float 
 |
| Top |
EAN: 2147483647
Pages: 141