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 /**
20 * <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)">Norm</a> functions.
21 *
22 * <p>The implementations provide increased numerical accuracy.
23 * Algorithms primary source is the 2005 paper
24 * <a href="https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
25 * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
26 * and Shin'ichi Oishi published in <em>SIAM J. Sci. Comput</em>.
27 */
28 public enum Norm {
29 /**
30 * <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Taxicab_norm_or_Manhattan_norm">
31 * Manhattan norm</a> (sum of the absolute values of the arguments).
32 */
33 L1(Norm::manhattan, Norm::manhattan, Norm::manhattan),
34 /** Alias for {@link #L1}. */
35 MANHATTAN(L1),
36 /** <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm">Euclidean norm</a>. */
37 L2(Norm::euclidean, Norm::euclidean, Norm::euclidean),
38 /** Alias for {@link #L2}. */
39 EUCLIDEAN(L2),
40 /**
41 * <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Maximum_norm_(special_case_of:_infinity_norm,_uniform_norm,_or_supremum_norm)">
42 * Maximum norm</a> (maximum of the absolute values of the arguments).
43 */
44 LINF(Norm::maximum, Norm::maximum, Norm::maximum),
45 /** Alias for {@link #LINF}. */
46 MAXIMUM(LINF);
47
48 /**
49 * Threshold for scaling small numbers. This value is chosen such that doubles
50 * set to this value can be squared without underflow. Values less than this must
51 * be scaled up.
52 */
53 private static final double SMALL_THRESH = 0x1.0p-511;
54 /**
55 * Threshold for scaling large numbers. This value is chosen such that 2^31 doubles
56 * set to this value can be squared and added without overflow. Values greater than
57 * this must be scaled down.
58 */
59 private static final double LARGE_THRESH = 0x1.0p+496;
60 /**
61 * Threshold for scaling up a single value by {@link #SCALE_UP} without risking
62 * overflow when the value is squared.
63 */
64 private static final double SAFE_SCALE_UP_THRESH = 0x1.0p-100;
65 /** Value used to scale down large numbers. */
66 private static final double SCALE_DOWN = 0x1.0p-600;
67 /** Value used to scale up small numbers. */
68 private static final double SCALE_UP = 0x1.0p+600;
69
70 /** Threshold for the difference between the exponents of two Euclidean 2D input values
71 * where the larger value dominates the calculation.
72 */
73 private static final int EXP_DIFF_THRESHOLD_2D = 54;
74
75 /** Function of 2 arguments. */
76 @FunctionalInterface
77 private interface Two {
78 /**
79 * @param x Argument.
80 * @param y Argument.
81 * @return the norm.
82 */
83 double of(double x, double y);
84 }
85 /** Function of 3 arguments. */
86 @FunctionalInterface
87 private interface Three {
88 /**
89 * @param x Argument.
90 * @param y Argument.
91 * @param z Argument.
92 * @return the norm.
93 */
94 double of(double x, double y, double z);
95 }
96 /** Function of array argument. */
97 @FunctionalInterface
98 private interface Array {
99 /**
100 * @param v Array of arguments.
101 * @return the norm.
102 */
103 double of(double[] v);
104 }
105
106 /** Function of 2 arguments. */
107 private final Two two;
108 /** Function of 3 arguments. */
109 private final Three three;
110 /** Function of array argument. */
111 private final Array array;
112
113 /**
114 * @param two Function of 2 arguments.
115 * @param three Function of 3 arguments.
116 * @param array Function of array argument.
117 */
118 Norm(Two two,
119 Three three,
120 Array array) {
121 this.two = two;
122 this.three = three;
123 this.array = array;
124 }
125
126 /**
127 * @param alias Alternative name.
128 */
129 Norm(Norm alias) {
130 this.two = alias.two;
131 this.three = alias.three;
132 this.array = alias.array;
133 }
134
135 /**
136 * Computes the norm.
137 *
138 * <p>Special cases:
139 * <ul>
140 * <li>If either value is {@link Double#NaN}, then the result is {@link Double#NaN}.</li>
141 * <li>If either value is infinite and the other value is not {@link Double#NaN}, then
142 * the result is {@link Double#POSITIVE_INFINITY}.</li>
143 * </ul>
144 *
145 * @param x Argument.
146 * @param y Argument.
147 * @return the norm.
148 */
149 public final double of(double x,
150 double y) {
151 return two.of(x, y);
152 }
153
154 /**
155 * Computes the norm.
156 *
157 * <p>Special cases:
158 * <ul>
159 * <li>If any value is {@link Double#NaN}, then the result is {@link Double#NaN}.</li>
160 * <li>If any value is infinite and no value is not {@link Double#NaN}, then the
161 * result is {@link Double#POSITIVE_INFINITY}.</li>
162 * </ul>
163 *
164 * @param x Argument.
165 * @param y Argument.
166 * @param z Argument.
167 * @return the norm.
168 */
169 public final double of(double x,
170 double y,
171 double z) {
172 return three.of(x, y, z);
173 }
174
175 /**
176 * Computes the norm.
177 *
178 * <p>Special cases:
179 * <ul>
180 * <li>If any value is {@link Double#NaN}, then the result is {@link Double#NaN}.</li>
181 * <li>If any value is infinite and no value is not {@link Double#NaN}, then the
182 * result is {@link Double#POSITIVE_INFINITY}.</li>
183 * </ul>
184 *
185 * @param v Argument.
186 * @return the norm.
187 * @throws IllegalArgumentException if the array is empty.
188 */
189 public final double of(double[] v) {
190 ensureNonEmpty(v);
191 return array.of(v);
192 }
193
194 /** Computes the Manhattan norm.
195 *
196 * @param x first input value
197 * @param y second input value
198 * @return \(|x| + |y|\).
199 *
200 * @see #L1
201 * @see #MANHATTAN
202 * @see #of(double,double)
203 */
204 private static double manhattan(final double x,
205 final double y) {
206 return Math.abs(x) + Math.abs(y);
207 }
208
209 /** Computes the Manhattan norm.
210 *
211 * @param x first input value
212 * @param y second input value
213 * @param z third input value
214 * @return \(|x| + |y| + |z|\)
215 *
216 * @see #L1
217 * @see #MANHATTAN
218 * @see #of(double,double,double)
219 */
220 private static double manhattan(final double x,
221 final double y,
222 final double z) {
223 return Sum.of(Math.abs(x))
224 .add(Math.abs(y))
225 .add(Math.abs(z))
226 .getAsDouble();
227 }
228
229 /** Computes the Manhattan norm.
230 *
231 * @param v input values
232 * @return \(|v_0| + ... + |v_i|\)
233 *
234 * @see #L1
235 * @see #MANHATTAN
236 * @see #of(double[])
237 */
238 private static double manhattan(final double[] v) {
239 final Sum sum = Sum.create();
240
241 for (int i = 0; i < v.length; ++i) {
242 sum.add(Math.abs(v[i]));
243 }
244
245 return sum.getAsDouble();
246 }
247
248 /** Computes the Euclidean norm.
249 * This implementation handles possible overflow or underflow.
250 *
251 * <p><strong>Comparison with Math.hypot()</strong>
252 * While not guaranteed to return the same result, this method provides
253 * similar error bounds as {@link Math#hypot(double, double)} (and may run faster on
254 * some JVM).
255 *
256 * @param x first input
257 * @param y second input
258 * @return \(\sqrt{x^2 + y^2}\).
259 *
260 * @see #L2
261 * @see #EUCLIDEAN
262 * @see #of(double,double)
263 */
264 private static double euclidean(final double x,
265 final double y) {
266 final double xabs = Math.abs(x);
267 final double yabs = Math.abs(y);
268
269 final double max;
270 final double min;
271 // the compare method considers NaN greater than other values, meaning that our
272 // check for if the max is finite later on will detect NaNs correctly
273 if (Double.compare(xabs, yabs) > 0) {
274 max = xabs;
275 min = yabs;
276 } else {
277 max = yabs;
278 min = xabs;
279 }
280
281 // if the max is not finite, then one of the inputs must not have
282 // been finite
283 if (!Double.isFinite(max)) {
284 // let the standard multiply operation determine whether to return NaN or infinite
285 return xabs * yabs;
286 } else if (Math.getExponent(max) - Math.getExponent(min) > EXP_DIFF_THRESHOLD_2D) {
287 // value is completely dominated by max; just return max
288 return max;
289 }
290
291 // compute the scale and rescale values
292 final double scale;
293 final double rescale;
294 if (max > LARGE_THRESH) {
295 scale = SCALE_DOWN;
296 rescale = SCALE_UP;
297 } else if (max < SAFE_SCALE_UP_THRESH) {
298 scale = SCALE_UP;
299 rescale = SCALE_DOWN;
300 } else {
301 scale = 1d;
302 rescale = 1d;
303 }
304
305 double sum = 0d;
306 double comp = 0d;
307
308 // add scaled x
309 double sx = xabs * scale;
310 final double px = sx * sx;
311 comp += ExtendedPrecision.squareLowUnscaled(sx, px);
312 final double sumPx = sum + px;
313 comp += ExtendedPrecision.twoSumLow(sum, px, sumPx);
314 sum = sumPx;
315
316 // add scaled y
317 double sy = yabs * scale;
318 final double py = sy * sy;
319 comp += ExtendedPrecision.squareLowUnscaled(sy, py);
320 final double sumPy = sum + py;
321 comp += ExtendedPrecision.twoSumLow(sum, py, sumPy);
322 sum = sumPy;
323
324 return Math.sqrt(sum + comp) * rescale;
325 }
326
327 /** Computes the Euclidean norm.
328 * This implementation handles possible overflow or underflow.
329 *
330 * @param x first input
331 * @param y second input
332 * @param z third input
333 * @return \(\sqrt{x^2 + y^2 + z^2}\)
334 *
335 * @see #L2
336 * @see #EUCLIDEAN
337 * @see #of(double,double,double)
338 */
339 private static double euclidean(final double x,
340 final double y,
341 final double z) {
342 final double xabs = Math.abs(x);
343 final double yabs = Math.abs(y);
344 final double zabs = Math.abs(z);
345
346 final double max = Math.max(Math.max(xabs, yabs), zabs);
347
348 // if the max is not finite, then one of the inputs must not have
349 // been finite
350 if (!Double.isFinite(max)) {
351 // let the standard multiply operation determine whether to
352 // return NaN or infinite
353 return xabs * yabs * zabs;
354 }
355
356 // compute the scale and rescale values
357 final double scale;
358 final double rescale;
359 if (max > LARGE_THRESH) {
360 scale = SCALE_DOWN;
361 rescale = SCALE_UP;
362 } else if (max < SAFE_SCALE_UP_THRESH) {
363 scale = SCALE_UP;
364 rescale = SCALE_DOWN;
365 } else {
366 scale = 1d;
367 rescale = 1d;
368 }
369
370 double sum = 0d;
371 double comp = 0d;
372
373 // add scaled x
374 double sx = xabs * scale;
375 final double px = sx * sx;
376 comp += ExtendedPrecision.squareLowUnscaled(sx, px);
377 final double sumPx = sum + px;
378 comp += ExtendedPrecision.twoSumLow(sum, px, sumPx);
379 sum = sumPx;
380
381 // add scaled y
382 double sy = yabs * scale;
383 final double py = sy * sy;
384 comp += ExtendedPrecision.squareLowUnscaled(sy, py);
385 final double sumPy = sum + py;
386 comp += ExtendedPrecision.twoSumLow(sum, py, sumPy);
387 sum = sumPy;
388
389 // add scaled z
390 final double sz = zabs * scale;
391 final double pz = sz * sz;
392 comp += ExtendedPrecision.squareLowUnscaled(sz, pz);
393 final double sumPz = sum + pz;
394 comp += ExtendedPrecision.twoSumLow(sum, pz, sumPz);
395 sum = sumPz;
396
397 return Math.sqrt(sum + comp) * rescale;
398 }
399
400 /** Computes the Euclidean norm.
401 * This implementation handles possible overflow or underflow.
402 *
403 * @param v input values
404 * @return \(\sqrt{v_0^2 + ... + v_{n-1}^2}\).
405 *
406 * @see #L2
407 * @see #EUCLIDEAN
408 * @see #of(double[])
409 */
410 private static double euclidean(final double[] v) {
411 // sum of big, normal and small numbers
412 double s1 = 0;
413 double s2 = 0;
414 double s3 = 0;
415
416 // sum compensation values
417 double c1 = 0;
418 double c2 = 0;
419 double c3 = 0;
420
421 for (int i = 0; i < v.length; ++i) {
422 final double x = Math.abs(v[i]);
423 if (!Double.isFinite(x)) {
424 // not finite; determine whether to return NaN or positive infinity
425 return euclideanNormSpecial(v, i);
426 } else if (x > LARGE_THRESH) {
427 // scale down
428 final double sx = x * SCALE_DOWN;
429
430 // compute the product and product compensation
431 final double p = sx * sx;
432 final double cp = ExtendedPrecision.squareLowUnscaled(sx, p);
433
434 // compute the running sum and sum compensation
435 final double s = s1 + p;
436 final double cs = ExtendedPrecision.twoSumLow(s1, p, s);
437
438 // update running totals
439 c1 += cp + cs;
440 s1 = s;
441 } else if (x < SMALL_THRESH) {
442 // scale up
443 final double sx = x * SCALE_UP;
444
445 // compute the product and product compensation
446 final double p = sx * sx;
447 final double cp = ExtendedPrecision.squareLowUnscaled(sx, p);
448
449 // compute the running sum and sum compensation
450 final double s = s3 + p;
451 final double cs = ExtendedPrecision.twoSumLow(s3, p, s);
452
453 // update running totals
454 c3 += cp + cs;
455 s3 = s;
456 } else {
457 // no scaling
458 // compute the product and product compensation
459 final double p = x * x;
460 final double cp = ExtendedPrecision.squareLowUnscaled(x, p);
461
462 // compute the running sum and sum compensation
463 final double s = s2 + p;
464 final double cs = ExtendedPrecision.twoSumLow(s2, p, s);
465
466 // update running totals
467 c2 += cp + cs;
468 s2 = s;
469 }
470 }
471
472 // The highest sum is the significant component. Add the next significant.
473 // Note that the "x * SCALE_DOWN * SCALE_DOWN" expressions must be executed
474 // in the order given. If the two scale factors are multiplied together first,
475 // they will underflow to zero.
476 if (s1 != 0) {
477 // add s1, s2, c1, c2
478 final double s2Adj = s2 * SCALE_DOWN * SCALE_DOWN;
479 final double sum = s1 + s2Adj;
480 final double comp = ExtendedPrecision.twoSumLow(s1, s2Adj, sum) +
481 c1 + (c2 * SCALE_DOWN * SCALE_DOWN);
482 return Math.sqrt(sum + comp) * SCALE_UP;
483 } else if (s2 != 0) {
484 // add s2, s3, c2, c3
485 final double s3Adj = s3 * SCALE_DOWN * SCALE_DOWN;
486 final double sum = s2 + s3Adj;
487 final double comp = ExtendedPrecision.twoSumLow(s2, s3Adj, sum) +
488 c2 + (c3 * SCALE_DOWN * SCALE_DOWN);
489 return Math.sqrt(sum + comp);
490 }
491 // add s3, c3
492 return Math.sqrt(s3 + c3) * SCALE_DOWN;
493 }
494
495 /** Special cases of non-finite input.
496 *
497 * @param v input vector
498 * @param start index to start examining the input vector from
499 * @return Euclidean norm special value
500 */
501 private static double euclideanNormSpecial(final double[] v,
502 final int start) {
503 for (int i = start; i < v.length; ++i) {
504 if (Double.isNaN(v[i])) {
505 return Double.NaN;
506 }
507 }
508 return Double.POSITIVE_INFINITY;
509 }
510
511 /** Computes the maximum norm.
512 *
513 * @param x first input
514 * @param y second input
515 * @return \(\max{(|x|, |y|)}\).
516 *
517 * @see #LINF
518 * @see #MAXIMUM
519 * @see #of(double,double)
520 */
521 private static double maximum(final double x,
522 final double y) {
523 return Math.max(Math.abs(x), Math.abs(y));
524 }
525
526 /** Computes the maximum norm.
527 *
528 * @param x first input
529 * @param y second input
530 * @param z third input
531 * @return \(\max{(|x|, |y|, |z|)}\).
532 *
533 * @see #LINF
534 * @see #MAXIMUM
535 * @see #of(double,double,double)
536 */
537 private static double maximum(final double x,
538 final double y,
539 final double z) {
540 return Math.max(Math.abs(x),
541 Math.max(Math.abs(y),
542 Math.abs(z)));
543 }
544
545 /** Computes the maximum norm.
546 *
547 * @param v input values
548 * @return \(\max{(|v_0|, \ldots, |v_{n-1}|)}\)
549 *
550 * @see #LINF
551 * @see #MAXIMUM
552 * @see #of(double[])
553 */
554 private static double maximum(final double[] v) {
555 double max = 0d;
556 for (int i = 0; i < v.length; ++i) {
557 max = Math.max(max, Math.abs(v[i]));
558 }
559 return max;
560 }
561
562 /**
563 * @param a Array.
564 * @throws IllegalArgumentException for zero-size array.
565 */
566 private static void ensureNonEmpty(double[] a) {
567 if (a.length == 0) {
568 throw new IllegalArgumentException("Empty array");
569 }
570 }
571 }