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