1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.numbers.core;
18
19 import java.math.BigDecimal;
20 import java.math.MathContext;
21 import java.util.Arrays;
22
23 import org.junit.jupiter.api.Assertions;
24 import org.junit.jupiter.api.Test;
25
26 class SumTest {
27
28 @Test
29 void testSum_simple() {
30
31 Assertions.assertEquals(0d, Sum.create().getAsDouble());
32
33 assertSum(Math.PI, Math.PI);
34 assertSum(Math.PI + Math.E, Math.PI, Math.E);
35
36 assertSum(0, 0, 0, 0);
37 assertSum(6, 1, 2, 3);
38 assertSum(2, 1, -2, 3);
39
40 assertSum(Double.NaN, Double.NaN, 0, 0);
41 assertSum(Double.NaN, 0, Double.NaN, 0);
42 assertSum(Double.NaN, 0, 0, Double.NaN);
43
44 assertSum(Double.NaN, Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, 0);
45
46 assertSum(Double.POSITIVE_INFINITY, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE);
47
48 assertSum(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, 1, 1);
49 assertSum(Double.NEGATIVE_INFINITY, 1, Double.NEGATIVE_INFINITY, 1);
50 }
51
52 @Test
53 void testSumAccuracy() {
54
55 final double a = 9.999999999;
56 final double b = Math.scalb(a, -53);
57 final double c = Math.scalb(a, -53);
58 final double d = Math.scalb(a, -27);
59 final double e = Math.scalb(a, -27);
60 final double f = Math.scalb(a, -50);
61
62
63 assertSumExact(a);
64
65 assertSumExact(a, b);
66 assertSumExact(b, a);
67
68 assertSumExact(a, b, c);
69 assertSumExact(c, b, a);
70
71 assertSumExact(a, b, c, d);
72 assertSumExact(d, c, b, a);
73
74 assertSumExact(a, -b, c, -d);
75 assertSumExact(d, -c, b, -a);
76
77 assertSumExact(a, b, c, d, e, f);
78 assertSumExact(f, e, d, c, b, a);
79
80 assertSumExact(a, -b, c, -d, e, f);
81 assertSumExact(f, -e, d, -c, b, -a);
82 }
83
84 @Test
85 void testAdd_sumInstance() {
86
87 final double a = Math.PI;
88 final double b = Math.scalb(a, -53);
89 final double c = Math.scalb(a, -53);
90 final double d = Math.scalb(a, -27);
91 final double e = Math.scalb(a, -27);
92 final double f = Math.scalb(a, -50);
93
94
95 Assertions.assertEquals(exactSum(a, b, c, d), Sum.of(a, b, c, d).add(Sum.create()).getAsDouble());
96 Assertions.assertEquals(exactSum(a, a, b, c, d, e, f),
97 Sum.of(a, b)
98 .add(Sum.of(a, c))
99 .add(Sum.of(d, e, f)).getAsDouble());
100
101 final Sum s = Sum.of(a, b);
102 Assertions.assertEquals(exactSum(a, b, a, b), s.add(s).getAsDouble());
103 }
104
105 @Test
106 void testSumOfProducts_dimensionMismatch() {
107
108 final Sum sum = Sum.create();
109 Assertions.assertThrows(IllegalArgumentException.class,
110 () -> sum.addProducts(new double[1], new double[2]));
111
112 Assertions.assertThrows(IllegalArgumentException.class,
113 () -> Sum.ofProducts(new double[1], new double[2]));
114 }
115
116 @Test
117 void testSumOfProducts_singleElement() {
118 final double[] a = {1.23456789};
119 final double[] b = {98765432.1};
120
121 Assertions.assertEquals(a[0] * b[0], Sum.ofProducts(a, b).getAsDouble());
122 }
123
124 @Test
125 void testSumOfProducts() {
126
127 final BigDecimal[] aFN = new BigDecimal[] {
128 BigDecimal.valueOf(-1321008684645961L),
129 BigDecimal.valueOf(-5774608829631843L),
130 BigDecimal.valueOf(-7645843051051357L),
131 };
132 final BigDecimal[] aFD = new BigDecimal[] {
133 BigDecimal.valueOf(268435456L),
134 BigDecimal.valueOf(268435456L),
135 BigDecimal.valueOf(8589934592L)
136 };
137 final BigDecimal[] bFN = new BigDecimal[] {
138 BigDecimal.valueOf(-5712344449280879L),
139 BigDecimal.valueOf(-4550117129121957L),
140 BigDecimal.valueOf(8846951984510141L)
141 };
142 final BigDecimal[] bFD = new BigDecimal[] {
143 BigDecimal.valueOf(2097152L),
144 BigDecimal.valueOf(2097152L),
145 BigDecimal.valueOf(131072L)
146 };
147
148 final int len = aFN.length;
149 final double[] a = new double[len];
150 final double[] b = new double[len];
151 for (int i = 0; i < len; i++) {
152 a[i] = aFN[i].doubleValue() / aFD[i].doubleValue();
153 b[i] = bFN[i].doubleValue() / bFD[i].doubleValue();
154 }
155
156
157 final double sum = Sum.ofProducts(a, b).getAsDouble();
158
159
160
161 BigDecimal result = BigDecimal.ZERO;
162 for (int i = 0; i < a.length; i++) {
163 result = result.add(aFN[i].divide(aFD[i]).multiply(bFN[i].divide(bFD[i])));
164 }
165 final double expected = result.doubleValue();
166 Assertions.assertEquals(expected, sum, 1e-15);
167
168 final double naive = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
169 Assertions.assertTrue(Math.abs(naive - sum) > 1.5);
170 }
171
172 @Test
173 void testSumOfProducts_huge() {
174
175 int scale = 971;
176 final double[] a = new double[] {
177 -1321008684645961.0 / 268435456.0,
178 -5774608829631843.0 / 268435456.0,
179 -7645843051051357.0 / 8589934592.0
180 };
181 final double[] b = new double[] {
182 -5712344449280879.0 / 2097152.0,
183 -4550117129121957.0 / 2097152.0,
184 8846951984510141.0 / 131072.0
185 };
186
187 final int len = a.length;
188 final double[] scaledA = new double[len];
189 final double[] scaledB = new double[len];
190 for (int i = 0; i < len; ++i) {
191 scaledA[i] = Math.scalb(a[i], -scale);
192 scaledB[i] = Math.scalb(b[i], scale);
193 }
194
195
196 final double sum = Sum.ofProducts(scaledA, scaledB).getAsDouble();
197
198
199 Assertions.assertEquals(-1.8551294182586248737720779899, sum, 1e-15);
200
201 final double naive = scaledA[0] * scaledB[0] + scaledA[1] * scaledB[1] + scaledA[2] * scaledB[2];
202 Assertions.assertTrue(Math.abs(naive - sum) > 1.5);
203 }
204
205 @Test
206 void testSumOfProducts_nonFinite() {
207
208 final double[][] a = new double[][] {
209 {1, 2, 3, 4},
210 {1, Double.POSITIVE_INFINITY, 3, 4},
211 {1, 2, Double.POSITIVE_INFINITY, 4},
212 {1, Double.POSITIVE_INFINITY, 3, Double.NEGATIVE_INFINITY},
213 {1, 2, 3, 4},
214 {1, 2, 3, 4},
215 {1, 2, 3, 4},
216 {1, 2, 3, 4},
217 {1, Double.MAX_VALUE, 3, 4},
218 {1, 2, Double.MAX_VALUE, 4},
219 {1, Double.MAX_VALUE / 2, 3, -Double.MAX_VALUE / 4},
220 };
221 final double[][] b = new double[][] {
222 {1, -2, 3, 4},
223 {1, -2, 3, 4},
224 {1, -2, 3, 4},
225 {1, -2, 3, 4},
226 {1, Double.POSITIVE_INFINITY, 3, 4},
227 {1, -2, Double.POSITIVE_INFINITY, 4},
228 {1, Double.POSITIVE_INFINITY, 3, Double.NEGATIVE_INFINITY},
229 {Double.NaN, -2, 3, 4},
230 {1, -2, 3, 4},
231 {1, -2, 3, 4},
232 {1, -2, 3, 4},
233 };
234
235
236 assertSumOfProducts(-3,
237 a[0][0], b[0][0],
238 a[0][1], b[0][1]);
239 assertSumOfProducts(6,
240 a[0][0], b[0][0],
241 a[0][1], b[0][1],
242 a[0][2], b[0][2]);
243 assertSumOfProducts(22, a[0], b[0]);
244
245 assertSumOfProducts(Double.NEGATIVE_INFINITY,
246 a[1][0], b[1][0],
247 a[1][1], b[1][1]);
248 assertSumOfProducts(Double.NEGATIVE_INFINITY,
249 a[1][0], b[1][0],
250 a[1][1], b[1][1],
251 a[1][2], b[1][2]);
252 assertSumOfProducts(Double.NEGATIVE_INFINITY, a[1], b[1]);
253
254 assertSumOfProducts(-3,
255 a[2][0], b[2][0],
256 a[2][1], b[2][1]);
257 assertSumOfProducts(Double.POSITIVE_INFINITY,
258 a[2][0], b[2][0],
259 a[2][1], b[2][1],
260 a[2][2], b[2][2]);
261 assertSumOfProducts(Double.POSITIVE_INFINITY, a[2], b[2]);
262
263 assertSumOfProducts(Double.NEGATIVE_INFINITY,
264 a[3][0], b[3][0],
265 a[3][1], b[3][1]);
266 assertSumOfProducts(Double.NEGATIVE_INFINITY,
267 a[3][0], b[3][0],
268 a[3][1], b[3][1],
269 a[3][2], b[3][2]);
270 assertSumOfProducts(Double.NEGATIVE_INFINITY, a[3], b[3]);
271
272 assertSumOfProducts(Double.POSITIVE_INFINITY,
273 a[4][0], b[4][0],
274 a[4][1], b[4][1]);
275 assertSumOfProducts(Double.POSITIVE_INFINITY,
276 a[4][0], b[4][0],
277 a[4][1], b[4][1],
278 a[4][2], b[4][2]);
279 assertSumOfProducts(Double.POSITIVE_INFINITY, a[4], b[4]);
280
281 assertSumOfProducts(-3,
282 a[5][0], b[5][0],
283 a[5][1], b[5][1]);
284 assertSumOfProducts(Double.POSITIVE_INFINITY,
285 a[5][0], b[5][0],
286 a[5][1], b[5][1],
287 a[5][2], b[5][2]);
288 assertSumOfProducts(Double.POSITIVE_INFINITY, a[5], b[5]);
289
290 assertSumOfProducts(Double.POSITIVE_INFINITY,
291 a[6][0], b[6][0],
292 a[6][1], b[6][1]);
293 assertSumOfProducts(Double.POSITIVE_INFINITY,
294 a[6][0], b[6][0],
295 a[6][1], b[6][1],
296 a[6][2], b[6][2]);
297 assertSumOfProducts(Double.NaN, a[6], b[6]);
298
299 assertSumOfProducts(Double.NaN,
300 a[7][0], b[7][0],
301 a[7][1], b[7][1]);
302 assertSumOfProducts(Double.NaN,
303 a[7][0], b[7][0],
304 a[7][1], b[7][1],
305 a[7][2], b[7][2]);
306 assertSumOfProducts(Double.NaN, a[7], b[7]);
307
308 assertSumOfProducts(Double.NEGATIVE_INFINITY,
309 a[8][0], b[8][0],
310 a[8][1], b[8][1]);
311 assertSumOfProducts(Double.NEGATIVE_INFINITY,
312 a[8][0], b[8][0],
313 a[8][1], b[8][1],
314 a[8][2], b[8][2]);
315 assertSumOfProducts(Double.NEGATIVE_INFINITY, a[8], b[8]);
316
317 assertSumOfProducts(-3,
318 a[9][0], b[9][0],
319 a[9][1], b[9][1]);
320 assertSumOfProducts(Double.POSITIVE_INFINITY,
321 a[9][0], b[9][0],
322 a[9][1], b[9][1],
323 a[9][2], b[9][2]);
324 assertSumOfProducts(Double.POSITIVE_INFINITY, a[9], b[9]);
325
326 assertSumOfProducts(-Double.MAX_VALUE,
327 a[10][0], b[10][0],
328 a[10][1], b[10][1]);
329 assertSumOfProducts(-Double.MAX_VALUE,
330 a[10][0], b[10][0],
331 a[10][1], b[10][1],
332 a[10][2], b[10][2]);
333 assertSumOfProducts(Double.NEGATIVE_INFINITY, a[10], b[10]);
334 }
335
336
337
338
339
340
341
342
343
344
345 @Test
346 void testSumOfProducts_overflow() {
347
348
349
350 final double x = Math.nextDown(2.0);
351 final double y = -Math.nextDown(x);
352 final double xxMxy = x * x + x * y;
353 final double xxMxyHighPrecision = Sum.create().addProduct(x, x).addProduct(x, y).getAsDouble();
354 Assertions.assertNotEquals(xxMxy, xxMxyHighPrecision, "High precision result should be different");
355
356
357
358
359 Assertions.assertEquals(0, Math.getExponent(x));
360 Assertions.assertEquals(0, Math.getExponent(y));
361
362 final double a1 = Math.scalb(x, 1022 - 30);
363 final double b1 = Math.scalb(x, 30);
364 final double a2 = a1;
365 final double b2 = Math.scalb(y, 30);
366
367 final double sxxMxy = Math.scalb(xxMxy, 1022);
368 Assertions.assertEquals(sxxMxy, a1 * b1 + a2 * b2);
369 Assertions.assertTrue(Double.isFinite(sxxMxy));
370
371
372 final double m = (1 << 27) + 1;
373
374 double c;
375 c = a1 * m;
376 final double ha1 = c - (c - a1);
377 c = b1 * m;
378 final double hb1 = c - (c - b1);
379 c = a2 * m;
380 final double ha2 = c - (c - a2);
381 c = b2 * m;
382 final double hb2 = c - (c - b2);
383 Assertions.assertTrue(Double.isFinite(ha1));
384 Assertions.assertTrue(Double.isFinite(hb1));
385 Assertions.assertTrue(Double.isFinite(ha2));
386 Assertions.assertTrue(Double.isFinite(hb2));
387
388 Assertions.assertTrue(Math.abs(ha1) > Math.abs(a1));
389 Assertions.assertTrue(Math.abs(hb1) > Math.abs(b1));
390 Assertions.assertTrue(Math.abs(ha2) > Math.abs(a2));
391 Assertions.assertTrue(Math.abs(hb2) > Math.abs(b2));
392 Assertions.assertEquals(Double.POSITIVE_INFINITY, ha1 * hb1, "Expected split high part to overflow");
393 Assertions.assertEquals(Double.NEGATIVE_INFINITY, ha2 * hb2, "Expected split high part to overflow");
394
395
396
397 final double expected = Math.scalb(xxMxyHighPrecision, 1022);
398 assertSumOfProducts(expected, a1, b1, a2, b2);
399 assertSumOfProducts(expected, a1, b1, a2, b2, 0, 0);
400 assertSumOfProducts(expected, a1, b1, a2, b2, 0, 0, 0, 0);
401 }
402
403 @Test
404 void testMixedSingleTermAndProduct() {
405
406 final double a = 9.999999999;
407 final double b = Math.scalb(a, -53);
408 final double c = Math.scalb(a, -53);
409 final double d = Math.scalb(a, -27);
410
411
412 Assertions.assertEquals(exactLinearCombination(1, a, -1, b, 2, c, 4, d),
413 Sum.create()
414 .add(a)
415 .add(-b)
416 .addProduct(2, c)
417 .addProduct(d, 4).getAsDouble());
418
419 Assertions.assertEquals(exactLinearCombination(1, a, -1, b, 2, c, 4, d),
420 Sum.create()
421 .addProduct(d, 4)
422 .add(a)
423 .addProduct(2, c)
424 .add(-b).getAsDouble());
425 }
426
427 @Test
428 void testUnityValuesInProduct() {
429
430 final double a = 9.999999999;
431 final double b = Math.scalb(a, -53);
432 final double c = Math.scalb(a, -53);
433 final double d = Math.scalb(a, -27);
434
435
436 Assertions.assertEquals(exactLinearCombination(1, a, -1, b, 2, c, 4, d),
437 Sum.create()
438 .addProduct(1, a)
439 .addProduct(-1, b)
440 .addProduct(2, c)
441 .addProduct(d, 4).getAsDouble());
442
443 Assertions.assertEquals(exactLinearCombination(1, a, -1, b, 2, c, 4, d),
444 Sum.create()
445 .addProduct(a, 1)
446 .addProduct(b, -1)
447 .addProduct(2, c)
448 .addProduct(d, 4).getAsDouble());
449 }
450
451 private static void assertSumExact(final double... values) {
452 final double exact = exactSum(values);
453 assertSum(exact, values);
454 }
455
456 private static void assertSum(final double expected, final double... values) {
457
458 final int len = values.length;
459 if (len == 1) {
460 Assertions.assertEquals(expected, Sum.of(values[0]).getAsDouble());
461 } else if (len == 2) {
462 Assertions.assertEquals(expected, Sum.of(values[0], values[1]).getAsDouble());
463 } else if (len == 3) {
464 Assertions.assertEquals(expected, Sum.of(values[0], values[1], values[2]).getAsDouble());
465 } else if (len == 4) {
466 Assertions.assertEquals(expected, Sum.of(values[0], values[1], values[2], values[3]).getAsDouble());
467 }
468
469
470 final Sum addAccumulator = Sum.create();
471 for (int i = 0; i < len; ++i) {
472 addAccumulator.add(values[i]);
473 }
474 Assertions.assertEquals(expected, addAccumulator.getAsDouble());
475
476
477 final Sum acceptAccumulator = Sum.create();
478 for (int i = 0; i < len; ++i) {
479 acceptAccumulator.accept(values[i]);
480 }
481 Assertions.assertEquals(expected, acceptAccumulator.getAsDouble());
482
483
484 final Sum streamAccumulator = Sum.create();
485 Arrays.stream(values).forEach(streamAccumulator);
486 Assertions.assertEquals(expected, streamAccumulator.getAsDouble());
487
488
489 Assertions.assertEquals(expected, Sum.create().add(values).getAsDouble());
490
491
492 Assertions.assertEquals(expected, Sum.of(values).getAsDouble());
493 }
494
495 private static void assertSumOfProducts(final double expected, final double... args) {
496 final int halfLen = args.length / 2;
497
498 final double[] a = new double[halfLen];
499 final double[] b = new double[halfLen];
500 for (int i = 0; i < halfLen; ++i) {
501 a[i] = args[2 * i];
502 b[i] = args[(2 * i) + 1];
503 }
504
505 assertSumOfProducts(expected, a, b);
506 }
507
508 private static void assertSumOfProducts(final double expected, final double[] a, final double[] b) {
509 final int len = a.length;
510
511
512 final Sum accumulator = Sum.create();
513 for (int i = 0; i < len; ++i) {
514 accumulator.addProduct(a[i], b[i]);
515 }
516 Assertions.assertEquals(expected, accumulator.getAsDouble());
517
518
519 Assertions.assertEquals(expected, Sum.create().addProducts(a, b).getAsDouble());
520
521
522 Assertions.assertEquals(expected, Sum.ofProducts(a, b).getAsDouble());
523 }
524
525
526
527
528
529 private static double exactSum(final double... values) {
530 BigDecimal sum = BigDecimal.ZERO;
531 for (double value : values) {
532 sum = sum.add(new BigDecimal(value), MathContext.UNLIMITED);
533 }
534
535 return sum.doubleValue();
536 }
537
538
539
540
541
542
543 private static double exactLinearCombination(final double... values) {
544 final int halfLen = values.length / 2;
545
546 BigDecimal sum = BigDecimal.ZERO;
547 for (int i = 0; i < halfLen; ++i) {
548 final BigDecimal a = new BigDecimal(values[2 * i]);
549 final BigDecimal b = new BigDecimal(values[(2 * i) + 1]);
550
551 sum = sum.add(a.multiply(b, MathContext.UNLIMITED));
552 }
553
554 return sum.doubleValue();
555 }
556 }