View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
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  import java.util.function.ToDoubleFunction;
23  
24  import org.apache.commons.rng.UniformRandomProvider;
25  import org.apache.commons.rng.simple.RandomSource;
26  import org.junit.jupiter.api.Assertions;
27  import org.junit.jupiter.api.Test;
28  
29  class NormTest {
30  
31      private static final int SMALL_THRESH_EXP = -511;
32  
33      private static final int LARGE_THRESH_EXP = +496;
34  
35      private static final int RAND_VECTOR_CNT = 1_000;
36  
37      private static final int MAX_ULP_ERR = 1;
38  
39      private static final double HYPOT_COMPARE_EPS = 1e-2;
40  
41      private static final BigDecimal BD_MAX_VALUE = new BigDecimal(Double.MAX_VALUE);
42      private static final BigDecimal BD_MIN_NORMAL = new BigDecimal(Double.MIN_NORMAL);
43  
44      /** The scale, used to scale the sqrt of the sum of squares. */
45      private static final double SCALE = 0x1.0p200;
46  
47      /** The scale squared, used to scale the sum of squares. */
48      private static final BigDecimal SCALE2 = new BigDecimal(SCALE * SCALE);
49  
50      @Test
51      void testManhattan_2d() {
52          // act/assert
53          Assertions.assertEquals(0d, Norm.L1.of(0d, -0d));
54          Assertions.assertEquals(3d, Norm.L1.of(-1d, 2d));
55          Assertions.assertEquals(Double.POSITIVE_INFINITY, Norm.L1.of(Double.MAX_VALUE, Double.MAX_VALUE));
56  
57          Assertions.assertEquals(Double.NaN, Norm.L1.of(Double.NaN, 1d));
58          Assertions.assertEquals(Double.NaN, Norm.L1.of(1d, Double.NaN));
59          Assertions.assertEquals(Double.NaN, Norm.L1.of(Double.POSITIVE_INFINITY, Double.NaN));
60  
61          Assertions.assertEquals(Double.POSITIVE_INFINITY,
62                  Norm.L1.of(Double.POSITIVE_INFINITY, 0d));
63          Assertions.assertEquals(Double.POSITIVE_INFINITY,
64                  Norm.L1.of(0d, Double.POSITIVE_INFINITY));
65          Assertions.assertEquals(Double.POSITIVE_INFINITY,
66                  Norm.L1.of(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
67      }
68  
69      @Test
70      void testManhattan_3d() {
71          // act/assert
72          Assertions.assertEquals(0d, Norm.L1.of(0d, -0d, 0d));
73          Assertions.assertEquals(6d, Norm.L1.of(-1d, 2d, -3d));
74          Assertions.assertEquals(Double.POSITIVE_INFINITY, Norm.L1.of(Double.MAX_VALUE, Double.MAX_VALUE, 0d));
75  
76          Assertions.assertEquals(Double.NaN, Norm.L1.of(Double.NaN, -2d, 1d));
77          Assertions.assertEquals(Double.NaN, Norm.L1.of(-2d, Double.NaN, 1d));
78          Assertions.assertEquals(Double.NaN, Norm.L1.of(-2d, 1d, Double.NaN));
79          Assertions.assertEquals(Double.NaN, Norm.L1.of(-2d, Double.POSITIVE_INFINITY, Double.NaN));
80  
81          Assertions.assertEquals(Double.POSITIVE_INFINITY,
82                  Norm.L1.of(Double.POSITIVE_INFINITY, 2d, -4d));
83          Assertions.assertEquals(Double.POSITIVE_INFINITY,
84                  Norm.L1.of(Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, -4d));
85          Assertions.assertEquals(Double.POSITIVE_INFINITY,
86                  Norm.L1.of(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
87      }
88  
89      @Test
90      void testManhattan_array() {
91          // act/assert
92          Assertions.assertThrows(IllegalArgumentException.class, () -> Norm.L1.of(new double[0]));
93  
94          Assertions.assertEquals(0d, Norm.L1.of(new double[] {0d, -0d}));
95          Assertions.assertEquals(6d, Norm.L1.of(new double[] {-1d, 2d, -3d}));
96          Assertions.assertEquals(10d, Norm.L1.of(new double[] {-1d, 2d, -3d, 4d}));
97          Assertions.assertEquals(Double.POSITIVE_INFINITY,
98                  Norm.L1.of(new double[] {Double.MAX_VALUE, Double.MAX_VALUE}));
99  
100         Assertions.assertEquals(Double.NaN, Norm.L1.of(new double[] {-2d, Double.NaN, 1d}));
101         Assertions.assertEquals(Double.NaN, Norm.L1.of(new double[] {Double.POSITIVE_INFINITY, Double.NaN, 1d}));
102 
103         Assertions.assertEquals(Double.POSITIVE_INFINITY,
104                 Norm.L1.of(new double[] {Double.POSITIVE_INFINITY, 0d}));
105         Assertions.assertEquals(Double.POSITIVE_INFINITY,
106                 Norm.L1.of(new double[] {Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY}));
107     }
108 
109     @Test
110     void testEuclidean_2d_simple() {
111         // act/assert
112         Assertions.assertEquals(0d, Norm.L2.of(0d, 0d));
113         Assertions.assertEquals(1d, Norm.L2.of(1d, 0d));
114         Assertions.assertEquals(1d, Norm.L2.of(0d, 1d));
115         Assertions.assertEquals(5d, Norm.L2.of(-3d, 4d));
116         Assertions.assertEquals(Double.MIN_VALUE, Norm.L2.of(0d, Double.MIN_VALUE));
117         Assertions.assertEquals(Double.MAX_VALUE, Norm.L2.of(Double.MAX_VALUE, 0d));
118         Assertions.assertEquals(Double.POSITIVE_INFINITY, Norm.L2.of(Double.MAX_VALUE, Double.MAX_VALUE));
119 
120         Assertions.assertEquals(Math.sqrt(2), Norm.L2.of(1d, -1d));
121 
122         Assertions.assertEquals(Double.NaN, Norm.L2.of(Double.NaN, -2d));
123         Assertions.assertEquals(Double.NaN, Norm.L2.of(Double.NaN, Double.POSITIVE_INFINITY));
124         Assertions.assertEquals(Double.NaN, Norm.L2.of(-2d, Double.NaN));
125         Assertions.assertEquals(Double.NaN,
126                 Norm.L2.of(Double.NaN, Double.NEGATIVE_INFINITY));
127 
128         Assertions.assertEquals(Double.POSITIVE_INFINITY,
129                 Norm.L2.of(1d, Double.NEGATIVE_INFINITY));
130         Assertions.assertEquals(Double.POSITIVE_INFINITY,
131                 Norm.L2.of(Double.POSITIVE_INFINITY, -1d));
132         Assertions.assertEquals(Double.POSITIVE_INFINITY,
133                 Norm.L2.of(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
134     }
135 
136     @Test
137     void testEuclidean_2d_scaled() {
138         // arrange
139         final double[] ones = new double[] {1, 1};
140         final double[] multiplesOfTen = new double[] {1, 10};
141 
142         // act/assert
143         checkScaledEuclideanNorm(ones, 0);
144         checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP);
145         checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP + 1);
146         checkScaledEuclideanNorm(ones, -100);
147         checkScaledEuclideanNorm(ones, -101);
148         checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP);
149         checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP - 1);
150 
151         checkScaledEuclideanNorm(multiplesOfTen, 0);
152         checkScaledEuclideanNorm(multiplesOfTen, -100);
153         checkScaledEuclideanNorm(multiplesOfTen, -101);
154         checkScaledEuclideanNorm(multiplesOfTen, LARGE_THRESH_EXP - 1);
155         checkScaledEuclideanNorm(multiplesOfTen, SMALL_THRESH_EXP);
156     }
157 
158     @Test
159     void testEuclidean_2d_dominantValue() {
160         // act/assert
161         Assertions.assertEquals(Math.PI, Norm.L2.of(-Math.PI, 0x1.0p-55));
162         Assertions.assertEquals(Math.PI, Norm.L2.of(0x1.0p-55, -Math.PI));
163     }
164 
165     @Test
166     void testEuclidean_2d_random() {
167         // arrange
168         final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create(1L);
169 
170         // act/assert
171         checkEuclideanRandom(2, rng);
172     }
173 
174     @Test
175     void testEuclidean_2d_vsArray() {
176         // arrange
177         final double[][] inputs = {
178             {-4.074598908124454E-9, 9.897869969944898E-28},
179             {1.3472131556526359E-27, -9.064577177323565E9},
180             {-3.9219339341360245E149, -7.132522817112096E148},
181             {-1.4888098520466735E153, -2.9099184907796666E150},
182             {-8.659395144898396E-152, -1.123275532302136E-150},
183             {-3.660198254902351E-152, -6.656524053354807E-153}
184         };
185 
186         // act/assert
187         for (final double[] input : inputs) {
188             Assertions.assertEquals(Norm.L2.of(input), Norm.L2.of(input[0], input[1]),
189                 () -> "Expected inline method result to equal array result for input " + Arrays.toString(input));
190         }
191     }
192 
193     @Test
194     void testEuclidean_2d_vsHypot() {
195         // arrange
196         final int samples = 1000;
197         final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create(3L);
198 
199         // act/assert
200         assertEuclidean2dVersusHypot(-10, +10, samples, rng);
201         assertEuclidean2dVersusHypot(0, +20, samples, rng);
202         assertEuclidean2dVersusHypot(-20, 0, samples, rng);
203         assertEuclidean2dVersusHypot(-20, +20, samples, rng);
204         assertEuclidean2dVersusHypot(-100, +100, samples, rng);
205         assertEuclidean2dVersusHypot(LARGE_THRESH_EXP - 10, LARGE_THRESH_EXP + 10, samples, rng);
206         assertEuclidean2dVersusHypot(SMALL_THRESH_EXP - 10, SMALL_THRESH_EXP + 10, samples, rng);
207         assertEuclidean2dVersusHypot(-600, +600, samples, rng);
208     }
209 
210     /** Assert that the Norms euclidean 2D computation produces similar error behavior to Math.hypot().
211      * @param minExp minimum exponent for random inputs
212      * @param maxExp maximum exponent for random inputs
213      * @param samples sample count
214      * @param rng random number generator
215      */
216     private static void assertEuclidean2dVersusHypot(final int minExp,
217                                                      final int maxExp,
218                                                      final int samples,
219                                                      final UniformRandomProvider rng) {
220         // generate random inputs
221         final double[][] inputs = new double[samples][];
222         for (int i = 0; i < samples; ++i) {
223             inputs[i] = DoubleTestUtils.randomArray(2, minExp, maxExp, rng);
224         }
225 
226         // compute exact results
227         final double[] exactResults = new double[samples];
228         for (int i = 0; i < samples; ++i) {
229             exactResults[i] = exactEuclideanNorm(inputs[i]);
230         }
231 
232         // compute the std devs
233         final UlpErrorStats hypotStats = computeUlpErrorStats(inputs, exactResults, v -> Math.hypot(v[0], v[1]));
234         final UlpErrorStats normStats = computeUlpErrorStats(inputs, exactResults, v -> Norm.L2.of(v[0], v[1]));
235 
236         // ensure that we are within the ballpark of Math.hypot
237         Assertions.assertTrue(normStats.getMean() <= (hypotStats.getMean() + HYPOT_COMPARE_EPS),
238             () -> "Expected 2D norm result to have similar error mean to Math.hypot(): hypot error mean= " +
239                     hypotStats.getMean() + ", norm error mean= " + normStats.getMean());
240 
241         Assertions.assertTrue(normStats.getStdDev() <= (hypotStats.getStdDev() + HYPOT_COMPARE_EPS),
242             () -> "Expected 2D norm result to have similar std deviation to Math.hypot(): hypot std dev= " +
243                     hypotStats.getStdDev() + ", norm std dev= " + normStats.getStdDev());
244     }
245 
246     @Test
247     void testEuclidean_3d_simple() {
248         // act/assert
249         Assertions.assertEquals(0d, Norm.L2.of(0d, 0d, 0d));
250         Assertions.assertEquals(1d, Norm.L2.of(1d, 0d, 0d));
251         Assertions.assertEquals(1d, Norm.L2.of(0d, 1d, 0d));
252         Assertions.assertEquals(1d, Norm.L2.of(0d, 0d, 1d));
253         Assertions.assertEquals(5 * Math.sqrt(2), Norm.L2.of(-3d, -4d, 5d));
254         Assertions.assertEquals(Double.MIN_VALUE, Norm.L2.of(0d, 0d, Double.MIN_VALUE));
255         Assertions.assertEquals(Double.MAX_VALUE, Norm.L2.of(Double.MAX_VALUE, 0d, 0d));
256         Assertions.assertEquals(Double.POSITIVE_INFINITY,
257                 Norm.L2.of(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE));
258 
259         Assertions.assertEquals(Math.sqrt(3), Norm.L2.of(1d, -1d, 1d));
260 
261         Assertions.assertEquals(Double.NaN, Norm.L2.of(Double.NaN, -2d, 0d));
262         Assertions.assertEquals(Double.NaN, Norm.L2.of(-2d, Double.NaN, 0d));
263         Assertions.assertEquals(Double.NaN, Norm.L2.of(-2d, 0d, Double.NaN));
264         Assertions.assertEquals(Double.NaN,
265                 Norm.L2.of(Double.POSITIVE_INFINITY, Double.NaN, 1d));
266 
267         Assertions.assertEquals(Double.POSITIVE_INFINITY,
268                 Norm.L2.of(Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, 1d));
269         Assertions.assertEquals(Double.POSITIVE_INFINITY,
270                 Norm.L2.of(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
271     }
272 
273     @Test
274     void testEuclidean_3d_scaled() {
275         // arrange
276         final double[] ones = new double[] {1, 1, 1};
277         final double[] multiplesOfTen = new double[] {1, 10, 100};
278 
279         // act/assert
280         checkScaledEuclideanNorm(ones, 0);
281         checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP);
282         checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP + 1);
283         checkScaledEuclideanNorm(ones, -100);
284         checkScaledEuclideanNorm(ones, -101);
285         checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP);
286         checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP - 1);
287 
288         checkScaledEuclideanNorm(multiplesOfTen, 0);
289         checkScaledEuclideanNorm(multiplesOfTen, -100);
290         checkScaledEuclideanNorm(multiplesOfTen, -101);
291         checkScaledEuclideanNorm(multiplesOfTen, LARGE_THRESH_EXP - 1);
292         checkScaledEuclideanNorm(multiplesOfTen, SMALL_THRESH_EXP - 1);
293     }
294 
295     @Test
296     void testEuclidean_3d_random() {
297         // arrange
298         final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create(1L);
299 
300         // act/assert
301         checkEuclideanRandom(3, rng);
302     }
303 
304     @Test
305     void testEuclidean_3d_vsArray() {
306         // arrange
307         final double[][] inputs = {
308             {-4.074598908124454E-9, 9.897869969944898E-28, 7.849935157082846E-14},
309             {1.3472131556526359E-27, -9.064577177323565E9, 323771.526282239},
310             {-3.9219339341360245E149, -7.132522817112096E148, -3.427334456813165E147},
311             {-1.4888098520466735E153, -2.9099184907796666E150, 1.0144962310234785E152},
312             {-8.659395144898396E-152, -1.123275532302136E-150, -2.151505326692001E-152},
313             {-3.660198254902351E-152, -6.656524053354807E-153, -3.198606556986218E-154}
314         };
315 
316         // act/assert
317         for (final double[] input : inputs) {
318             Assertions.assertEquals(Norm.L2.of(input), Norm.L2.of(input[0], input[1], input[2]),
319                 () -> "Expected inline method result to equal array result for input " + Arrays.toString(input));
320         }
321     }
322 
323     @Test
324     void testEuclidean_array_simple() {
325         // act/assert
326         Assertions.assertThrows(IllegalArgumentException.class, () -> Norm.L2.of(new double[0]));
327 
328         Assertions.assertEquals(5d, Norm.L2.of(new double[] {-3d, 4d}));
329 
330         Assertions.assertEquals(Math.sqrt(2), Norm.L2.of(new double[] {1d, -1d}));
331         Assertions.assertEquals(Math.sqrt(3), Norm.L2.of(new double[] {1d, -1d, 1d}));
332         Assertions.assertEquals(2, Norm.L2.of(new double[] {1d, -1d, 1d, -1d}));
333 
334         final double[] longVec = new double[] {-0.9, 8.7, -6.5, -4.3, -2.1, 0, 1.2, 3.4, -5.6, 7.8, 9.0};
335         Assertions.assertEquals(directEuclideanNorm(longVec), Norm.L2.of(longVec));
336 
337         Assertions.assertEquals(Double.MIN_VALUE, Norm.L2.of(new double[] {0d, Double.MIN_VALUE}));
338         Assertions.assertEquals(Double.MAX_VALUE, Norm.L2.of(new double[] {Double.MAX_VALUE, 0d}));
339 
340         final double[] maxVec = new double[1000];
341         Arrays.fill(maxVec, Double.MAX_VALUE);
342         Assertions.assertEquals(Double.POSITIVE_INFINITY, Norm.L2.of(maxVec));
343 
344         final double[] largeThreshVec = new double[1000];
345         Arrays.fill(largeThreshVec, 0x1.0p496);
346         Assertions.assertEquals(Math.sqrt(largeThreshVec.length) * largeThreshVec[0], Norm.L2.of(largeThreshVec));
347 
348         Assertions.assertEquals(Double.NaN, Norm.L2.of(new double[] {-2d, Double.NaN, 1d}));
349         Assertions.assertEquals(Double.NaN,
350                 Norm.L2.of(new double[] {Double.POSITIVE_INFINITY, Double.NaN}));
351         Assertions.assertEquals(Double.POSITIVE_INFINITY,
352                 Norm.L2.of(new double[] {Double.POSITIVE_INFINITY, 1, 0}));
353         Assertions.assertEquals(Double.POSITIVE_INFINITY,
354                 Norm.L2.of(new double[] {Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY}));
355         Assertions.assertEquals(Double.POSITIVE_INFINITY,
356                 Norm.L2.of(new double[] {Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY}));
357     }
358 
359     @Test
360     void testEuclidean_array_scaled() {
361         // arrange
362         final double[] ones = new double[] {1, 1, 1, 1};
363         final double[] multiplesOfTen = new double[] {1, 10, 100, 1000};
364 
365         // act/assert
366         checkScaledEuclideanNorm(ones, 0);
367         checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP);
368         checkScaledEuclideanNorm(ones, LARGE_THRESH_EXP + 1);
369         checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP);
370         checkScaledEuclideanNorm(ones, SMALL_THRESH_EXP - 1);
371 
372         checkScaledEuclideanNorm(multiplesOfTen, 1);
373         checkScaledEuclideanNorm(multiplesOfTen, LARGE_THRESH_EXP - 1);
374         checkScaledEuclideanNorm(multiplesOfTen, SMALL_THRESH_EXP - 1);
375     }
376 
377     @Test
378     void testEuclidean_array_random() {
379         // arrange
380         final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create(1L);
381 
382         // act/assert
383         checkEuclideanRandom(2, rng);
384         checkEuclideanRandom(3, rng);
385         checkEuclideanRandom(4, rng);
386         checkEuclideanRandom(10, rng);
387         checkEuclideanRandom(100, rng);
388     }
389 
390     @Test
391     void testMaximum_2d() {
392         // act/assert
393         Assertions.assertEquals(0d, Norm.LINF.of(0d, -0d));
394         Assertions.assertEquals(2d, Norm.LINF.of(1d, -2d));
395         Assertions.assertEquals(3d, Norm.LINF.of(3d, 1d));
396         Assertions.assertEquals(Double.MAX_VALUE, Norm.LINF.of(Double.MAX_VALUE, Double.MAX_VALUE));
397 
398         Assertions.assertEquals(Double.NaN, Norm.LINF.of(Double.NaN, 0d));
399         Assertions.assertEquals(Double.NaN, Norm.LINF.of(0d, Double.NaN));
400         Assertions.assertEquals(Double.NaN, Norm.LINF.of(Double.POSITIVE_INFINITY, Double.NaN));
401 
402         Assertions.assertEquals(Double.POSITIVE_INFINITY, Norm.LINF.of(Double.POSITIVE_INFINITY, 0d));
403         Assertions.assertEquals(Double.POSITIVE_INFINITY, Norm.LINF.of(Double.NEGATIVE_INFINITY, 0d));
404         Assertions.assertEquals(Double.POSITIVE_INFINITY, Norm.LINF.of(0d, Double.NEGATIVE_INFINITY));
405         Assertions.assertEquals(Double.POSITIVE_INFINITY,
406                 Norm.LINF.of(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
407     }
408 
409     @Test
410     void testMaximum_3d() {
411         // act/assert
412         Assertions.assertEquals(0d, Norm.LINF.of(0d, -0d, 0d));
413         Assertions.assertEquals(3d, Norm.LINF.of(1d, -2d, 3d));
414         Assertions.assertEquals(4d, Norm.LINF.of(-4d, -2d, 3d));
415         Assertions.assertEquals(Double.MAX_VALUE, Norm.LINF.of(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE));
416 
417         Assertions.assertEquals(Double.NaN, Norm.LINF.of(Double.NaN, 3d, 0d));
418         Assertions.assertEquals(Double.NaN, Norm.LINF.of(3d, Double.NaN, 0d));
419         Assertions.assertEquals(Double.NaN, Norm.LINF.of(3d, 0d, Double.NaN));
420         Assertions.assertEquals(Double.NaN, Norm.LINF.of(Double.POSITIVE_INFINITY, 0d, Double.NaN));
421 
422         Assertions.assertEquals(Double.POSITIVE_INFINITY, Norm.LINF.of(Double.POSITIVE_INFINITY, 0d, 1d));
423         Assertions.assertEquals(Double.POSITIVE_INFINITY, Norm.LINF.of(0d, Double.POSITIVE_INFINITY, 1d));
424         Assertions.assertEquals(Double.POSITIVE_INFINITY, Norm.LINF.of(0d, 1d, Double.NEGATIVE_INFINITY));
425         Assertions.assertEquals(Double.POSITIVE_INFINITY,
426                 Norm.LINF.of(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY));
427     }
428 
429     @Test
430     void testMaximum_array() {
431         // act/assert
432         Assertions.assertThrows(IllegalArgumentException.class, () -> Norm.LINF.of(new double[0]));
433 
434         Assertions.assertEquals(0d, Norm.LINF.of(new double[] {0d, -0d}));
435         Assertions.assertEquals(3d, Norm.LINF.of(new double[] {-1d, 2d, -3d}));
436         Assertions.assertEquals(4d, Norm.LINF.of(new double[] {-1d, 2d, -3d, 4d}));
437         Assertions.assertEquals(Double.MAX_VALUE, Norm.LINF.of(new double[] {Double.MAX_VALUE, Double.MAX_VALUE}));
438 
439         Assertions.assertEquals(Double.NaN, Norm.LINF.of(new double[] {-2d, Double.NaN, 1d}));
440         Assertions.assertEquals(Double.NaN, Norm.LINF.of(new double[] {Double.POSITIVE_INFINITY, Double.NaN}));
441 
442         Assertions.assertEquals(Double.POSITIVE_INFINITY,
443                 Norm.LINF.of(new double[] {0d, Double.POSITIVE_INFINITY}));
444         Assertions.assertEquals(Double.POSITIVE_INFINITY,
445                 Norm.LINF.of(new double[] {Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY}));
446     }
447 
448     /** Check a number of random vectors of length {@code len} with various exponent
449      * ranges.
450      * @param len vector array length
451      * @param rng random number generator
452      * @param fn euclidean norm test function
453      */
454     private static void checkEuclideanRandom(final int len,
455                                              final UniformRandomProvider rng) {
456         checkEuclideanRandom(len, +600, +620, rng);
457         checkEuclideanRandom(len, LARGE_THRESH_EXP - 10, LARGE_THRESH_EXP + 10, rng);
458         checkEuclideanRandom(len, +400, +420, rng);
459         checkEuclideanRandom(len, +100, +120, rng);
460         checkEuclideanRandom(len, -10, +10, rng);
461         checkEuclideanRandom(len, -120, -100, rng);
462         checkEuclideanRandom(len, -420, -400, rng);
463         checkEuclideanRandom(len, SMALL_THRESH_EXP - 10, SMALL_THRESH_EXP + 10, rng);
464         checkEuclideanRandom(len, -620, -600, rng);
465 
466         checkEuclideanRandom(len, -600, +600, rng);
467     }
468 
469     /** Check a number of random vectors of length {@code len} with elements containing
470      * exponents in the range {@code [minExp, maxExp]}.
471      * @param len vector array length
472      * @param minExp min exponent
473      * @param maxExp max exponent
474      * @param rng random number generator
475      * @param fn euclidean norm test function
476      */
477     private static void checkEuclideanRandom(final int len,
478                                              final int minExp,
479                                              final int maxExp,
480                                              final UniformRandomProvider rng) {
481         for (int i = 0; i < RAND_VECTOR_CNT; ++i) {
482             // arrange
483             final double[] v = DoubleTestUtils.randomArray(len, minExp, maxExp, rng);
484 
485             final double exact = exactEuclideanNorm(v);
486             final double direct = directEuclideanNorm(v);
487 
488             // act
489             final double actual = Norm.L2.of(v);
490 
491             // assert
492             Assertions.assertTrue(Double.isFinite(actual), () ->
493                 "Computed norm was not finite; vector= " + Arrays.toString(v) + ", exact= " + exact +
494                 ", direct= " + direct + ", actual= " + actual);
495 
496             final int ulpError = Math.abs(DoubleTestUtils.computeUlpDifference(exact, actual));
497 
498             Assertions.assertTrue(ulpError <= MAX_ULP_ERR, () ->
499                 "Computed norm ulp error exceeds bounds; vector= " + Arrays.toString(v) +
500                 ", exact= " + exact + ", actual= " + actual + ", ulpError= " + ulpError);
501         }
502     }
503 
504     /** Assert that {@code directNorm(v) * 2^scaleExp = fn(v * 2^scaleExp)}.
505      * @param v unscaled vector
506      * @param scaleExp scale factor exponent
507      */
508     private static void checkScaledEuclideanNorm(final double[] v,
509                                                  final int scaleExp) {
510 
511         final double scale = Math.scalb(1d, scaleExp);
512         final double[] scaledV = new double[v.length];
513         for (int i = 0; i < v.length; ++i) {
514             scaledV[i] = v[i] * scale;
515         }
516 
517         final double norm = directEuclideanNorm(v);
518         final double scaledNorm = Norm.L2.of(scaledV);
519 
520         Assertions.assertEquals(norm * scale, scaledNorm);
521     }
522 
523     /** Direct euclidean norm computation.
524      * @param v array
525      * @return euclidean norm using direct summation.
526      */
527     private static double directEuclideanNorm(final double[] v) {
528         double n = 0;
529         for (int i = 0; i < v.length; i++) {
530             n += v[i] * v[i];
531         }
532         return Math.sqrt(n);
533     }
534 
535     /** Compute the exact double value of the vector norm using BigDecimals
536      * with a math context of {@link MathContext#DECIMAL128}.
537      * @param v array
538      * @return euclidean norm using BigDecimal with MathContext.DECIMAL128
539      */
540     private static double exactEuclideanNorm(final double[] v) {
541         final MathContext ctx = MathContext.DECIMAL128;
542 
543         BigDecimal sum = BigDecimal.ZERO;
544         for (final double d : v) {
545             sum = sum.add(new BigDecimal(d).pow(2), ctx);
546         }
547         if (sum.equals(BigDecimal.ZERO)) {
548             return 0;
549         }
550 
551         // Java 9+:
552         // sum.sqrt(ctx).doubleValue()
553 
554         // Require the sum to be in the range of a double for conversion before sqrt().
555         // We scale by a power of 2. Rescaling uses the square root of this which is also
556         // a power of 2 and can be accumulated for exact rescaling.
557         double rescale = 1.0;
558         if (sum.compareTo(BD_MIN_NORMAL) < 0) {
559             while (sum.compareTo(BD_MIN_NORMAL) < 0) {
560                 sum = sum.multiply(SCALE2);
561                 rescale /= SCALE;
562             }
563         } else if (sum.compareTo(BD_MAX_VALUE) > 0) {
564             while (sum.compareTo(BD_MAX_VALUE) > 0) {
565                 sum = sum.divide(SCALE2);
566                 rescale *= SCALE;
567             }
568         }
569 
570         return Math.sqrt(sum.doubleValue()) * rescale;
571     }
572 
573     /** Compute statistics for the ulp error of {@code fn} for the given inputs and
574      * array of exact results.
575      * @param inputs sample inputs
576      * @param exactResults array containing the exact expected results
577      * @param fn function to perform the computation
578      * @return ulp error statistics
579      */
580     private static UlpErrorStats computeUlpErrorStats(final double[][] inputs,
581                                                       final double[] exactResults,
582                                                       ToDoubleFunction<double[]> fn) {
583 
584         // compute the ulp errors for each input
585         final int[] ulpErrors = new int[inputs.length];
586         int sum = 0;
587         for (int i = 0; i < inputs.length; ++i) {
588             final double exact = exactResults[i];
589             final double actual = fn.applyAsDouble(inputs[i]);
590 
591             final int error = DoubleTestUtils.computeUlpDifference(exact, actual);
592             ulpErrors[i] = error;
593             sum += error;
594         }
595 
596         // compute the mean
597         final double mean = sum / (double) ulpErrors.length;
598 
599         // compute the std dev
600         double diffSumSq = 0d;
601         double diff;
602         for (int ulpError : ulpErrors) {
603             diff = ulpError - mean;
604             diffSumSq += diff * diff;
605         }
606 
607         final double stdDev = Math.sqrt(diffSumSq / (inputs.length - 1));
608 
609         return new UlpErrorStats(mean, stdDev);
610     }
611 
612     /** Class containing ULP error statistics. */
613     private static final class UlpErrorStats {
614 
615         private final double mean;
616 
617         private final double stdDev;
618 
619         UlpErrorStats(final double mean, final double stdDev) {
620             this.mean = mean;
621             this.stdDev = stdDev;
622         }
623 
624         public double getMean() {
625             return mean;
626         }
627 
628         public double getStdDev() {
629             return stdDev;
630         }
631     }
632 }