--- old/src/share/classes/java/util/DoubleSummaryStatistics.java 2013-11-20 00:17:18.000000000 -0800 +++ new/src/share/classes/java/util/DoubleSummaryStatistics.java 2013-11-20 00:17:18.000000000 -0800 @@ -63,6 +63,7 @@ public class DoubleSummaryStatistics implements DoubleConsumer { private long count; private double sum; + private double sumCompensation; // Low order bits of sum private double min = Double.POSITIVE_INFINITY; private double max = Double.NEGATIVE_INFINITY; @@ -81,7 +82,7 @@ @Override public void accept(double value) { ++count; - sum += value; + sumWithCompensation(value); min = Math.min(min, value); max = Math.max(max, value); } @@ -95,12 +96,24 @@ */ public void combine(DoubleSummaryStatistics other) { count += other.count; - sum += other.sum; + sumWithCompensation(other.sum); + sumWithCompensation(other.sumCompensation); min = Math.min(min, other.min); max = Math.max(max, other.max); } /** + * Incorporate a new double value using Kahan summation / + * compensated summation. + */ + private void sumWithCompensation(double value) { + double tmp = value - sumCompensation; + double velvel = sum + tmp; // Little wolf of rounding error + sumCompensation = (velvel - sum) - tmp; + sum = velvel; + } + + /** * Return the count of values recorded. * * @return the count of values --- old/src/share/classes/java/util/stream/Collectors.java 2013-11-20 00:17:19.000000000 -0800 +++ new/src/share/classes/java/util/stream/Collectors.java 2013-11-20 00:17:18.000000000 -0800 @@ -506,12 +506,34 @@ public static Collector summingDouble(ToDoubleFunction mapper) { return new CollectorImpl<>( - () -> new double[1], - (a, t) -> { a[0] += mapper.applyAsDouble(t); }, - (a, b) -> { a[0] += b[0]; return a; }, - a -> a[0], CH_NOID); + () -> new double[2], + (a, t) -> { sumWithCompensation(a, mapper.applyAsDouble(t)); }, + (a, b) -> { sumWithCompensation(a, b[0]); return sumWithCompensation(a, b[1]); }, + a -> a[0], + CH_NOID); + } + + /** + * Incorporate a new double value using Kahan summation / + * compensation summation. + * + * High-order bits of the sum are in intermediateSum[0], low-order + * bits of the sum are in intermediateSum[1], any additional + * elements are application-specific. + * + * @param intermediateSum the high-order and low-order words of the intermediate sum + * @param value the name value to be included in the running sum + */ + static double[] sumWithCompensation(double[] intermediateSum, double value) { + double tmp = value - intermediateSum[1]; + double sum = intermediateSum[0]; + double velvel = sum + tmp; // Little wolf of rounding error + intermediateSum[1] = (velvel - sum) - tmp; + intermediateSum[0] = velvel; + return intermediateSum; } + /** * Returns a {@code Collector} that produces the arithmetic mean of an integer-valued * function applied to the input elements. If no elements are present, @@ -567,10 +589,11 @@ public static Collector averagingDouble(ToDoubleFunction mapper) { return new CollectorImpl<>( - () -> new double[2], - (a, t) -> { a[0] += mapper.applyAsDouble(t); a[1]++; }, - (a, b) -> { a[0] += b[0]; a[1] += b[1]; return a; }, - a -> (a[1] == 0) ? 0.0d : a[0] / a[1], CH_NOID); + () -> new double[3], + (a, t) -> { sumWithCompensation(a, mapper.applyAsDouble(t)); a[2]++; }, + (a, b) -> { sumWithCompensation(a, b[0]); sumWithCompensation(a, b[1]); a[2] += b[2]; return a; }, + a -> (a[2] == 0) ? 0.0d : (a[0] / a[2]), + CH_NOID); } /** --- old/src/share/classes/java/util/stream/DoublePipeline.java 2013-11-20 00:17:19.000000000 -0800 +++ new/src/share/classes/java/util/stream/DoublePipeline.java 2013-11-20 00:17:19.000000000 -0800 @@ -377,8 +377,15 @@ @Override public final double sum() { - // TODO: better algorithm to compensate for errors - return reduce(0.0, Double::sum); + double[] summation = collect(() -> new double[2], + (ll, d) -> { + Collectors.sumWithCompensation(ll, d); + }, + (ll, rr) -> { + Collectors.sumWithCompensation(ll, rr[0]); + Collectors.sumWithCompensation(ll, rr[1]); + }); + return summation[0]; } @Override @@ -391,19 +398,29 @@ return reduce(Math::max); } + /** + * {@inheritDoc} + * + * @implNote The {@code double} format can represent all + * consecutive integers in the range -253 to + * 253. If the pipeline has more than 253 + * values, the divisor in the average computation will saturate at + * 253, leading to additional numerical errors. + */ @Override public final OptionalDouble average() { - double[] avg = collect(() -> new double[2], - (ll, i) -> { - ll[0]++; - ll[1] += i; + double[] avg = collect(() -> new double[3], + (ll, d) -> { + ll[2]++; + Collectors.sumWithCompensation(ll, d); }, (ll, rr) -> { - ll[0] += rr[0]; - ll[1] += rr[1]; + Collectors.sumWithCompensation(ll, rr[0]); + Collectors.sumWithCompensation(ll, rr[1]); + ll[2] += rr[2]; }); return avg[0] > 0 - ? OptionalDouble.of(avg[1] / avg[0]) + ? OptionalDouble.of(avg[0] / avg[2]) : OptionalDouble.empty(); } --- /dev/null 2013-11-19 22:55:24.836127557 -0800 +++ new/test/java/util/stream/TestDoubleSumAverage.java 2013-11-20 00:17:20.000000000 -0800 @@ -0,0 +1,134 @@ +/* + * Copyright (c) 2013, Oracle and/or its affiliates. All rights reserved. + * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. + * + * This code is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License version 2 only, as + * published by the Free Software Foundation. + * + * This code is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * version 2 for more details (a copy is included in the LICENSE file that + * accompanied this code). + * + * You should have received a copy of the GNU General Public License version + * 2 along with this work; if not, write to the Free Software Foundation, + * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. + * + * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA + * or visit www.oracle.com if you need additional information or have any + * questions. + */ + +import java.util.*; +import java.util.stream.*; + +/* + * @test + * @bug 8006572 + */ +public class TestDoubleSumAverage { + public static void main(String... args) { + int failures = 0; + + double base = 1.0; + double increment = Math.ulp(base)/2.0; + int count = 1_000_001; + /* + * The exact sum of the test stream is 1 + 1e6*ulp(1.0) but a + * naive summation algorithm will return 1.0 since (1.0 + + * ulp(1.0)/2) will round to 1.0 again. + */ + DoubleSummaryStatistics stats = + testStream(base, increment, count).collect(DoubleSummaryStatistics::new, + DoubleSummaryStatistics::accept, + DoubleSummaryStatistics::combine); + + double expectedSum = base + (increment * (count - 1)); + double expectedAvg = expectedSum / count; + + failures += compareUlpDifference(expectedSum, stats.getSum(), 3); + failures += compareUlpDifference(expectedAvg, stats.getAverage(), 3); + + failures += compareUlpDifference(expectedSum, testStream(base, increment, count).sum(), 3); + failures += compareUlpDifference(expectedAvg, testStream(base, increment, count).average().getAsDouble(), 3); + + double collectorSum = testBoxedStream(base, increment, count). + collect(Collectors.summingDouble(d -> d)); + failures += compareUlpDifference(expectedSum, collectorSum, 3); + + double collectorAvg = testBoxedStream(base, increment, count). + collect(Collectors.averagingDouble(d -> d)); + failures += compareUlpDifference(expectedAvg, collectorAvg, 3); + + if (failures > 0) { + throw new RuntimeException("Found " + failures + " numerical failures."); + } + } + + private static DoubleStream testStream(double base, double increment, int count) { + return StreamSupport. + doubleStream(Spliterators.spliteratorUnknownSize(new TestDoubleIterator(base, + increment, + count), + Spliterator.IMMUTABLE | + Spliterator.NONNULL), + false); + } + + private static Stream testBoxedStream(double base, double increment, int count) { + TestDoubleIterator tdi = new TestDoubleIterator(base, increment, count); + Double[] tmp = new Double[count]; + int i = 0; + while(tdi.hasNext()) { + tmp[i] = tdi.next(); + i++; + } + return Stream.of(tmp); + } + + private static int compareUlpDifference(double expected, double computed, double threshold) { + double ulpDifference = Math.abs(expected - computed) / Math.ulp(expected); + + if (ulpDifference > threshold) { + System.err.printf("Numerical summation error too large, %g ulps rather than %g.%n", + ulpDifference, threshold); + return 1; + } else + return 0; + } + + static class TestDoubleIterator implements PrimitiveIterator.OfDouble { + private double base; + private double subsequent; + private long count; + private long iteration = 0; + + /** + * Construct a {@code TestDoubleStream} of {@code count} + * elements whose first value is {@code base} and whose + * following values are all {@code subsequent}. + */ + TestDoubleIterator(double base, double subsequent, long count) { + this.base = base; + this.subsequent = subsequent; + this.count = count; + } + + @Override + public boolean hasNext() { + return (iteration < count); + } + + @Override + public double nextDouble() { + if (iteration >= count) + throw new NoSuchElementException(); + else { + iteration++; + return (iteration == 1) ? base : subsequent; + } + } + } +}