001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.numbers.core; 018 019/** 020 * <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)">Norm</a> functions. 021 * 022 * <p>The implementations provide increased numerical accuracy. 023 * Algorithms primary source is the 2005 paper 024 * <a href="https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547"> 025 * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump, 026 * and Shin'ichi Oishi published in <em>SIAM J. Sci. Comput</em>. 027 */ 028public enum Norm { 029 /** 030 * <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Taxicab_norm_or_Manhattan_norm"> 031 * Manhattan norm</a> (sum of the absolute values of the arguments). 032 */ 033 L1(Norm::manhattan, Norm::manhattan, Norm::manhattan), 034 /** Alias for {@link #L1}. */ 035 MANHATTAN(L1), 036 /** @see <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm">Euclidean norm</a>. */ 037 L2(Norm::euclidean, Norm::euclidean, Norm::euclidean), 038 /** Alias for {@link #L2}. */ 039 EUCLIDEAN(L2), 040 /** 041 * <a href="https://en.wikipedia.org/wiki/Norm_(mathematics)#Maximum_norm_(special_case_of:_infinity_norm,_uniform_norm,_or_supremum_norm)"> 042 * Maximum norm</a> (maximum of the absolute values of the arguments). 043 */ 044 LINF(Norm::maximum, Norm::maximum, Norm::maximum), 045 /** Alias for {@link #LINF}. */ 046 MAXIMUM(LINF); 047 048 /** 049 * Threshold for scaling small numbers. This value is chosen such that doubles 050 * set to this value can be squared without underflow. Values less than this must 051 * be scaled up. 052 */ 053 private static final double SMALL_THRESH = 0x1.0p-511; 054 /** 055 * Threshold for scaling large numbers. This value is chosen such that 2^31 doubles 056 * set to this value can be squared and added without overflow. Values greater than 057 * this must be scaled down. 058 */ 059 private static final double LARGE_THRESH = 0x1.0p+496; 060 /** 061 * Threshold for scaling up a single value by {@link #SCALE_UP} without risking 062 * overflow when the value is squared. 063 */ 064 private static final double SAFE_SCALE_UP_THRESH = 0x1.0p-100; 065 /** Value used to scale down large numbers. */ 066 private static final double SCALE_DOWN = 0x1.0p-600; 067 /** Value used to scale up small numbers. */ 068 private static final double SCALE_UP = 0x1.0p+600; 069 070 /** Threshold for the difference between the exponents of two Euclidean 2D input values 071 * where the larger value dominates the calculation. 072 */ 073 private static final int EXP_DIFF_THRESHOLD_2D = 54; 074 075 /** Function of 2 arguments. */ 076 @FunctionalInterface 077 private interface Two { 078 /** 079 * @param x Argument. 080 * @param y Argument. 081 * @return the norm. 082 */ 083 double of(double x, double y); 084 } 085 /** Function of 3 arguments. */ 086 @FunctionalInterface 087 private interface Three { 088 /** 089 * @param x Argument. 090 * @param y Argument. 091 * @param z Argument. 092 * @return the norm. 093 */ 094 double of(double x, double y, double z); 095 } 096 /** Function of array argument. */ 097 @FunctionalInterface 098 private interface Array { 099 /** 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()} (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}