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 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
45 private static final double SCALE = 0x1.0p200;
46
47
48 private static final BigDecimal SCALE2 = new BigDecimal(SCALE * SCALE);
49
50 @Test
51 void testManhattan_2d() {
52
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
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
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
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
139 final double[] ones = new double[] {1, 1};
140 final double[] multiplesOfTen = new double[] {1, 10};
141
142
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
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
168 final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create(1L);
169
170
171 checkEuclideanRandom(2, rng);
172 }
173
174 @Test
175 void testEuclidean_2d_vsArray() {
176
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
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
196 final int samples = 1000;
197 final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create(3L);
198
199
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
211
212
213
214
215
216 private static void assertEuclidean2dVersusHypot(final int minExp,
217 final int maxExp,
218 final int samples,
219 final UniformRandomProvider rng) {
220
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
227 final double[] exactResults = new double[samples];
228 for (int i = 0; i < samples; ++i) {
229 exactResults[i] = exactEuclideanNorm(inputs[i]);
230 }
231
232
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
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
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
276 final double[] ones = new double[] {1, 1, 1};
277 final double[] multiplesOfTen = new double[] {1, 10, 100};
278
279
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
298 final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create(1L);
299
300
301 checkEuclideanRandom(3, rng);
302 }
303
304 @Test
305 void testEuclidean_3d_vsArray() {
306
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
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
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
362 final double[] ones = new double[] {1, 1, 1, 1};
363 final double[] multiplesOfTen = new double[] {1, 10, 100, 1000};
364
365
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
380 final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create(1L);
381
382
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
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
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
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
449
450
451
452
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
470
471
472
473
474
475
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
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
489 final double actual = Norm.L2.of(v);
490
491
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
505
506
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
524
525
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
536
537
538
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
552
553
554
555
556
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
574
575
576
577
578
579
580 private static UlpErrorStats computeUlpErrorStats(final double[][] inputs,
581 final double[] exactResults,
582 ToDoubleFunction<double[]> fn) {
583
584
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
597 final double mean = sum / (double) ulpErrors.length;
598
599
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
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 }