| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.commons.numbers.examples.jmh.arrays; |
| |
| import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.FourD; |
| import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.ND; |
| import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.ThreeD; |
| import org.apache.commons.numbers.examples.jmh.arrays.LinearCombination.TwoD; |
| import org.apache.commons.numbers.fraction.BigFraction; |
| import org.apache.commons.rng.UniformRandomProvider; |
| import org.apache.commons.rng.simple.RandomSource; |
| import org.junit.jupiter.api.Assertions; |
| import org.junit.jupiter.api.Test; |
| import org.junit.jupiter.params.ParameterizedTest; |
| import org.junit.jupiter.params.provider.Arguments; |
| import org.junit.jupiter.params.provider.MethodSource; |
| |
| import java.math.BigDecimal; |
| import java.math.MathContext; |
| import java.math.RoundingMode; |
| import java.util.SplittableRandom; |
| import java.util.stream.Stream; |
| |
| /** |
| * Test each implementation of the LinearCombination interface. |
| */ |
| class LinearCombinationsTest { |
| /** Double.MIN_VALUE as a BigDecimal. Use string constructor to truncate precision to 4.9e-324. */ |
| private static final BigDecimal MIN = BigDecimal.valueOf(Double.MIN_VALUE); |
| |
| /** |
| * Provide instances of the LinearCombination interface as arguments. |
| * |
| * @return the stream |
| */ |
| static Stream<Arguments> provideLinearCombination() { |
| return Stream.of( |
| Arguments.of(LinearCombinations.Dekker.INSTANCE), |
| Arguments.of(LinearCombinations.Dot2s.INSTANCE), |
| Arguments.of(LinearCombinations.DotK.DOT_3), |
| Arguments.of(LinearCombinations.DotK.DOT_4), |
| Arguments.of(LinearCombinations.DotK.DOT_5), |
| Arguments.of(LinearCombinations.DotK.DOT_6), |
| Arguments.of(LinearCombinations.DotK.DOT_7), |
| Arguments.of(LinearCombinations.ExtendedPrecision.INSTANCE), |
| Arguments.of(LinearCombinations.ExtendedPrecision.DOUBLE), |
| Arguments.of(LinearCombinations.ExtendedPrecision.EXACT), |
| Arguments.of(LinearCombinations.ExtendedPrecision.EXACT2), |
| Arguments.of(LinearCombinations.Exact.INSTANCE) |
| ); |
| } |
| |
| @ParameterizedTest |
| @MethodSource("provideLinearCombination") |
| void testSingleElementArray(ND fun) { |
| final double[] a = {1.23456789}; |
| final double[] b = {98765432.1}; |
| |
| Assertions.assertEquals(a[0] * b[0], fun.value(a, b)); |
| } |
| |
| @ParameterizedTest |
| @MethodSource("provideLinearCombination") |
| void testTwoSums(ND fun) { |
| final BigFraction[] aF = new BigFraction[] { |
| BigFraction.of(-1321008684645961L, 268435456L), |
| BigFraction.of(-5774608829631843L, 268435456L), |
| BigFraction.of(-7645843051051357L, 8589934592L) |
| }; |
| final BigFraction[] bF = new BigFraction[] { |
| BigFraction.of(-5712344449280879L, 2097152L), |
| BigFraction.of(-4550117129121957L, 2097152L), |
| BigFraction.of(8846951984510141L, 131072L) |
| }; |
| final ThreeD fun3 = (ThreeD) fun; |
| |
| final int len = aF.length; |
| final double[] a = new double[len]; |
| final double[] b = new double[len]; |
| for (int i = 0; i < len; i++) { |
| a[i] = aF[i].getNumerator().doubleValue() / aF[i].getDenominator().doubleValue(); |
| b[i] = bF[i].getNumerator().doubleValue() / bF[i].getDenominator().doubleValue(); |
| } |
| |
| // Ensure "array" and "inline" implementations give the same result. |
| final double abSumInline = fun3.value(a[0], b[0], |
| a[1], b[1], |
| a[2], b[2]); |
| final double abSumArray = fun.value(a, b); |
| Assertions.assertEquals(abSumInline, abSumArray); |
| |
| // Compare with arbitrary precision computation. |
| BigFraction result = BigFraction.ZERO; |
| for (int i = 0; i < a.length; i++) { |
| result = result.add(aF[i].multiply(bF[i])); |
| } |
| final double expected = result.doubleValue(); |
| Assertions.assertEquals(expected, abSumInline, "Expecting exact result"); |
| |
| final double naive = a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; |
| Assertions.assertTrue(Math.abs(naive - abSumInline) > 1.5); |
| } |
| |
| @ParameterizedTest |
| @MethodSource("provideLinearCombination") |
| void testHuge(ND fun) { |
| final int scale = 971; |
| final double[] a = new double[] { |
| -1321008684645961.0 / 268435456.0, |
| -5774608829631843.0 / 268435456.0, |
| -7645843051051357.0 / 8589934592.0 |
| }; |
| final double[] b = new double[] { |
| -5712344449280879.0 / 2097152.0, |
| -4550117129121957.0 / 2097152.0, |
| 8846951984510141.0 / 131072.0 |
| }; |
| final ThreeD fun3 = (ThreeD) fun; |
| |
| final int len = a.length; |
| final double[] scaledA = new double[len]; |
| final double[] scaledB = new double[len]; |
| for (int i = 0; i < len; ++i) { |
| scaledA[i] = Math.scalb(a[i], -scale); |
| scaledB[i] = Math.scalb(b[i], scale); |
| } |
| final double abSumInline = fun3.value(scaledA[0], scaledB[0], |
| scaledA[1], scaledB[1], |
| scaledA[2], scaledB[2]); |
| final double abSumArray = fun.value(scaledA, scaledB); |
| |
| Assertions.assertEquals(abSumInline, abSumArray); |
| Assertions.assertEquals(-1.8551294182586248737720779899, abSumInline, "Expecting exact result"); |
| |
| final double naive = scaledA[0] * scaledB[0] + scaledA[1] * scaledB[1] + scaledA[2] * scaledB[2]; |
| Assertions.assertTrue(Math.abs(naive - abSumInline) > 1.5); |
| } |
| |
| @ParameterizedTest |
| @MethodSource("provideLinearCombination") |
| void testArrayVsInline(ND fun) { |
| // Assume the instance implements the inline functions |
| final TwoD fun2 = (TwoD) fun; |
| final ThreeD fun3 = (ThreeD) fun; |
| final FourD fun4 = (FourD) fun; |
| |
| final SplittableRandom rng = new SplittableRandom(); |
| |
| double sInline; |
| double sArray; |
| final double scale = 1e17; |
| for (int i = 0; i < 1000; ++i) { |
| final double u1 = scale * rng.nextDouble(); |
| final double u2 = scale * rng.nextDouble(); |
| final double u3 = scale * rng.nextDouble(); |
| final double u4 = scale * rng.nextDouble(); |
| final double v1 = scale * rng.nextDouble(); |
| final double v2 = scale * rng.nextDouble(); |
| final double v3 = scale * rng.nextDouble(); |
| final double v4 = scale * rng.nextDouble(); |
| |
| // One sum. |
| sInline = fun2.value(u1, v1, u2, v2); |
| sArray = fun.value(new double[] {u1, u2}, |
| new double[] {v1, v2}); |
| Assertions.assertEquals(sInline, sArray); |
| |
| // Two sums. |
| sInline = fun3.value(u1, v1, u2, v2, u3, v3); |
| sArray = fun.value(new double[] {u1, u2, u3}, |
| new double[] {v1, v2, v3}); |
| Assertions.assertEquals(sInline, sArray); |
| |
| // Three sums. |
| sInline = fun4.value(u1, v1, u2, v2, u3, v3, u4, v4); |
| sArray = fun.value(new double[] {u1, u2, u3, u4}, |
| new double[] {v1, v2, v3, v4}); |
| Assertions.assertEquals(sInline, sArray); |
| } |
| } |
| |
| @ParameterizedTest |
| @MethodSource("provideLinearCombination") |
| void testNonFinite(ND fun) { |
| final double[][] a = new double[][] { |
| {1, 2, 3, 4}, |
| {1, Double.POSITIVE_INFINITY, 3, 4}, |
| {1, 2, Double.POSITIVE_INFINITY, 4}, |
| {1, Double.POSITIVE_INFINITY, 3, Double.NEGATIVE_INFINITY}, |
| {1, 2, 3, 4}, |
| {1, 2, 3, 4}, |
| {1, 2, 3, 4}, |
| {1, 2, 3, 4}, |
| {1, Double.MAX_VALUE, 3, 4}, |
| {1, 2, Double.MAX_VALUE, 4}, |
| {1, Double.MAX_VALUE / 2, 3, -Double.MAX_VALUE / 4}, |
| }; |
| final double[][] b = new double[][] { |
| {1, -2, 3, 4}, |
| {1, -2, 3, 4}, |
| {1, -2, 3, 4}, |
| {1, -2, 3, 4}, |
| {1, Double.POSITIVE_INFINITY, 3, 4}, |
| {1, -2, Double.POSITIVE_INFINITY, 4}, |
| {1, Double.POSITIVE_INFINITY, 3, Double.NEGATIVE_INFINITY}, |
| {Double.NaN, -2, 3, 4}, |
| {1, -2, 3, 4}, |
| {1, -2, 3, 4}, |
| {1, -2, 3, 4}, |
| }; |
| |
| final TwoD fun2 = (TwoD) fun; |
| final ThreeD fun3 = (ThreeD) fun; |
| final FourD fun4 = (FourD) fun; |
| |
| Assertions.assertEquals(-3, fun2.value(a[0][0], b[0][0], |
| a[0][1], b[0][1])); |
| Assertions.assertEquals(6, fun3.value(a[0][0], b[0][0], |
| a[0][1], b[0][1], |
| a[0][2], b[0][2])); |
| Assertions.assertEquals(22, fun4.value(a[0][0], b[0][0], |
| a[0][1], b[0][1], |
| a[0][2], b[0][2], |
| a[0][3], b[0][3])); |
| Assertions.assertEquals(22, fun.value(a[0], b[0])); |
| |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun2.value(a[1][0], b[1][0], |
| a[1][1], b[1][1])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun3.value(a[1][0], b[1][0], |
| a[1][1], b[1][1], |
| a[1][2], b[1][2])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun4.value(a[1][0], b[1][0], |
| a[1][1], b[1][1], |
| a[1][2], b[1][2], |
| a[1][3], b[1][3])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun.value(a[1], b[1])); |
| |
| Assertions.assertEquals(-3, fun2.value(a[2][0], b[2][0], |
| a[2][1], b[2][1])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun3.value(a[2][0], b[2][0], |
| a[2][1], b[2][1], |
| a[2][2], b[2][2])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun4.value(a[2][0], b[2][0], |
| a[2][1], b[2][1], |
| a[2][2], b[2][2], |
| a[2][3], b[2][3])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun.value(a[2], b[2])); |
| |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun2.value(a[3][0], b[3][0], |
| a[3][1], b[3][1])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun3.value(a[3][0], b[3][0], |
| a[3][1], b[3][1], |
| a[3][2], b[3][2])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun4.value(a[3][0], b[3][0], |
| a[3][1], b[3][1], |
| a[3][2], b[3][2], |
| a[3][3], b[3][3])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun.value(a[3], b[3])); |
| |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun2.value(a[4][0], b[4][0], |
| a[4][1], b[4][1])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun3.value(a[4][0], b[4][0], |
| a[4][1], b[4][1], |
| a[4][2], b[4][2])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun4.value(a[4][0], b[4][0], |
| a[4][1], b[4][1], |
| a[4][2], b[4][2], |
| a[4][3], b[4][3])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun.value(a[4], b[4])); |
| |
| Assertions.assertEquals(-3, fun2.value(a[5][0], b[5][0], |
| a[5][1], b[5][1])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun3.value(a[5][0], b[5][0], |
| a[5][1], b[5][1], |
| a[5][2], b[5][2])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun4.value(a[5][0], b[5][0], |
| a[5][1], b[5][1], |
| a[5][2], b[5][2], |
| a[5][3], b[5][3])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun.value(a[5], b[5])); |
| |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun2.value(a[6][0], b[6][0], |
| a[6][1], b[6][1])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun3.value(a[6][0], b[6][0], |
| a[6][1], b[6][1], |
| a[6][2], b[6][2])); |
| Assertions.assertEquals(Double.NaN, fun4.value(a[6][0], b[6][0], |
| a[6][1], b[6][1], |
| a[6][2], b[6][2], |
| a[6][3], b[6][3])); |
| Assertions.assertEquals(Double.NaN, fun.value(a[6], b[6])); |
| |
| Assertions.assertEquals(Double.NaN, fun2.value(a[7][0], b[7][0], |
| a[7][1], b[7][1])); |
| Assertions.assertEquals(Double.NaN, fun3.value(a[7][0], b[7][0], |
| a[7][1], b[7][1], |
| a[7][2], b[7][2])); |
| Assertions.assertEquals(Double.NaN, fun4.value(a[7][0], b[7][0], |
| a[7][1], b[7][1], |
| a[7][2], b[7][2], |
| a[7][3], b[7][3])); |
| Assertions.assertEquals(Double.NaN, fun.value(a[7], b[7])); |
| |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun2.value(a[8][0], b[8][0], |
| a[8][1], b[8][1])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun3.value(a[8][0], b[8][0], |
| a[8][1], b[8][1], |
| a[8][2], b[8][2])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun4.value(a[8][0], b[8][0], |
| a[8][1], b[8][1], |
| a[8][2], b[8][2], |
| a[8][3], b[8][3])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun.value(a[8], b[8])); |
| |
| Assertions.assertEquals(-3, fun2.value(a[9][0], b[9][0], |
| a[9][1], b[9][1])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun3.value(a[9][0], b[9][0], |
| a[9][1], b[9][1], |
| a[9][2], b[9][2])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun4.value(a[9][0], b[9][0], |
| a[9][1], b[9][1], |
| a[9][2], b[9][2], |
| a[9][3], b[9][3])); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, fun.value(a[9], b[9])); |
| |
| Assertions.assertEquals(-Double.MAX_VALUE, fun2.value(a[10][0], b[10][0], |
| a[10][1], b[10][1])); |
| Assertions.assertEquals(-Double.MAX_VALUE, fun3.value(a[10][0], b[10][0], |
| a[10][1], b[10][1], |
| a[10][2], b[10][2])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun4.value(a[10][0], b[10][0], |
| a[10][1], b[10][1], |
| a[10][2], b[10][2], |
| a[10][3], b[10][3])); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, fun.value(a[10], b[10])); |
| } |
| |
| /** |
| * This creates a scenario where the split product will overflow but the standard |
| * precision computation will not. The result is expected to be in extended precision, |
| * i.e. the method correctly detects and handles intermediate overflow. |
| * |
| * <p>Note: This test assumes that LinearCombination computes a split number |
| * using Dekker's method. This can result in the high part of the number being |
| * greater in magnitude than the the original number due to round-off in the split. |
| */ |
| @Test |
| void testOverflow() { |
| // Create a simple dot product that is different in high precision and has |
| // values that create a high part above the original number. This can be done using |
| // a mantissa with almost all bits set to 1. |
| final double x = Math.nextDown(2.0); |
| final double y = -Math.nextDown(x); |
| final double xxMxy = x * x + x * y; |
| final double xxMxyHighPrecision = LinearCombinations.DotK.DOT_3.value(x, x, x, y); |
| Assertions.assertNotEquals(xxMxy, xxMxyHighPrecision, "High precision result should be different"); |
| |
| // Scale it close to max value. |
| // The current exponent is 0 so the combined scale must be 1023-1 as the |
| // partial product x*x and x*y have an exponent 1 higher |
| Assertions.assertEquals(0, Math.getExponent(x)); |
| Assertions.assertEquals(0, Math.getExponent(y)); |
| |
| final double a1 = Math.scalb(x, 1022 - 30); |
| final double b1 = Math.scalb(x, 30); |
| final double a2 = a1; |
| final double b2 = Math.scalb(y, 30); |
| // Verify low precision result is scaled and finite |
| final double sxxMxy = Math.scalb(xxMxy, 1022); |
| Assertions.assertEquals(sxxMxy, a1 * b1 + a2 * b2); |
| Assertions.assertTrue(Double.isFinite(sxxMxy)); |
| |
| // High precision result. |
| // First demonstrate that Dekker's split will create overflow in the high part. |
| final double m = (1 << 27) + 1; |
| double c; |
| c = a1 * m; |
| final double ha1 = c - (c - a1); |
| c = b1 * m; |
| final double hb1 = c - (c - b1); |
| c = a2 * m; |
| final double ha2 = c - (c - a2); |
| c = b2 * m; |
| final double hb2 = c - (c - b2); |
| Assertions.assertTrue(Double.isFinite(ha1)); |
| Assertions.assertTrue(Double.isFinite(hb1)); |
| Assertions.assertTrue(Double.isFinite(ha2)); |
| Assertions.assertTrue(Double.isFinite(hb2)); |
| // High part should be bigger in magnitude |
| Assertions.assertTrue(Math.abs(ha1) > Math.abs(a1)); |
| Assertions.assertTrue(Math.abs(hb1) > Math.abs(b1)); |
| Assertions.assertTrue(Math.abs(ha2) > Math.abs(a2)); |
| Assertions.assertTrue(Math.abs(hb2) > Math.abs(b2)); |
| Assertions.assertEquals(Double.POSITIVE_INFINITY, ha1 * hb1, "Expected split high part to overflow"); |
| Assertions.assertEquals(Double.NEGATIVE_INFINITY, ha2 * hb2, "Expected split high part to overflow"); |
| |
| // LinearCombination should detect and handle intermediate overflow and return the |
| // high precision result. |
| final double expected = Math.scalb(xxMxyHighPrecision, 1022); |
| Assertions.assertEquals(expected, LinearCombinations.DotK.DOT_3.value(a1, b1, a2, b2)); |
| Assertions.assertEquals(expected, LinearCombinations.DotK.DOT_3.value(a1, b1, a2, b2, 0, 0)); |
| Assertions.assertEquals(expected, LinearCombinations.DotK.DOT_3.value(a1, b1, a2, b2, 0, 0, 0, 0)); |
| Assertions.assertEquals(expected, LinearCombinations.DotK.DOT_3.value(new double[] {a1, a2}, new double[] {b1, b2})); |
| } |
| |
| /** |
| * This is an extreme case of the sum x^2 + y^2 - 1 when x^2 + y^2 are 1.0 within |
| * floating-point error but if performed using high precision subtracting 1.0 is not 0.0. |
| * This case is derived from computations on a complex cis number. |
| */ |
| @Test |
| void testCisNumber() { |
| final double theta = 5.992112452678286E-7; |
| final double x = Math.cos(theta); |
| final double y = Math.sin(theta); |
| assertValue(LinearCombinations.DotK.DOT_3.value(x, x, y, y, 1, -1), |
| new double[] {x, y, 1}, |
| new double[] {x, y, -1}); |
| } |
| |
| /** |
| * Test the sum of vectors composed of sub-vectors of [a1, a2, a3, a4] * [a1, a2, -a3, -a4] |
| * where a1^2 + a2^2 = 1 and a3^2 + a4^2 = 1 such that the sum is approximately 0 every |
| * 4 products. This is a test that is failed by various implementations that accumulate the |
| * round-off sum in single or 2-fold precision. |
| */ |
| @Test |
| void testSumZero() { |
| // Fixed seed for stability |
| final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP, 876543L); |
| final int size = 10; |
| // Create random doublets of pairs of numbers that sum to 1 or -1. |
| for (int length = 4; length <= 12; length += 4) { |
| final double[] a = new double[length]; |
| final double[] b = new double[length]; |
| for (int i = 0; i < size; i++) { |
| // Flip-flop above and below zero |
| double sign = 1; |
| for (int k = 0; k < length; k += 4) { |
| // Create 2 complex cis numbers |
| final double theta1 = rng.nextDouble() * Math.PI / 2; |
| final double theta2 = rng.nextDouble() * Math.PI / 2; |
| a[k + 0] = b[k + 0] = Math.cos(theta1); |
| a[k + 1] = b[k + 1] = Math.sin(theta1); |
| a[k + 2] = b[k + 2] = Math.cos(theta2); |
| a[k + 3] = b[k + 3] = Math.sin(theta2); |
| a[k + 0] *= sign; |
| a[k + 1] *= sign; |
| a[k + 2] *= sign; |
| a[k + 3] *= sign; |
| // Invert second pair. |
| // The sum of the pairs should be zero +/- floating point error. |
| a[k + 2] = -a[k + 2]; |
| a[k + 3] = -a[k + 3]; |
| sign = -sign; |
| } |
| assertValue(LinearCombinations.DotK.DOT_3.value(a, b), a, b); |
| } |
| } |
| } |
| |
| /** |
| * Compute the sum of the product of factors in arbitrary precision and compare it to the |
| * given value. |
| * |
| * @param value the value |
| * @param a factors |
| * @param b factors |
| */ |
| private static void assertValue(double value, double[] a, double[] b) { |
| final double expected = computeValue(a, b); |
| Assertions.assertEquals(expected, value, Math.ulp(expected), |
| () -> "Difference in Ulps = " + ulps(expected, value)); |
| } |
| |
| /** |
| * Compute the sum of the product of pairs of input data using BigDecimal. |
| * The BigDecimal is not allowed to underflow Double.MIN_VALUE. |
| * |
| * @param data the data |
| * @return the sum of products |
| */ |
| private static double computeValue(double[] a, double[] b) { |
| BigDecimal sum = new BigDecimal(a[0]).multiply(new BigDecimal(b[0])); |
| for (int i = 1; i < a.length; i++) { |
| sum = clip(sum.add(clip(new BigDecimal(a[i]).multiply(new BigDecimal(b[i]))))); |
| } |
| return sum.doubleValue(); |
| } |
| |
| /** |
| * Compute the units of least precision (ulps) between the two numbers. |
| * |
| * @param a first number |
| * @param b second number |
| * @return the ulps |
| */ |
| private static long ulps(double a, double b) { |
| long x = Double.doubleToLongBits(a); |
| long y = Double.doubleToLongBits(b); |
| if (x != y) { |
| if ((x ^ y) < 0L) { |
| // Opposite signs. Measure the combined distance to zero. |
| if (x < 0) { |
| final long tmp = x; |
| x = y; |
| y = tmp; |
| } |
| return (x - Double.doubleToLongBits(0.0)) + (y - Double.doubleToLongBits(-0.0)) + 1; |
| } |
| return Math.abs(x - y); |
| } |
| return 0; |
| } |
| |
| /** |
| * Clip the value to the minimum value that can be stored by a double. |
| * Ideally this should round BigDecimal to values occupied by sub-normal numbers. |
| * That is non-trivial so this just removes excess precision in the significand and |
| * clips it to Double.MIN_VALUE or zero if the value is very small. The ultimate use for |
| * the BigDecimal is rounded to the closest double so this method is adequate. It would |
| * take many summations of extended precision sub-normal numbers to create more |
| * than a few ULP difference to the final double value |
| * |
| * <p>In data output by the various tests the values have never been known to require |
| * clipping so this is just a safety threshold. |
| * |
| * @param a the value |
| * @return the clipped value |
| */ |
| private static BigDecimal clip(BigDecimal a) { |
| // Min value is approx 4.9e-324. Anything with fewer decimal digits to the right of the |
| // decimal point is OK. |
| if (a.scale() < 324) { |
| return a; |
| } |
| // Reduce the scale |
| final BigDecimal b = a.setScale(MIN.scale(), RoundingMode.HALF_UP); |
| // Clip to min value |
| final BigDecimal bb = b.abs(); |
| if (bb.compareTo(MIN) < 0) { |
| // Note the number may be closer to MIN than zero so do rounding |
| if (MIN.subtract(bb).compareTo(bb) < 0) { |
| // Closer to MIN |
| return a.signum() == -1 ? MIN.negate() : MIN; |
| } |
| // Closer to zero |
| return BigDecimal.ZERO; |
| } |
| // Anything above min is allowed. |
| return b; |
| } |
| |
| /** |
| * Test the clip method does what it specifies. |
| */ |
| @Test |
| void testClip() { |
| // min value is not affected |
| Assertions.assertEquals(Double.MIN_VALUE, clip(MIN).doubleValue()); |
| Assertions.assertEquals(-Double.MIN_VALUE, clip(MIN.negate()).doubleValue()); |
| // Round-up to min |
| Assertions.assertEquals(Double.MIN_VALUE, clip(MIN.divide(new BigDecimal(2))).doubleValue()); |
| Assertions.assertEquals(-Double.MIN_VALUE, clip(MIN.negate().divide(new BigDecimal(2))).doubleValue()); |
| // Round down to zero |
| Assertions.assertEquals(0, clip(MIN.divide(new BigDecimal(2.1), MathContext.DECIMAL64)).doubleValue()); |
| Assertions.assertEquals(0, clip(MIN.negate().divide(new BigDecimal(2.1), MathContext.DECIMAL64)).doubleValue()); |
| // It does not matter if BigDecimal is more precise than a sub-normal number |
| // when the output is ultimately rounded to a double. |
| Assertions.assertEquals(Double.MIN_VALUE, clip(MIN.multiply(new BigDecimal(1.1))).doubleValue()); |
| Assertions.assertEquals(Double.MIN_VALUE, clip(MIN.multiply(new BigDecimal(1.5))).doubleValue()); |
| Assertions.assertEquals(Double.MIN_VALUE * 2, clip(MIN.multiply(new BigDecimal(1.6))).doubleValue()); |
| } |
| } |
