--- old/src/share/classes/java/util/DoubleSummaryStatistics.java 2013-11-19 14:16:50.000000000 -0800 +++ new/src/share/classes/java/util/DoubleSummaryStatistics.java 2013-11-19 14:16:50.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 / + * compensation 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-19 14:16:50.000000000 -0800 +++ new/src/share/classes/java/util/stream/Collectors.java 2013-11-19 14:16:50.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] + a[1], + 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-19 14:16:51.000000000 -0800 +++ new/src/share/classes/java/util/stream/DoublePipeline.java 2013-11-19 14:16:51.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 @@ -393,17 +400,18 @@ @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-18 16:43:35.025865568 -0800 +++ new/test/java/util/DoubleSummaryStatistics/TestDoubleSummaryStatisticsSummation.java 2013-11-19 14:16:51.000000000 -0800 @@ -0,0 +1,102 @@ +/* + * 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 TestDoubleSummaryStatisticsSummation { + public static void main(String... args) { + double base = 1.0; + double increment = Math.ulp(base)/2.0; + int count = 1_000_001; + /* + * The exact sum of this 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. + */ + TestDoubleIterator tdi = new TestDoubleIterator(base, + increment, + count); + DoubleSummaryStatistics stats = + StreamSupport.doubleStream(Spliterators.spliteratorUnknownSize(tdi, + Spliterator.IMMUTABLE | + Spliterator.NONNULL), + false).collect(DoubleSummaryStatistics::new, + DoubleSummaryStatistics::accept, + DoubleSummaryStatistics::combine); + + double expectedSum = base + (increment * (count - 1)); + compareUlpDifference(expectedSum, stats.getSum(), 3); + + double expectedAvg = expectedSum / count; + compareUlpDifference(expectedAvg, stats.getAverage(), 3); + } + + private static void compareUlpDifference(double expected, double computed, double threshold) { + double ulpDifference = Math.abs(expected - computed) / Math.ulp(expected); + + if (ulpDifference > threshold) { + System.out.printf("Numerical summation error too large, %g ulps rather than %g.%n", + ulpDifference, threshold); + throw new RuntimeException(); + } + + } + + 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; + } + } + } +}