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 */
017
018 package org.apache.commons.math3.util;
019
020 import java.util.List;
021 import java.util.ArrayList;
022 import java.util.Comparator;
023 import java.util.Collections;
024
025 import org.apache.commons.math3.exception.DimensionMismatchException;
026 import org.apache.commons.math3.exception.MathInternalError;
027 import org.apache.commons.math3.exception.NonMonotonicSequenceException;
028 import org.apache.commons.math3.exception.NotPositiveException;
029 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
030 import org.apache.commons.math3.exception.NullArgumentException;
031 import org.apache.commons.math3.exception.MathIllegalArgumentException;
032 import org.apache.commons.math3.exception.util.LocalizedFormats;
033 import org.apache.commons.math3.exception.MathArithmeticException;
034
035 /**
036 * Arrays utilities.
037 *
038 * @since 3.0
039 * @version $Id: MathArrays.java 1422313 2012-12-15 18:53:41Z psteitz $
040 */
041 public class MathArrays {
042 /** Factor used for splitting double numbers: n = 2^27 + 1 (i.e. {@value}). */
043 private static final int SPLIT_FACTOR = 0x8000001;
044
045 /**
046 * Private constructor.
047 */
048 private MathArrays() {}
049
050 /**
051 * Real-valued function that operate on an array or a part of it.
052 * @since 3.1
053 */
054 public interface Function {
055 /**
056 * Operates on an entire array.
057 *
058 * @param array Array to operate on.
059 * @return the result of the operation.
060 */
061 double evaluate(double[] array);
062 /**
063 * @param array Array to operate on.
064 * @param startIndex Index of the first element to take into account.
065 * @param numElements Number of elements to take into account.
066 * @return the result of the operation.
067 */
068 double evaluate(double[] array,
069 int startIndex,
070 int numElements);
071 }
072
073 /**
074 * Creates an array whose contents will be the element-by-element
075 * addition of the arguments.
076 *
077 * @param a First term of the addition.
078 * @param b Second term of the addition.
079 * @return a new array {@code r} where {@code r[i] = a[i] + b[i]}.
080 * @throws DimensionMismatchException if the array lengths differ.
081 * @since 3.1
082 */
083 public static double[] ebeAdd(double[] a,
084 double[] b) {
085 if (a.length != b.length) {
086 throw new DimensionMismatchException(a.length, b.length);
087 }
088
089 final double[] result = a.clone();
090 for (int i = 0; i < a.length; i++) {
091 result[i] += b[i];
092 }
093 return result;
094 }
095 /**
096 * Creates an array whose contents will be the element-by-element
097 * subtraction of the second argument from the first.
098 *
099 * @param a First term.
100 * @param b Element to be subtracted.
101 * @return a new array {@code r} where {@code r[i] = a[i] - b[i]}.
102 * @throws DimensionMismatchException if the array lengths differ.
103 * @since 3.1
104 */
105 public static double[] ebeSubtract(double[] a,
106 double[] b) {
107 if (a.length != b.length) {
108 throw new DimensionMismatchException(a.length, b.length);
109 }
110
111 final double[] result = a.clone();
112 for (int i = 0; i < a.length; i++) {
113 result[i] -= b[i];
114 }
115 return result;
116 }
117 /**
118 * Creates an array whose contents will be the element-by-element
119 * multiplication of the arguments.
120 *
121 * @param a First factor of the multiplication.
122 * @param b Second factor of the multiplication.
123 * @return a new array {@code r} where {@code r[i] = a[i] * b[i]}.
124 * @throws DimensionMismatchException if the array lengths differ.
125 * @since 3.1
126 */
127 public static double[] ebeMultiply(double[] a,
128 double[] b) {
129 if (a.length != b.length) {
130 throw new DimensionMismatchException(a.length, b.length);
131 }
132
133 final double[] result = a.clone();
134 for (int i = 0; i < a.length; i++) {
135 result[i] *= b[i];
136 }
137 return result;
138 }
139 /**
140 * Creates an array whose contents will be the element-by-element
141 * division of the first argument by the second.
142 *
143 * @param a Numerator of the division.
144 * @param b Denominator of the division.
145 * @return a new array {@code r} where {@code r[i] = a[i] / b[i]}.
146 * @throws DimensionMismatchException if the array lengths differ.
147 * @since 3.1
148 */
149 public static double[] ebeDivide(double[] a,
150 double[] b) {
151 if (a.length != b.length) {
152 throw new DimensionMismatchException(a.length, b.length);
153 }
154
155 final double[] result = a.clone();
156 for (int i = 0; i < a.length; i++) {
157 result[i] /= b[i];
158 }
159 return result;
160 }
161
162 /**
163 * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
164 *
165 * @param p1 the first point
166 * @param p2 the second point
167 * @return the L<sub>1</sub> distance between the two points
168 */
169 public static double distance1(double[] p1, double[] p2) {
170 double sum = 0;
171 for (int i = 0; i < p1.length; i++) {
172 sum += FastMath.abs(p1[i] - p2[i]);
173 }
174 return sum;
175 }
176
177 /**
178 * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
179 *
180 * @param p1 the first point
181 * @param p2 the second point
182 * @return the L<sub>1</sub> distance between the two points
183 */
184 public static int distance1(int[] p1, int[] p2) {
185 int sum = 0;
186 for (int i = 0; i < p1.length; i++) {
187 sum += FastMath.abs(p1[i] - p2[i]);
188 }
189 return sum;
190 }
191
192 /**
193 * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
194 *
195 * @param p1 the first point
196 * @param p2 the second point
197 * @return the L<sub>2</sub> distance between the two points
198 */
199 public static double distance(double[] p1, double[] p2) {
200 double sum = 0;
201 for (int i = 0; i < p1.length; i++) {
202 final double dp = p1[i] - p2[i];
203 sum += dp * dp;
204 }
205 return FastMath.sqrt(sum);
206 }
207
208 /**
209 * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
210 *
211 * @param p1 the first point
212 * @param p2 the second point
213 * @return the L<sub>2</sub> distance between the two points
214 */
215 public static double distance(int[] p1, int[] p2) {
216 double sum = 0;
217 for (int i = 0; i < p1.length; i++) {
218 final double dp = p1[i] - p2[i];
219 sum += dp * dp;
220 }
221 return FastMath.sqrt(sum);
222 }
223
224 /**
225 * Calculates the L<sub>∞</sub> (max of abs) distance between two points.
226 *
227 * @param p1 the first point
228 * @param p2 the second point
229 * @return the L<sub>∞</sub> distance between the two points
230 */
231 public static double distanceInf(double[] p1, double[] p2) {
232 double max = 0;
233 for (int i = 0; i < p1.length; i++) {
234 max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
235 }
236 return max;
237 }
238
239 /**
240 * Calculates the L<sub>∞</sub> (max of abs) distance between two points.
241 *
242 * @param p1 the first point
243 * @param p2 the second point
244 * @return the L<sub>∞</sub> distance between the two points
245 */
246 public static int distanceInf(int[] p1, int[] p2) {
247 int max = 0;
248 for (int i = 0; i < p1.length; i++) {
249 max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
250 }
251 return max;
252 }
253
254 /**
255 * Specification of ordering direction.
256 */
257 public static enum OrderDirection {
258 /** Constant for increasing direction. */
259 INCREASING,
260 /** Constant for decreasing direction. */
261 DECREASING
262 }
263
264 /**
265 * Check that an array is monotonically increasing or decreasing.
266 *
267 * @param <T> the type of the elements in the specified array
268 * @param val Values.
269 * @param dir Ordering direction.
270 * @param strict Whether the order should be strict.
271 * @return {@code true} if sorted, {@code false} otherwise.
272 */
273 public static <T extends Comparable<? super T>> boolean isMonotonic(T[] val,
274 OrderDirection dir,
275 boolean strict) {
276 T previous = val[0];
277 final int max = val.length;
278 for (int i = 1; i < max; i++) {
279 final int comp;
280 switch (dir) {
281 case INCREASING:
282 comp = previous.compareTo(val[i]);
283 if (strict) {
284 if (comp >= 0) {
285 return false;
286 }
287 } else {
288 if (comp > 0) {
289 return false;
290 }
291 }
292 break;
293 case DECREASING:
294 comp = val[i].compareTo(previous);
295 if (strict) {
296 if (comp >= 0) {
297 return false;
298 }
299 } else {
300 if (comp > 0) {
301 return false;
302 }
303 }
304 break;
305 default:
306 // Should never happen.
307 throw new MathInternalError();
308 }
309
310 previous = val[i];
311 }
312 return true;
313 }
314
315 /**
316 * Check that an array is monotonically increasing or decreasing.
317 *
318 * @param val Values.
319 * @param dir Ordering direction.
320 * @param strict Whether the order should be strict.
321 * @return {@code true} if sorted, {@code false} otherwise.
322 */
323 public static boolean isMonotonic(double[] val,
324 OrderDirection dir,
325 boolean strict) {
326 return checkOrder(val, dir, strict, false);
327 }
328
329 /**
330 * Check that the given array is sorted.
331 *
332 * @param val Values.
333 * @param dir Ordering direction.
334 * @param strict Whether the order should be strict.
335 * @param abort Whether to throw an exception if the check fails.
336 * @return {@code true} if the array is sorted.
337 * @throws NonMonotonicSequenceException if the array is not sorted
338 * and {@code abort} is {@code true}.
339 */
340 public static boolean checkOrder(double[] val, OrderDirection dir,
341 boolean strict, boolean abort)
342 throws NonMonotonicSequenceException {
343 double previous = val[0];
344 final int max = val.length;
345
346 int index;
347 ITEM:
348 for (index = 1; index < max; index++) {
349 switch (dir) {
350 case INCREASING:
351 if (strict) {
352 if (val[index] <= previous) {
353 break ITEM;
354 }
355 } else {
356 if (val[index] < previous) {
357 break ITEM;
358 }
359 }
360 break;
361 case DECREASING:
362 if (strict) {
363 if (val[index] >= previous) {
364 break ITEM;
365 }
366 } else {
367 if (val[index] > previous) {
368 break ITEM;
369 }
370 }
371 break;
372 default:
373 // Should never happen.
374 throw new MathInternalError();
375 }
376
377 previous = val[index];
378 }
379
380 if (index == max) {
381 // Loop completed.
382 return true;
383 }
384
385 // Loop early exit means wrong ordering.
386 if (abort) {
387 throw new NonMonotonicSequenceException(val[index], previous, index, dir, strict);
388 } else {
389 return false;
390 }
391 }
392
393 /**
394 * Check that the given array is sorted.
395 *
396 * @param val Values.
397 * @param dir Ordering direction.
398 * @param strict Whether the order should be strict.
399 * @throws NonMonotonicSequenceException if the array is not sorted.
400 * @since 2.2
401 */
402 public static void checkOrder(double[] val, OrderDirection dir,
403 boolean strict) throws NonMonotonicSequenceException {
404 checkOrder(val, dir, strict, true);
405 }
406
407 /**
408 * Check that the given array is sorted in strictly increasing order.
409 *
410 * @param val Values.
411 * @throws NonMonotonicSequenceException if the array is not sorted.
412 * @since 2.2
413 */
414 public static void checkOrder(double[] val) throws NonMonotonicSequenceException {
415 checkOrder(val, OrderDirection.INCREASING, true);
416 }
417
418 /**
419 * Throws DimensionMismatchException if the input array is not rectangular.
420 *
421 * @param in array to be tested
422 * @throws NullArgumentException if input array is null
423 * @throws DimensionMismatchException if input array is not rectangular
424 * @since 3.1
425 */
426 public static void checkRectangular(final long[][] in)
427 throws NullArgumentException, DimensionMismatchException {
428 MathUtils.checkNotNull(in);
429 for (int i = 1; i < in.length; i++) {
430 if (in[i].length != in[0].length) {
431 throw new DimensionMismatchException(
432 LocalizedFormats.DIFFERENT_ROWS_LENGTHS,
433 in[i].length, in[0].length);
434 }
435 }
436 }
437
438 /**
439 * Check that all entries of the input array are strictly positive.
440 *
441 * @param in Array to be tested
442 * @throws NotStrictlyPositiveException if any entries of the array are not
443 * strictly positive.
444 * @since 3.1
445 */
446 public static void checkPositive(final double[] in)
447 throws NotStrictlyPositiveException {
448 for (int i = 0; i < in.length; i++) {
449 if (in[i] <= 0) {
450 throw new NotStrictlyPositiveException(in[i]);
451 }
452 }
453 }
454
455 /**
456 * Check that all entries of the input array are >= 0.
457 *
458 * @param in Array to be tested
459 * @throws NotPositiveException if any array entries are less than 0.
460 * @since 3.1
461 */
462 public static void checkNonNegative(final long[] in)
463 throws NotPositiveException {
464 for (int i = 0; i < in.length; i++) {
465 if (in[i] < 0) {
466 throw new NotPositiveException(in[i]);
467 }
468 }
469 }
470
471 /**
472 * Check all entries of the input array are >= 0.
473 *
474 * @param in Array to be tested
475 * @throws NotPositiveException if any array entries are less than 0.
476 * @since 3.1
477 */
478 public static void checkNonNegative(final long[][] in)
479 throws NotPositiveException {
480 for (int i = 0; i < in.length; i ++) {
481 for (int j = 0; j < in[i].length; j++) {
482 if (in[i][j] < 0) {
483 throw new NotPositiveException(in[i][j]);
484 }
485 }
486 }
487 }
488
489 /**
490 * Returns the Cartesian norm (2-norm), handling both overflow and underflow.
491 * Translation of the minpack enorm subroutine.
492 *
493 * The redistribution policy for MINPACK is available
494 * <a href="http://www.netlib.org/minpack/disclaimer">here</a>, for
495 * convenience, it is reproduced below.</p>
496 *
497 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
498 * <tr><td>
499 * Minpack Copyright Notice (1999) University of Chicago.
500 * All rights reserved
501 * </td></tr>
502 * <tr><td>
503 * Redistribution and use in source and binary forms, with or without
504 * modification, are permitted provided that the following conditions
505 * are met:
506 * <ol>
507 * <li>Redistributions of source code must retain the above copyright
508 * notice, this list of conditions and the following disclaimer.</li>
509 * <li>Redistributions in binary form must reproduce the above
510 * copyright notice, this list of conditions and the following
511 * disclaimer in the documentation and/or other materials provided
512 * with the distribution.</li>
513 * <li>The end-user documentation included with the redistribution, if any,
514 * must include the following acknowledgment:
515 * {@code This product includes software developed by the University of
516 * Chicago, as Operator of Argonne National Laboratory.}
517 * Alternately, this acknowledgment may appear in the software itself,
518 * if and wherever such third-party acknowledgments normally appear.</li>
519 * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
520 * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
521 * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
522 * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
523 * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
524 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
525 * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
526 * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
527 * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
528 * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
529 * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
530 * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
531 * BE CORRECTED.</strong></li>
532 * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
533 * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
534 * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
535 * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
536 * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
537 * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
538 * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
539 * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
540 * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
541 * POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
542 * <ol></td></tr>
543 * </table>
544 *
545 * @param v Vector of doubles.
546 * @return the 2-norm of the vector.
547 * @since 2.2
548 */
549 public static double safeNorm(double[] v) {
550 double rdwarf = 3.834e-20;
551 double rgiant = 1.304e+19;
552 double s1 = 0;
553 double s2 = 0;
554 double s3 = 0;
555 double x1max = 0;
556 double x3max = 0;
557 double floatn = v.length;
558 double agiant = rgiant / floatn;
559 for (int i = 0; i < v.length; i++) {
560 double xabs = Math.abs(v[i]);
561 if (xabs < rdwarf || xabs > agiant) {
562 if (xabs > rdwarf) {
563 if (xabs > x1max) {
564 double r = x1max / xabs;
565 s1= 1 + s1 * r * r;
566 x1max = xabs;
567 } else {
568 double r = xabs / x1max;
569 s1 += r * r;
570 }
571 } else {
572 if (xabs > x3max) {
573 double r = x3max / xabs;
574 s3= 1 + s3 * r * r;
575 x3max = xabs;
576 } else {
577 if (xabs != 0) {
578 double r = xabs / x3max;
579 s3 += r * r;
580 }
581 }
582 }
583 } else {
584 s2 += xabs * xabs;
585 }
586 }
587 double norm;
588 if (s1 != 0) {
589 norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max);
590 } else {
591 if (s2 == 0) {
592 norm = x3max * Math.sqrt(s3);
593 } else {
594 if (s2 >= x3max) {
595 norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
596 } else {
597 norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
598 }
599 }
600 }
601 return norm;
602 }
603
604 /**
605 * Sort an array in ascending order in place and perform the same reordering
606 * of entries on other arrays. For example, if
607 * {@code x = [3, 1, 2], y = [1, 2, 3]} and {@code z = [0, 5, 7]}, then
608 * {@code sortInPlace(x, y, z)} will update {@code x} to {@code [1, 2, 3]},
609 * {@code y} to {@code [2, 3, 1]} and {@code z} to {@code [5, 7, 0]}.
610 *
611 * @param x Array to be sorted and used as a pattern for permutation
612 * of the other arrays.
613 * @param yList Set of arrays whose permutations of entries will follow
614 * those performed on {@code x}.
615 * @throws DimensionMismatchException if any {@code y} is not the same
616 * size as {@code x}.
617 * @throws NullArgumentException if {@code x} or any {@code y} is null.
618 * @since 3.0
619 */
620 public static void sortInPlace(double[] x, double[] ... yList)
621 throws DimensionMismatchException, NullArgumentException {
622 sortInPlace(x, OrderDirection.INCREASING, yList);
623 }
624
625 /**
626 * Sort an array in place and perform the same reordering of entries on
627 * other arrays. This method works the same as the other
628 * {@link #sortInPlace(double[], double[][]) sortInPlace} method, but
629 * allows the order of the sort to be provided in the {@code dir}
630 * parameter.
631 *
632 * @param x Array to be sorted and used as a pattern for permutation
633 * of the other arrays.
634 * @param dir Order direction.
635 * @param yList Set of arrays whose permutations of entries will follow
636 * those performed on {@code x}.
637 * @throws DimensionMismatchException if any {@code y} is not the same
638 * size as {@code x}.
639 * @throws NullArgumentException if {@code x} or any {@code y} is null
640 * @since 3.0
641 */
642 public static void sortInPlace(double[] x,
643 final OrderDirection dir,
644 double[] ... yList)
645 throws NullArgumentException, DimensionMismatchException {
646 if (x == null) {
647 throw new NullArgumentException();
648 }
649
650 final int len = x.length;
651 final List<Pair<Double, double[]>> list
652 = new ArrayList<Pair<Double, double[]>>(len);
653
654 final int yListLen = yList.length;
655 for (int i = 0; i < len; i++) {
656 final double[] yValues = new double[yListLen];
657 for (int j = 0; j < yListLen; j++) {
658 double[] y = yList[j];
659 if (y == null) {
660 throw new NullArgumentException();
661 }
662 if (y.length != len) {
663 throw new DimensionMismatchException(y.length, len);
664 }
665 yValues[j] = y[i];
666 }
667 list.add(new Pair<Double, double[]>(x[i], yValues));
668 }
669
670 final Comparator<Pair<Double, double[]>> comp
671 = new Comparator<Pair<Double, double[]>>() {
672 public int compare(Pair<Double, double[]> o1,
673 Pair<Double, double[]> o2) {
674 int val;
675 switch (dir) {
676 case INCREASING:
677 val = o1.getKey().compareTo(o2.getKey());
678 break;
679 case DECREASING:
680 val = o2.getKey().compareTo(o1.getKey());
681 break;
682 default:
683 // Should never happen.
684 throw new MathInternalError();
685 }
686 return val;
687 }
688 };
689
690 Collections.sort(list, comp);
691
692 for (int i = 0; i < len; i++) {
693 final Pair<Double, double[]> e = list.get(i);
694 x[i] = e.getKey();
695 final double[] yValues = e.getValue();
696 for (int j = 0; j < yListLen; j++) {
697 yList[j][i] = yValues[j];
698 }
699 }
700 }
701
702 /**
703 * Creates a copy of the {@code source} array.
704 *
705 * @param source Array to be copied.
706 * @return the copied array.
707 */
708 public static int[] copyOf(int[] source) {
709 return copyOf(source, source.length);
710 }
711
712 /**
713 * Creates a copy of the {@code source} array.
714 *
715 * @param source Array to be copied.
716 * @return the copied array.
717 */
718 public static double[] copyOf(double[] source) {
719 return copyOf(source, source.length);
720 }
721
722 /**
723 * Creates a copy of the {@code source} array.
724 *
725 * @param source Array to be copied.
726 * @param len Number of entries to copy. If smaller then the source
727 * length, the copy will be truncated, if larger it will padded with
728 * zeroes.
729 * @return the copied array.
730 */
731 public static int[] copyOf(int[] source, int len) {
732 final int[] output = new int[len];
733 System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length));
734 return output;
735 }
736
737 /**
738 * Creates a copy of the {@code source} array.
739 *
740 * @param source Array to be copied.
741 * @param len Number of entries to copy. If smaller then the source
742 * length, the copy will be truncated, if larger it will padded with
743 * zeroes.
744 * @return the copied array.
745 */
746 public static double[] copyOf(double[] source, int len) {
747 final double[] output = new double[len];
748 System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length));
749 return output;
750 }
751
752 /**
753 * Compute a linear combination accurately.
754 * This method computes the sum of the products
755 * <code>a<sub>i</sub> b<sub>i</sub></code> to high accuracy.
756 * It does so by using specific multiplication and addition algorithms to
757 * preserve accuracy and reduce cancellation effects.
758 * <br/>
759 * It is based on the 2005 paper
760 * <a href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
761 * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
762 * and Shin'ichi Oishi published in SIAM J. Sci. Comput.
763 *
764 * @param a Factors.
765 * @param b Factors.
766 * @return <code>Σ<sub>i</sub> a<sub>i</sub> b<sub>i</sub></code>.
767 * @throws DimensionMismatchException if arrays dimensions don't match
768 */
769 public static double linearCombination(final double[] a, final double[] b)
770 throws DimensionMismatchException {
771 final int len = a.length;
772 if (len != b.length) {
773 throw new DimensionMismatchException(len, b.length);
774 }
775
776 final double[] prodHigh = new double[len];
777 double prodLowSum = 0;
778
779 for (int i = 0; i < len; i++) {
780 final double ai = a[i];
781 final double ca = SPLIT_FACTOR * ai;
782 final double aHigh = ca - (ca - ai);
783 final double aLow = ai - aHigh;
784
785 final double bi = b[i];
786 final double cb = SPLIT_FACTOR * bi;
787 final double bHigh = cb - (cb - bi);
788 final double bLow = bi - bHigh;
789 prodHigh[i] = ai * bi;
790 final double prodLow = aLow * bLow - (((prodHigh[i] -
791 aHigh * bHigh) -
792 aLow * bHigh) -
793 aHigh * bLow);
794 prodLowSum += prodLow;
795 }
796
797
798 final double prodHighCur = prodHigh[0];
799 double prodHighNext = prodHigh[1];
800 double sHighPrev = prodHighCur + prodHighNext;
801 double sPrime = sHighPrev - prodHighNext;
802 double sLowSum = (prodHighNext - (sHighPrev - sPrime)) + (prodHighCur - sPrime);
803
804 final int lenMinusOne = len - 1;
805 for (int i = 1; i < lenMinusOne; i++) {
806 prodHighNext = prodHigh[i + 1];
807 final double sHighCur = sHighPrev + prodHighNext;
808 sPrime = sHighCur - prodHighNext;
809 sLowSum += (prodHighNext - (sHighCur - sPrime)) + (sHighPrev - sPrime);
810 sHighPrev = sHighCur;
811 }
812
813 double result = sHighPrev + (prodLowSum + sLowSum);
814
815 if (Double.isNaN(result)) {
816 // either we have split infinite numbers or some coefficients were NaNs,
817 // just rely on the naive implementation and let IEEE754 handle this
818 result = 0;
819 for (int i = 0; i < len; ++i) {
820 result += a[i] * b[i];
821 }
822 }
823
824 return result;
825 }
826
827 /**
828 * Compute a linear combination accurately.
829 * <p>
830 * This method computes a<sub>1</sub>×b<sub>1</sub> +
831 * a<sub>2</sub>×b<sub>2</sub> to high accuracy. It does
832 * so by using specific multiplication and addition algorithms to
833 * preserve accuracy and reduce cancellation effects. It is based
834 * on the 2005 paper <a
835 * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
836 * Accurate Sum and Dot Product</a> by Takeshi Ogita,
837 * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
838 * </p>
839 * @param a1 first factor of the first term
840 * @param b1 second factor of the first term
841 * @param a2 first factor of the second term
842 * @param b2 second factor of the second term
843 * @return a<sub>1</sub>×b<sub>1</sub> +
844 * a<sub>2</sub>×b<sub>2</sub>
845 * @see #linearCombination(double, double, double, double, double, double)
846 * @see #linearCombination(double, double, double, double, double, double, double, double)
847 */
848 public static double linearCombination(final double a1, final double b1,
849 final double a2, final double b2) {
850
851 // the code below is split in many additions/subtractions that may
852 // appear redundant. However, they should NOT be simplified, as they
853 // use IEEE754 floating point arithmetic rounding properties.
854 // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
855 // The variable naming conventions are that xyzHigh contains the most significant
856 // bits of xyz and xyzLow contains its least significant bits. So theoretically
857 // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
858 // be represented in only one double precision number so we preserve two numbers
859 // to hold it as long as we can, combining the high and low order bits together
860 // only at the end, after cancellation may have occurred on high order bits
861
862 // split a1 and b1 as two 26 bits numbers
863 final double ca1 = SPLIT_FACTOR * a1;
864 final double a1High = ca1 - (ca1 - a1);
865 final double a1Low = a1 - a1High;
866 final double cb1 = SPLIT_FACTOR * b1;
867 final double b1High = cb1 - (cb1 - b1);
868 final double b1Low = b1 - b1High;
869
870 // accurate multiplication a1 * b1
871 final double prod1High = a1 * b1;
872 final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
873
874 // split a2 and b2 as two 26 bits numbers
875 final double ca2 = SPLIT_FACTOR * a2;
876 final double a2High = ca2 - (ca2 - a2);
877 final double a2Low = a2 - a2High;
878 final double cb2 = SPLIT_FACTOR * b2;
879 final double b2High = cb2 - (cb2 - b2);
880 final double b2Low = b2 - b2High;
881
882 // accurate multiplication a2 * b2
883 final double prod2High = a2 * b2;
884 final double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
885
886 // accurate addition a1 * b1 + a2 * b2
887 final double s12High = prod1High + prod2High;
888 final double s12Prime = s12High - prod2High;
889 final double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
890
891 // final rounding, s12 may have suffered many cancellations, we try
892 // to recover some bits from the extra words we have saved up to now
893 double result = s12High + (prod1Low + prod2Low + s12Low);
894
895 if (Double.isNaN(result)) {
896 // either we have split infinite numbers or some coefficients were NaNs,
897 // just rely on the naive implementation and let IEEE754 handle this
898 result = a1 * b1 + a2 * b2;
899 }
900
901 return result;
902 }
903
904 /**
905 * Compute a linear combination accurately.
906 * <p>
907 * This method computes a<sub>1</sub>×b<sub>1</sub> +
908 * a<sub>2</sub>×b<sub>2</sub> + a<sub>3</sub>×b<sub>3</sub>
909 * to high accuracy. It does so by using specific multiplication and
910 * addition algorithms to preserve accuracy and reduce cancellation effects.
911 * It is based on the 2005 paper <a
912 * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
913 * Accurate Sum and Dot Product</a> by Takeshi Ogita,
914 * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
915 * </p>
916 * @param a1 first factor of the first term
917 * @param b1 second factor of the first term
918 * @param a2 first factor of the second term
919 * @param b2 second factor of the second term
920 * @param a3 first factor of the third term
921 * @param b3 second factor of the third term
922 * @return a<sub>1</sub>×b<sub>1</sub> +
923 * a<sub>2</sub>×b<sub>2</sub> + a<sub>3</sub>×b<sub>3</sub>
924 * @see #linearCombination(double, double, double, double)
925 * @see #linearCombination(double, double, double, double, double, double, double, double)
926 */
927 public static double linearCombination(final double a1, final double b1,
928 final double a2, final double b2,
929 final double a3, final double b3) {
930
931 // the code below is split in many additions/subtractions that may
932 // appear redundant. However, they should NOT be simplified, as they
933 // do use IEEE754 floating point arithmetic rounding properties.
934 // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
935 // The variables naming conventions are that xyzHigh contains the most significant
936 // bits of xyz and xyzLow contains its least significant bits. So theoretically
937 // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
938 // be represented in only one double precision number so we preserve two numbers
939 // to hold it as long as we can, combining the high and low order bits together
940 // only at the end, after cancellation may have occurred on high order bits
941
942 // split a1 and b1 as two 26 bits numbers
943 final double ca1 = SPLIT_FACTOR * a1;
944 final double a1High = ca1 - (ca1 - a1);
945 final double a1Low = a1 - a1High;
946 final double cb1 = SPLIT_FACTOR * b1;
947 final double b1High = cb1 - (cb1 - b1);
948 final double b1Low = b1 - b1High;
949
950 // accurate multiplication a1 * b1
951 final double prod1High = a1 * b1;
952 final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
953
954 // split a2 and b2 as two 26 bits numbers
955 final double ca2 = SPLIT_FACTOR * a2;
956 final double a2High = ca2 - (ca2 - a2);
957 final double a2Low = a2 - a2High;
958 final double cb2 = SPLIT_FACTOR * b2;
959 final double b2High = cb2 - (cb2 - b2);
960 final double b2Low = b2 - b2High;
961
962 // accurate multiplication a2 * b2
963 final double prod2High = a2 * b2;
964 final double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
965
966 // split a3 and b3 as two 26 bits numbers
967 final double ca3 = SPLIT_FACTOR * a3;
968 final double a3High = ca3 - (ca3 - a3);
969 final double a3Low = a3 - a3High;
970 final double cb3 = SPLIT_FACTOR * b3;
971 final double b3High = cb3 - (cb3 - b3);
972 final double b3Low = b3 - b3High;
973
974 // accurate multiplication a3 * b3
975 final double prod3High = a3 * b3;
976 final double prod3Low = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);
977
978 // accurate addition a1 * b1 + a2 * b2
979 final double s12High = prod1High + prod2High;
980 final double s12Prime = s12High - prod2High;
981 final double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
982
983 // accurate addition a1 * b1 + a2 * b2 + a3 * b3
984 final double s123High = s12High + prod3High;
985 final double s123Prime = s123High - prod3High;
986 final double s123Low = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime);
987
988 // final rounding, s123 may have suffered many cancellations, we try
989 // to recover some bits from the extra words we have saved up to now
990 double result = s123High + (prod1Low + prod2Low + prod3Low + s12Low + s123Low);
991
992 if (Double.isNaN(result)) {
993 // either we have split infinite numbers or some coefficients were NaNs,
994 // just rely on the naive implementation and let IEEE754 handle this
995 result = a1 * b1 + a2 * b2 + a3 * b3;
996 }
997
998 return result;
999 }
1000
1001 /**
1002 * Compute a linear combination accurately.
1003 * <p>
1004 * This method computes a<sub>1</sub>×b<sub>1</sub> +
1005 * a<sub>2</sub>×b<sub>2</sub> + a<sub>3</sub>×b<sub>3</sub> +
1006 * a<sub>4</sub>×b<sub>4</sub>
1007 * to high accuracy. It does so by using specific multiplication and
1008 * addition algorithms to preserve accuracy and reduce cancellation effects.
1009 * It is based on the 2005 paper <a
1010 * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
1011 * Accurate Sum and Dot Product</a> by Takeshi Ogita,
1012 * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
1013 * </p>
1014 * @param a1 first factor of the first term
1015 * @param b1 second factor of the first term
1016 * @param a2 first factor of the second term
1017 * @param b2 second factor of the second term
1018 * @param a3 first factor of the third term
1019 * @param b3 second factor of the third term
1020 * @param a4 first factor of the third term
1021 * @param b4 second factor of the third term
1022 * @return a<sub>1</sub>×b<sub>1</sub> +
1023 * a<sub>2</sub>×b<sub>2</sub> + a<sub>3</sub>×b<sub>3</sub> +
1024 * a<sub>4</sub>×b<sub>4</sub>
1025 * @see #linearCombination(double, double, double, double)
1026 * @see #linearCombination(double, double, double, double, double, double)
1027 */
1028 public static double linearCombination(final double a1, final double b1,
1029 final double a2, final double b2,
1030 final double a3, final double b3,
1031 final double a4, final double b4) {
1032
1033 // the code below is split in many additions/subtractions that may
1034 // appear redundant. However, they should NOT be simplified, as they
1035 // do use IEEE754 floating point arithmetic rounding properties.
1036 // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
1037 // The variables naming conventions are that xyzHigh contains the most significant
1038 // bits of xyz and xyzLow contains its least significant bits. So theoretically
1039 // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
1040 // be represented in only one double precision number so we preserve two numbers
1041 // to hold it as long as we can, combining the high and low order bits together
1042 // only at the end, after cancellation may have occurred on high order bits
1043
1044 // split a1 and b1 as two 26 bits numbers
1045 final double ca1 = SPLIT_FACTOR * a1;
1046 final double a1High = ca1 - (ca1 - a1);
1047 final double a1Low = a1 - a1High;
1048 final double cb1 = SPLIT_FACTOR * b1;
1049 final double b1High = cb1 - (cb1 - b1);
1050 final double b1Low = b1 - b1High;
1051
1052 // accurate multiplication a1 * b1
1053 final double prod1High = a1 * b1;
1054 final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
1055
1056 // split a2 and b2 as two 26 bits numbers
1057 final double ca2 = SPLIT_FACTOR * a2;
1058 final double a2High = ca2 - (ca2 - a2);
1059 final double a2Low = a2 - a2High;
1060 final double cb2 = SPLIT_FACTOR * b2;
1061 final double b2High = cb2 - (cb2 - b2);
1062 final double b2Low = b2 - b2High;
1063
1064 // accurate multiplication a2 * b2
1065 final double prod2High = a2 * b2;
1066 final double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
1067
1068 // split a3 and b3 as two 26 bits numbers
1069 final double ca3 = SPLIT_FACTOR * a3;
1070 final double a3High = ca3 - (ca3 - a3);
1071 final double a3Low = a3 - a3High;
1072 final double cb3 = SPLIT_FACTOR * b3;
1073 final double b3High = cb3 - (cb3 - b3);
1074 final double b3Low = b3 - b3High;
1075
1076 // accurate multiplication a3 * b3
1077 final double prod3High = a3 * b3;
1078 final double prod3Low = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);
1079
1080 // split a4 and b4 as two 26 bits numbers
1081 final double ca4 = SPLIT_FACTOR * a4;
1082 final double a4High = ca4 - (ca4 - a4);
1083 final double a4Low = a4 - a4High;
1084 final double cb4 = SPLIT_FACTOR * b4;
1085 final double b4High = cb4 - (cb4 - b4);
1086 final double b4Low = b4 - b4High;
1087
1088 // accurate multiplication a4 * b4
1089 final double prod4High = a4 * b4;
1090 final double prod4Low = a4Low * b4Low - (((prod4High - a4High * b4High) - a4Low * b4High) - a4High * b4Low);
1091
1092 // accurate addition a1 * b1 + a2 * b2
1093 final double s12High = prod1High + prod2High;
1094 final double s12Prime = s12High - prod2High;
1095 final double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
1096
1097 // accurate addition a1 * b1 + a2 * b2 + a3 * b3
1098 final double s123High = s12High + prod3High;
1099 final double s123Prime = s123High - prod3High;
1100 final double s123Low = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime);
1101
1102 // accurate addition a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4
1103 final double s1234High = s123High + prod4High;
1104 final double s1234Prime = s1234High - prod4High;
1105 final double s1234Low = (prod4High - (s1234High - s1234Prime)) + (s123High - s1234Prime);
1106
1107 // final rounding, s1234 may have suffered many cancellations, we try
1108 // to recover some bits from the extra words we have saved up to now
1109 double result = s1234High + (prod1Low + prod2Low + prod3Low + prod4Low + s12Low + s123Low + s1234Low);
1110
1111 if (Double.isNaN(result)) {
1112 // either we have split infinite numbers or some coefficients were NaNs,
1113 // just rely on the naive implementation and let IEEE754 handle this
1114 result = a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4;
1115 }
1116
1117 return result;
1118 }
1119
1120 /**
1121 * Returns true iff both arguments are null or have same dimensions and all
1122 * their elements are equal as defined by
1123 * {@link Precision#equals(float,float)}.
1124 *
1125 * @param x first array
1126 * @param y second array
1127 * @return true if the values are both null or have same dimension
1128 * and equal elements.
1129 */
1130 public static boolean equals(float[] x, float[] y) {
1131 if ((x == null) || (y == null)) {
1132 return !((x == null) ^ (y == null));
1133 }
1134 if (x.length != y.length) {
1135 return false;
1136 }
1137 for (int i = 0; i < x.length; ++i) {
1138 if (!Precision.equals(x[i], y[i])) {
1139 return false;
1140 }
1141 }
1142 return true;
1143 }
1144
1145 /**
1146 * Returns true iff both arguments are null or have same dimensions and all
1147 * their elements are equal as defined by
1148 * {@link Precision#equalsIncludingNaN(double,double) this method}.
1149 *
1150 * @param x first array
1151 * @param y second array
1152 * @return true if the values are both null or have same dimension and
1153 * equal elements
1154 * @since 2.2
1155 */
1156 public static boolean equalsIncludingNaN(float[] x, float[] y) {
1157 if ((x == null) || (y == null)) {
1158 return !((x == null) ^ (y == null));
1159 }
1160 if (x.length != y.length) {
1161 return false;
1162 }
1163 for (int i = 0; i < x.length; ++i) {
1164 if (!Precision.equalsIncludingNaN(x[i], y[i])) {
1165 return false;
1166 }
1167 }
1168 return true;
1169 }
1170
1171 /**
1172 * Returns {@code true} iff both arguments are {@code null} or have same
1173 * dimensions and all their elements are equal as defined by
1174 * {@link Precision#equals(double,double)}.
1175 *
1176 * @param x First array.
1177 * @param y Second array.
1178 * @return {@code true} if the values are both {@code null} or have same
1179 * dimension and equal elements.
1180 */
1181 public static boolean equals(double[] x, double[] y) {
1182 if ((x == null) || (y == null)) {
1183 return !((x == null) ^ (y == null));
1184 }
1185 if (x.length != y.length) {
1186 return false;
1187 }
1188 for (int i = 0; i < x.length; ++i) {
1189 if (!Precision.equals(x[i], y[i])) {
1190 return false;
1191 }
1192 }
1193 return true;
1194 }
1195
1196 /**
1197 * Returns {@code true} iff both arguments are {@code null} or have same
1198 * dimensions and all their elements are equal as defined by
1199 * {@link Precision#equalsIncludingNaN(double,double) this method}.
1200 *
1201 * @param x First array.
1202 * @param y Second array.
1203 * @return {@code true} if the values are both {@code null} or have same
1204 * dimension and equal elements.
1205 * @since 2.2
1206 */
1207 public static boolean equalsIncludingNaN(double[] x, double[] y) {
1208 if ((x == null) || (y == null)) {
1209 return !((x == null) ^ (y == null));
1210 }
1211 if (x.length != y.length) {
1212 return false;
1213 }
1214 for (int i = 0; i < x.length; ++i) {
1215 if (!Precision.equalsIncludingNaN(x[i], y[i])) {
1216 return false;
1217 }
1218 }
1219 return true;
1220 }
1221
1222 /**
1223 * Normalizes an array to make it sum to a specified value.
1224 * Returns the result of the transformation <pre>
1225 * x |-> x * normalizedSum / sum
1226 * </pre>
1227 * applied to each non-NaN element x of the input array, where sum is the
1228 * sum of the non-NaN entries in the input array.</p>
1229 *
1230 * <p>Throws IllegalArgumentException if {@code normalizedSum} is infinite
1231 * or NaN and ArithmeticException if the input array contains any infinite elements
1232 * or sums to 0.</p>
1233 *
1234 * <p>Ignores (i.e., copies unchanged to the output array) NaNs in the input array.</p>
1235 *
1236 * @param values Input array to be normalized
1237 * @param normalizedSum Target sum for the normalized array
1238 * @return the normalized array.
1239 * @throws MathArithmeticException if the input array contains infinite
1240 * elements or sums to zero.
1241 * @throws MathIllegalArgumentException if the target sum is infinite or {@code NaN}.
1242 * @since 2.1
1243 */
1244 public static double[] normalizeArray(double[] values, double normalizedSum)
1245 throws MathIllegalArgumentException, MathArithmeticException {
1246 if (Double.isInfinite(normalizedSum)) {
1247 throw new MathIllegalArgumentException(LocalizedFormats.NORMALIZE_INFINITE);
1248 }
1249 if (Double.isNaN(normalizedSum)) {
1250 throw new MathIllegalArgumentException(LocalizedFormats.NORMALIZE_NAN);
1251 }
1252 double sum = 0d;
1253 final int len = values.length;
1254 double[] out = new double[len];
1255 for (int i = 0; i < len; i++) {
1256 if (Double.isInfinite(values[i])) {
1257 throw new MathIllegalArgumentException(LocalizedFormats.INFINITE_ARRAY_ELEMENT, values[i], i);
1258 }
1259 if (!Double.isNaN(values[i])) {
1260 sum += values[i];
1261 }
1262 }
1263 if (sum == 0) {
1264 throw new MathArithmeticException(LocalizedFormats.ARRAY_SUMS_TO_ZERO);
1265 }
1266 for (int i = 0; i < len; i++) {
1267 if (Double.isNaN(values[i])) {
1268 out[i] = Double.NaN;
1269 } else {
1270 out[i] = values[i] * normalizedSum / sum;
1271 }
1272 }
1273 return out;
1274 }
1275 }