001/* A library for mathematical operations on arrays of complex numbers.
002
003 Copyright (c) 1998-2013 The Regents of the University of California.
004 All rights reserved.
005
006 Permission is hereby granted, without written agreement and without
007 license or royalty fees, to use, copy, modify, and distribute this
008 software and its documentation for any purpose, provided that the above
009 copyright notice and the following two paragraphs appear in all copies
010 of this software.
011
012 IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
013 FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
014 ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
015 THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF
016 SUCH DAMAGE.
017
018 THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
019 INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
020 MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE
021 PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
022 CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES,
023 ENHANCEMENTS, OR MODIFICATIONS.
024
025 PT_COPYRIGHT_VERSION_2
026 COPYRIGHTENDKEY
027
028 */
029package ptolemy.math;
030
031///////////////////////////////////////////////////////////////////
032//// ComplexArrayMath
033
034/**
035 This class a provides a library for mathematical operations on arrays of
036 complex numbers, in particular arrays of instances of class
037 ptolemy.math.Complex.
038 Unless explicitly noted otherwise, all array arguments are assumed to be
039 non-null. If a null array is passed to a method, a NullPointerException
040 will be thrown in the method or called methods.
041
042 @author Albert Chen, William Wu, Edward A. Lee, Jeff Tsay
043 @version $Id$
044 @since Ptolemy II 0.3
045 @Pt.ProposedRating Yellow (ctsay)
046 @Pt.AcceptedRating Yellow (ctsay)
047 */
048public class ComplexArrayMath {
049    // Protected constructor prevents construction of this class.
050    protected ComplexArrayMath() {
051    }
052
053    /** Return a new array that is constructed from <i>array</i> by
054     *  adding the complex number <i>z</i> to every element of <i>array</i>.
055     *
056     *  @param array An array of complex numbers.
057     *  @param z The complex number to be added to each element of <i>array</i>.
058     *  @return A new array of complex numbers equal to <i>array</i> with
059     *  <i>z</i> added to each element.
060     */
061    public static final Complex[] add(Complex[] array, Complex z) {
062        Complex[] result = new Complex[array.length];
063
064        for (int i = 0; i < array.length; i++) {
065            result[i] = array[i].add(z);
066        }
067
068        return result;
069    }
070
071    /** Return a new array that is the element-by-element sum of the two
072     *  input arrays.
073     *  If the sizes of both arrays are 0, return a new array of size 0.
074     *
075     *  @param array1 The first array.
076     *  @param array2 The second array.
077     *  @return A new array that is the element-by-element sum of the two
078     *  input arrays.
079     *  @exception IllegalArgumentException If the arrays are not of the same
080     *   length.
081     */
082    public static final Complex[] add(Complex[] array1, Complex[] array2) {
083        int length = _commonLength(array1, array2, "ComplexArrayMath.add");
084        Complex[] returnValue = new Complex[length];
085
086        for (int i = 0; i < length; i++) {
087            returnValue[i] = array1[i].add(array2[i]);
088        }
089
090        return returnValue;
091    }
092
093    /** Return a new array that is the result of appending <i>array2</i>
094     *  to the end of <i>array1</i>. This method simply calls
095     *  append(array1, 0, array1.length, array2, 0, array2.length)
096     *
097     *  @param array1 The first array of complex numbers.
098     *  @param array2 The second array of complex numbers.
099     *  @return A new array formed by appending <i>array2</i>
100     *  to the end of <i>array1</i>.
101     */
102    public static final Complex[] append(Complex[] array1, Complex[] array2) {
103        return append(array1, 0, array1.length, array2, 0, array2.length);
104    }
105
106    /** Return a new array that is the result of appending
107     *  <i>length2</i> elements of <i>array2</i>, starting from the
108     *  <i>idx2</i>th element, to <i>length1</i> elements of
109     *  <i>array1</i>, starting from the <i>idx1</i>th element.
110     *
111     *  Appending empty arrays is supported. In that case, the corresponding
112     *  idx may be any number. Allow System.arraycopy() to throw array access
113     *  exceptions if idx .. idx + length - 1 are not all valid array indices,
114     *  for both of the arrays.
115     *
116     *  @param array1 The first array of complex numbers.
117     *  @param idx1 The starting index for array1.
118     *  @param length1 The number of elements of array1 to use.
119     *  @param array2 The second array of complex numbers, which is appended.
120     *  @param idx2 The starting index for array2.
121     *  @param length2 The number of elements of array2 to append.
122     *  @return A new array of Complex.
123     */
124    public static final Complex[] append(Complex[] array1, int idx1,
125            int length1, Complex[] array2, int idx2, int length2) {
126        Complex[] returnValue = new Complex[length1 + length2];
127
128        if (length1 > 0) {
129            System.arraycopy(array1, idx1, returnValue, 0, length1);
130        }
131
132        if (length2 > 0) {
133            System.arraycopy(array2, idx2, returnValue, length1, length2);
134        }
135
136        return returnValue;
137    }
138
139    /** Return a new array that is formed by applying an instance of a
140     *  ComplexBinaryOperation to each element in the input array,
141     *  using <i>z</i> as the left argument to <i>op</i> in all cases and
142     *  the array elements as the right arguments (op.operate(z, array[i])).
143     *  If the length of the array is 0, return a new array of length 0.
144     *
145     *  @param op The complex binary operation to be applied to the given
146     *  complex number and complex array.
147     *  @param z The complex number that is the first argument to <i>op</i>.
148     *  @param array The array of complex numbers that is the second argument
149     *  to <i>op</i>.
150     *  @return A new array containing elements equal to
151     *  (op.operate(z, array[i])).
152     */
153    public static final Complex[] applyBinaryOperation(
154            ComplexBinaryOperation op, final Complex z, final Complex[] array) {
155        int length = array.length;
156        Complex[] returnValue = new Complex[length];
157
158        for (int i = 0; i < length; i++) {
159            returnValue[i] = op.operate(z, array[i]);
160        }
161
162        return returnValue;
163    }
164
165    /** Return a new array that is formed by applying an instance of a
166     *  ComplexBinaryOperation to each element in the input array,
167     *  using <i>z</i> as the right operand in all cases and the array elements
168     *  as the left operands (op.operate(array[i], z)).
169     *  If the length of the array is 0, return a new array of length 0.
170     *
171     *  @param op The complex binary operation to be applied to the given
172     *  complex number and complex array.
173     *  @param z The complex number that is the second argument to <i>op</i>.
174     *  @param array The array of complex numbers that is the first argument
175     *  to <i>op</i>.
176     *  @return A new array containing elements equal to
177     *  (op.operate(array[i], z)).
178     */
179    public static final Complex[] applyBinaryOperation(
180            ComplexBinaryOperation op, final Complex[] array, final Complex z) {
181        int length = array.length;
182        Complex[] returnValue = new Complex[length];
183
184        for (int i = 0; i < length; i++) {
185            returnValue[i] = op.operate(array[i], z);
186        }
187
188        return returnValue;
189    }
190
191    /** Return a new array that is formed by applying an instance of a
192     *  ComplexBinaryOperation to the two arrays, element by element,
193     *  using the elements of the first array as the left operands and the
194     *  elements of the second array as the right operands.
195     *  (op.operate(array[i], array2[i])).
196     *  If the lengths of both arrays are 0, return a new array of length 0.
197     *
198     *  @param op The complex binary operation to be applied to each pair
199     *  of corresponding elements of each complex array.
200     *  @param array1 The first array.
201     *  @param array2 The second array.
202     *  @return A new array that with elements equal to
203     *  (op.operate(array[i], array2[i])).
204     *  @exception IllegalArgumentException If the arrays are not of the same
205     *   length.
206     */
207    public static final Complex[] applyBinaryOperation(
208            ComplexBinaryOperation op, final Complex[] array1,
209            final Complex[] array2) {
210        int length = _commonLength(array1, array2,
211                "ComplexArrayMath.applyBinaryOperation");
212        Complex[] returnValue = new Complex[length];
213
214        for (int i = 0; i < length; i++) {
215            returnValue[i] = op.operate(array1[i], array2[i]);
216        }
217
218        return returnValue;
219    }
220
221    /** Return a new array that is formed by applying an instance of a
222     *  ComplexUnaryOperation to each element in the input array
223     *  (op.operate(array[i])).
224     *  If the length of the array is 0, return a new array of length 0.
225     *
226     *  @param op The complex unary operation to be applied to each
227     *  element of the given complex array.
228     *  @param array An array of complex numbers.
229     *  @return A new array of complex numbers with each element
230     *  equal to (op.operate(array[i])).
231     */
232    public static final Complex[] applyUnaryOperation(
233            final ComplexUnaryOperation op, final Complex[] array) {
234        int length = array.length;
235        Complex[] returnValue = new Complex[length];
236
237        for (int i = 0; i < length; i++) {
238            returnValue[i] = op.operate(array[i]);
239        }
240
241        return returnValue;
242    }
243
244    /** Return a new array of complex numbers that is formed by taking the
245     *  complex-conjugate of each element in the argument array.
246     *  If the argument has length 0, return a new array of complex numbers,
247     *  with length 0.
248     *
249     *  @param array The given array of complex numbers.
250     *  @return A new array of complex numbers formed by taking the
251     *  complex-conjugate of each element in the argument array.
252     */
253    public static final Complex[] conjugate(Complex[] array) {
254        Complex[] result = new Complex[array.length];
255
256        for (int i = array.length - 1; i >= 0; i--) {
257            result[i] = array[i].conjugate();
258        }
259
260        return result;
261    }
262
263    /** Return a new array that is the element-by-element division of
264     *  the first array by the second array.
265     *  If the sizes of both arrays are 0, return a new array of size 0.
266     *
267     *  @param array1 The first array of complex numbers.
268     *  @param array2 The second array of complex numbers.
269     *  @return A new array of complex numbers equal to the
270     *  element-by-element division of the first array by the second array.
271     *  @exception IllegalArgumentException If the arrays are not of the same
272     *   length.
273     */
274    public static final Complex[] divideElements(Complex[] array1,
275            Complex[] array2) {
276        int length = _commonLength(array1, array2,
277                "ComplexArrayMath.divideElements");
278        Complex[] returnValue = new Complex[length];
279
280        for (int i = 0; i < length; i++) {
281            returnValue[i] = array1[i].divide(array2[i]);
282        }
283
284        return returnValue;
285    }
286
287    /** Return a new array that is the result of dividing each element of
288     *  the given array by the given value.
289     *
290     *  @param array An array of complex numbers.
291     *  @param divisor The number by which to divide each element of the array.
292     *  @return A new array of complex numbers that is the result of dividing
293     *  each element of <i>array</i> by <i>divisor</i>.
294     */
295    public static final Complex[] divide(Complex[] array, Complex divisor) {
296        int length = array.length;
297        Complex[] returnValue = new Complex[length];
298
299        for (int i = 0; i < length; i++) {
300            returnValue[i] = array[i].divide(divisor);
301        }
302
303        return returnValue;
304    }
305
306    /** Return a complex number that is the dot product of the two argument
307     *  arrays. The dot product is computed by the sum of the
308     *  element-by-element products of <i>array1</i> and the complex
309     *  conjugates of <i>array2</i>.
310     *  If the size of each array is 0, return Complex.ZERO.
311     *
312     *  @param array1 The first array of complex numbers.
313     *  @param array2 The second array of complex numbers.
314     *  @return A new array of complex numbers equal to the
315     *  dot product of the two given arrays.
316     *  @exception IllegalArgumentException If the arrays are not of the same
317     *   length.
318     */
319    public static final Complex dotProduct(final Complex[] array1,
320            final Complex[] array2) {
321        int length = _commonLength(array1, array2,
322                "ComplexArrayMath.dotProduct");
323        Complex returnValue = Complex.ZERO;
324
325        for (int i = 0; i < length; i++) {
326            returnValue = returnValue
327                    .add(array1[i].multiply(array2[i].conjugate()));
328        }
329
330        return returnValue;
331    }
332
333    /** Return a new array of Complex numbers using two arrays for the
334     *  real and imaginary parts. If <i>realPart</i> is null, it is
335     *  treated as if it were an array of zeros and an array of
336     *  complex numbers that are purely imaginary is constructed. If
337     *  <i>imagPart</i> is null, it is treated as if it were an array
338     *  of zeroes and an array of complex numbers that are purely real
339     *  is constructed.  If both arrays are of length 0, or one array
340     *  is null and the other array is of length 0, return a new array
341     *  of complex numbers with length 0.  If both arrays are null,
342     *  allow a NullPointerException to be thrown by the array access
343     *  code.
344     *
345     *  @param realPart An array of doubles, used for the real parts.
346     *  @param imagPart An array of doubles, used for the imaginary parts.
347     *  @return A new array of complex numbers each containing real
348     *  and imaginary parts equal to the elements in <i>realPart</i>
349     *  and <i>imagPart</i>, respectively.
350     */
351    public static final Complex[] formComplexArray(final double[] realPart,
352            final double[] imagPart) {
353        Complex[] returnValue;
354        int size;
355
356        if (realPart != null && imagPart != null) {
357            size = DoubleArrayMath._commonLength(realPart, imagPart,
358                    "ComplexArrayMath.formComplexArray");
359            returnValue = new Complex[size];
360
361            for (int i = 0; i < size; i++) {
362                returnValue[i] = new Complex(realPart[i], imagPart[i]);
363            }
364        } else if (realPart == null) {
365            // NullPointerException will be thrown here if both
366            // arrays are null.
367            size = imagPart.length;
368
369            returnValue = new Complex[size];
370
371            for (int i = 0; i < size; i++) {
372                returnValue[i] = new Complex(0.0, imagPart[i]);
373            }
374        } else { // imagPart == null
375            size = realPart.length;
376
377            returnValue = new Complex[size];
378
379            for (int i = 0; i < size; i++) {
380                returnValue[i] = new Complex(realPart[i], 0.0);
381            }
382        }
383
384        return returnValue;
385    }
386
387    /** Return a new array of doubles with the imaginary parts of the array of
388     *  complex numbers.
389     *
390     *  @param x The input array of complex numbers.
391     *  @return A new array of doubles with the imaginary parts of the array of
392     *  complex numbers.
393     */
394    public static final double[] imagParts(final Complex[] x) {
395        int size = x.length;
396
397        double[] returnValue = new double[size];
398
399        for (int i = 0; i < size; i++) {
400            returnValue[i] = x[i].imag;
401        }
402
403        return returnValue;
404    }
405
406    /** Return a double that is the L2-norm of the array.  If the
407     *  length of the array is zero, return 0.0.
408     *
409     *  @param array The given array of complex numbers.
410     *  @return A double that is the L2-norm of <i>array</i>.
411     */
412    public static final double l2norm(Complex[] array) {
413        return Math.sqrt(l2normSquared(array));
414    }
415
416    /** Return a double that is the sum of the squared magnitudes of
417     *  the elements of <i>array</i>. This is equal to the square of the
418     *  L2-norm of <i>array</i>. If the length of the array is zero,
419     *  return 0.0.
420     *
421     *  @param array The given array of complex numbers.
422     *  @return A double that is square of the L2-norm of <i>array</i>.
423     */
424    public static final double l2normSquared(Complex[] array) {
425        int length = array.length;
426
427        if (length <= 0) {
428            return 0.0;
429        }
430
431        double returnValue = 0.0;
432
433        for (int i = 0; i < length; i++) {
434            returnValue += array[i].magnitudeSquared();
435        }
436
437        return returnValue;
438    }
439
440    /** Return a new array that is a copy of the first argument except
441     *  that the elements are limited to lie within the specified range.
442     *  The specified range is given by two complex numbers, <i>bottom</i>
443     *  and <i>top</i>, where both the real and imaginary parts of
444     *  <i>bottom</i> are expected to be less than the real and imaginary
445     *  parts of <i>top</i>. Thus, <i>bottom</i> and <i>top</i> define a
446     *  rectangle in the complex plane, with <i>bottom</i> at the lower
447     *  left and <i>top</i> at the upper right.
448     *  If any value in the array is infinite
449     *  then it is replaced by the corresponding real or imaginary part
450     *  of <i>top</i> or <i>bottom</i>, depending on the sign of infinity.
451     *  If any value is NaN (not a number), then the result will be NaN.
452     *  To leave either the bottom or the top unconstrained,
453     *  specify Complex.NEGATIVE_INFINITY or Complex.POSITIVE_INFINITY.
454     *  If the length of the array is 0, return a new array of length 0.
455     *
456     *  @param array An array of complex numbers.
457     *  @param bottom The bottom limit.
458     *  @param top The top limit.
459     *  @return A new array with values in the rectangle defined by
460     *   <i>bottom</i> and <i>top</i>.
461     *  @exception IllegalArgumentException If <i>bottom</i> has either a
462     *   real or imaginary part larger than the corresponding part of
463     *   <i>top</i>.
464     */
465    public static final Complex[] limit(Complex[] array, Complex bottom,
466            Complex top) throws IllegalArgumentException {
467        Complex[] returnValue = new Complex[array.length];
468
469        // Check validity of the rectangle.
470        if (bottom.real > top.real || bottom.imag > top.imag) {
471            throw new IllegalArgumentException(
472                    "Complex.limit requires that bottom lie below and "
473                            + "to the left of top.");
474        }
475
476        for (int i = 0; i < array.length; i++) {
477            double realPart;
478            double imagPart;
479
480            // NOTE: Assume here that if array[i].real is NaN, then
481            // this test returns false.
482            if (array[i].real > top.real) {
483                realPart = top.real;
484            } else if (array[i].real < bottom.real) {
485                realPart = bottom.real;
486            } else {
487                realPart = array[i].real;
488            }
489
490            // NOTE: Assume here that if array[i].imag is NaN, then
491            // this test returns false.
492            if (array[i].imag > top.imag) {
493                imagPart = top.imag;
494            } else if (array[i].imag < bottom.imag) {
495                imagPart = bottom.imag;
496            } else {
497                imagPart = array[i].imag;
498            }
499
500            returnValue[i] = new Complex(realPart, imagPart);
501        }
502
503        return returnValue;
504    }
505
506    /** Return a new array of doubles containing the magnitudes of the elements
507     *  of the specified array of complex numbers.
508     *
509     *  @param array The given array of complex numbers.
510     *  @return A new array of doubles containing the magnitudes of
511     *  the elements of the specified array of complex numbers.
512     */
513    public static final double[] magnitude(Complex[] array) {
514        double[] mags = new double[array.length];
515
516        for (int i = array.length - 1; i >= 0; i--) {
517            mags[i] = array[i].magnitude();
518        }
519
520        return mags;
521    }
522
523    /** Return a new array that is the element-by-element multiplication of
524     *  the two input arrays.
525     *  If the size of each array is 0, return a new array of size 0.
526     *
527     *  @param array1 The first array of complex numbers.
528     *  @param array2 The second array of complex numbers.
529     *  @return A new array of complex numbers equal to the
530     *  element-by-element multiplication of the two given arrays.
531     *  @exception IllegalArgumentException If the arrays are not of the same
532     *   length.
533     */
534    public static final Complex[] multiply(Complex[] array1, Complex[] array2) {
535        int length = _commonLength(array1, array2, "ComplexArrayMath.multiply");
536        Complex[] returnValue = new Complex[length];
537
538        for (int i = 0; i < length; i++) {
539            returnValue[i] = array1[i].multiply(array2[i]);
540        }
541
542        return returnValue;
543    }
544
545    /** Return a new array that is constructed from the argument by
546     *  multiplying each element in <i>array</i> by the second argument,
547     *  which is a complex number.
548     *  If the size of the array is 0, return a new array of size 0.
549     *
550     *  @param array An array of complex numbers.
551     *  @param factor The complex number by which each element of <i>array</i>
552     *  is multiplied by.
553     *  @return A new array of complex numbers formed by
554     *  multiplying each element in <i>array</i> by <i>factor</i>.
555     */
556    public static final Complex[] multiply(Complex[] array, Complex factor) {
557        int length = array.length;
558        Complex[] returnValue = new Complex[length];
559
560        for (int i = 0; i < length; i++) {
561            returnValue[i] = array[i].multiply(factor);
562        }
563
564        return returnValue;
565    }
566
567    /** Return a new array that is the formed from the additive inverse of each
568     *  element of the input array (-array[i]).
569     *
570     *  @param array The given array of complex numbers.
571     *  @return A new array of complex numbers formed from the
572     *  additive inverse of each element of the input array (-array[i]).
573     */
574    public static final Complex[] negative(final Complex[] array) {
575        int length = array.length;
576        Complex[] returnValue = new Complex[length];
577
578        for (int i = 0; i < length; i++) {
579            returnValue[i] = new Complex(-array[i].real, -array[i].imag);
580        }
581
582        return returnValue;
583    }
584
585    /** Return a new array of Complex numbers that is formed by padding the
586     *  middle of the array with 0's. If the length of the
587     *  input array is odd, the sample with index ceil(L/2) will be
588     *  repeated in the output array, where L is the length of the input array.
589     *  If the lengths of the input and output arrays are equal, return
590     *  a copy of the input array.
591     *  This method is useful for preparing data for an IFFT.
592     *
593     *  @param array An array of complex numbers.
594     *  @param newLength The desired length of the returned array.
595     *  @return A new array of complex numbers formed by padding the
596     *  middle of <i>array</i> with 0's.
597     */
598    public static final Complex[] padMiddle(final Complex[] array,
599            final int newLength) {
600        int length = array.length;
601
602        int entriesNeeded = newLength - length;
603
604        if (entriesNeeded < 0) {
605            throw new IllegalArgumentException("ptolemy.math."
606                    + "ComplexArrayMath.padMiddle() : newLength must be "
607                    + ">= length of array.");
608        } else if (entriesNeeded == 0) {
609            return resize(array, newLength); // allocates a new array
610        }
611
612        double halfLength = length * 0.5;
613        int halfLengthFloor = (int) Math.floor(halfLength);
614        int halfLengthCeil = (int) Math.ceil(halfLength);
615        Complex[] returnValue = new Complex[newLength];
616
617        System.arraycopy(array, 0, returnValue, 0, halfLengthCeil);
618
619        System.arraycopy(array, halfLengthFloor, returnValue,
620                newLength - halfLengthCeil, halfLengthCeil);
621
622        for (int i = halfLengthCeil; i < newLength - halfLengthCeil; i++) {
623            returnValue[i] = Complex.ZERO;
624        }
625
626        return returnValue;
627    }
628
629    /** Return a new array containing the angles of the elements of the
630     *  specified complex array.
631     *
632     *  @param array An array of complex numbers.
633     *  @return An array of angles in the range of
634     *  <em>-pi</em> to <em>pi</em>.
635     */
636    public static final double[] phase(Complex[] array) {
637        double[] angles = new double[array.length];
638
639        for (int i = array.length - 1; i >= 0; i--) {
640            angles[i] = array[i].angle();
641        }
642
643        return angles;
644    }
645
646    /** Given the roots of a polynomial, return a polynomial that has
647     *  has such roots.  If the roots are
648     *  [<em>r</em><sub>0</sub>, ..., <em>r</em><sub>N-1</sub>],
649     *  then the polynomial is given by
650     *  [<em>a</em><sub>0</sub>, ..., <em>a</em><sub>N</sub>], where
651     *  <p>
652     *  <em>a</em><sub>0</sub> +
653     *  <em>a</em><sub>1</sub><em>z</em><sup>-1</sup> + ... +
654     *  <em>a</em><sub>N</sub><em>z</em><sup>-N</sup> =
655     *  (1 - <em>r</em><sub>0</sub><em>z</em><sup>-1</sup>)
656     *  (1 - <em>r</em><sub>1</sub><em>z</em><sup>-1</sup>) ...
657     *  (1 - <em>r</em><sub>N-1</sub><em>z</em><sup>-1</sup>).
658     *  The returned polynomial will always be monic, meaning that
659     *  <em>a</em><sub>0</sub> = 1.
660     *
661     *  @param roots An array of roots of a polynomial.
662     *  @return A new array representing a monic polynomial with the given
663     *  roots.
664     */
665    public static final Complex[] polynomial(final Complex[] roots) {
666        if (roots.length <= 1) {
667            return new Complex[] { Complex.ONE };
668        }
669
670        Complex[] result = new Complex[2];
671        result[0] = Complex.ONE;
672
673        if (roots.length >= 1) {
674            result[1] = roots[0].negate();
675
676            if (roots.length > 1) {
677                for (int i = 1; i < roots.length; i++) {
678                    Complex[] factor = { new Complex(1), roots[i].negate() };
679                    result = SignalProcessing.convolve(result, factor);
680                }
681            }
682        }
683
684        return result;
685    }
686
687    /** Return a new array of complex numbers that is formed by raising each
688     *  element to the specified exponent, a double.
689     *  If the size of the array is 0, return a new array of size 0.
690     *
691     *  @param array The input array of complex numbers.
692     *  @param exponent A double, which is the exponent.
693     *  @return A new array that is formed by raising each
694     *  element of <i>array</i> to the <i>element</i>th power.
695     */
696    public static final Complex[] pow(final Complex[] array,
697            final double exponent) {
698        int length = array.length;
699        Complex[] returnValue = new Complex[length];
700
701        for (int i = 0; i < length; i++) {
702            returnValue[i] = array[i].pow(exponent);
703        }
704
705        return returnValue;
706    }
707
708    /** Return the product of the elements in the array.
709     *  If there are no elements in the array, return a Complex number
710     *  with value zero.
711     *
712     *  @param array An array of complex numbers.
713     *  @return A new complex number equal to the product of
714     *  the elements in the array.
715     */
716    public static final Complex product(final Complex[] array) {
717        if (array.length == 0) {
718            return Complex.ZERO;
719        }
720
721        double real = 1.0;
722        double imag = 0.0;
723
724        for (Complex element : array) {
725            double tmp = real * element.real - imag * element.imag;
726            imag = real * element.imag + imag * element.real;
727            real = tmp;
728        }
729
730        return new Complex(real, imag);
731    }
732
733    /** Return a new array of doubles that includes the real
734     *  parts of the array of complex numbers.
735     *
736     *  @param x An array of complex numbers
737     *  @return A new array of doubles that includes the real
738     *  parts of the array of complex numbers.
739     */
740    public static final double[] realParts(final Complex[] x) {
741        int size = x.length;
742
743        double[] returnValue = new double[size];
744
745        for (int i = 0; i < size; i++) {
746            returnValue[i] = x[i].real;
747        }
748
749        return returnValue;
750    }
751
752    /** Return a new array of length <i>newLength</i> that is formed by
753     *  either truncating or padding the input array.
754     *  This method simply calls :
755     *  resize(array, newLength, 0)
756     *
757     *  @param array An array of complex numbers.
758     *  @param newLength The desired size of the output array.
759     *  @return A new array of length <i>newLength</i> that is formed by
760     *  either truncating or padding <i>array</i>.
761     */
762    public static final Complex[] resize(final Complex[] array,
763            final int newLength) {
764        return resize(array, newLength, 0);
765    }
766
767    /** Return a new array of length <i>newLength</i> that is formed by
768     *  either truncating or padding the input array.
769     *  Elements from the input array are copied to the output array,
770     *  starting from array[startIdx] until one of the following conditions
771     *  is met :
772     *  1) The input array has no more elements to copy.
773     *  2) The output array has been completely filled.
774     *  <i>startIdx</i> must index a valid entry in <i>array</i> unless
775     *  the input array is of zero length or the output array is of zero length.
776     *  If case 1) is met, the remainder of the output array is filled with
777     *  new complex numbers with value 0.
778     *  Copying here means shallow copying, i.e. pointers to Complex objects
779     *  are copied instead of allocation of new copies. This works because
780     *  Complex objects are immutable.
781     *
782     *  @param array An array of complex numbers.
783     *  @param newLength The desired size of the output array.
784     *  @param startIdx The starting index for the input array.
785     *  @return A new array of length <i>newLength</i> that is formed by
786     *  either truncating or padding <i>array</i>.
787     */
788    public static final Complex[] resize(final Complex[] array,
789            final int newLength, final int startIdx) {
790        Complex[] returnValue = new Complex[newLength];
791        int copySize = Math.min(newLength, array.length - startIdx);
792
793        if (startIdx >= array.length && copySize > 0) {
794            throw new IllegalArgumentException("resize():  the start index '"
795                    + startIdx + "' is greater than equal to the array length '"
796                    + array.length + "' and the number of items to be copied '"
797                    + copySize + "' is greater than zero.");
798        }
799
800        if (copySize > 0) {
801            System.arraycopy(array, startIdx, returnValue, 0, copySize);
802        }
803
804        for (int i = copySize; i < newLength; i++) {
805            returnValue[i] = Complex.ZERO;
806        }
807
808        return returnValue;
809    }
810
811    /** Return a new array that is constructed from the argument by
812     *  scaling each element in <i>array</i> by <i>factor</i>, which is a
813     *  complex number. If the array argument is of length 0, return a new
814     *  array of length 0.
815     *
816     *  @param array An array of complex numbers.
817     *  @param factor A complex number used to multiply each element
818     *  of <i>array</i> by.
819     *  @return A new array of complex numbers that is constructed from
820     *  the argument by scaling each element in <i>array</i> by <i>factor</i>.
821     */
822    public static final Complex[] scale(final Complex[] array,
823            final Complex factor) {
824        int len = array.length;
825        Complex[] returnValue = new Complex[len];
826
827        for (int i = 0; i < len; i++) {
828            returnValue[i] = array[i].multiply(factor);
829        }
830
831        return returnValue;
832    }
833
834    /** Return a new array that is constructed from the argument by
835     *  scaling each element in the array by factor, which is a
836     *  double. If the array argument is of length 0, return a new
837     *  array of length 0.
838     *
839     *  @param array An array of complex numbers.
840     *  @param factor A double used to multiply each element
841     *  of <i>array</i> by.
842     *  @return A new array of complex numbers that is constructed from
843     *  the argument by scaling each element in <i>array</i> by <i>factor</i>.
844     */
845    public static final Complex[] scale(final Complex[] array,
846            final double factor) {
847        int len = array.length;
848        Complex[] returnValue = new Complex[len];
849
850        for (int i = 0; i < len; i++) {
851            returnValue[i] = array[i].scale(factor);
852        }
853
854        return returnValue;
855    }
856
857    /** Return a new array that is constructed by subtracting the complex
858     *  number z from every element in the given array. If the array argument
859     *  is of length 0, return a new array of length 0.
860     *
861     *  @param array An array of complex numbers.
862     *  @param z A complex number subtracted from each element of <i>array</i>.
863     *  @return A new array that is constructed by subtracting <i>z</i>
864     *  from every element in <i>array</i>.
865     */
866    public static final Complex[] subtract(final Complex[] array,
867            final Complex z) {
868        Complex[] result = new Complex[array.length];
869
870        for (int i = array.length - 1; i >= 0; i--) {
871            result[i] = array[i].subtract(z);
872        }
873
874        return result;
875    }
876
877    /** Return a new array that is the element-by-element
878     *  subtraction of the second array from the first array.
879     *  If the size of each array is 0, return a new array of size 0.
880     *
881     *  @param array1 An array of complex numbers from which to subtract.
882     *  @param array2 An array of complex numbers to subtract.
883     *  @return A new array of complex numbers equal to the element-by-element
884     *  subtraction of the second array from the first array.
885     *  @exception IllegalArgumentException If the arrays are not of the same
886     *   length.
887     */
888    public static final Complex[] subtract(Complex[] array1,
889            final Complex[] array2) {
890        int length = _commonLength(array1, array2, "ComplexArrayMath.subtract");
891        Complex[] result = new Complex[length];
892
893        for (int i = 0; i < length; i++) {
894            result[i] = array1[i].subtract(array2[i]);
895        }
896
897        return result;
898    }
899
900    /** Return a new String representing the array, formatted as
901     *  in Java array initializers.
902     *
903     *  @param array An array of complex numbers.
904     *  @return A new string representing the array.
905     */
906    public static final String toString(Complex[] array) {
907        return toString(array, ", ", "{", "}");
908    }
909
910    /** Return a new String representing the array, formatted
911     *  specified starting with vectorBegin, where each
912     *  successive element is separated by elementDelimiter
913     *  and ending with vectorEnd.
914     *
915     *  @param array An array of complex numbers.
916     *  @param elementDelimiter The delimiter between elements,
917     *  typically ", ".
918     *  @param vectorBegin The start of the array, typically "{".
919     *  @param vectorEnd The end of the array, typically "}".
920     *  @return A new string representing the array in the format specified
921     *  by the ArrayStringFormat argument.
922     */
923    public static final String toString(final Complex[] array,
924            String elementDelimiter, String vectorBegin, String vectorEnd) {
925        int length = array.length;
926        StringBuffer sb = new StringBuffer();
927
928        sb.append(vectorBegin);
929
930        for (int i = 0; i < length; i++) {
931            sb.append(array[i].toString());
932
933            if (i < length - 1) {
934                sb.append(elementDelimiter);
935            }
936        }
937
938        sb.append(vectorEnd);
939
940        return sb.toString();
941    }
942
943    /** Return true if all the distances between corresponding elements
944     *  <i>array1</i> and <i>array2</i> are all less than or equal to
945     *  the magnitude of <i>maxError</i>. If both arrays are empty,
946     *  return true.
947     *
948     *  @param array1 The first array.
949     *  @param array2 The second array.
950     *  @param maxError A complex number whose magnitude is taken to
951     *   be the distance threshold.
952     *  @return True if all the distances between corresponding elements
953     *  <i>array1</i> and <i>array2</i> are all less than or equal to
954     *  the magnitude of <i>maxError</i>.
955     *  @exception IllegalArgumentException If the arrays are not of the same
956     *   length.
957     */
958    public static final boolean within(Complex[] array1, Complex[] array2,
959            Complex maxError) {
960        return within(array1, array2, maxError.magnitude());
961    }
962
963    /** Return true if all the distances between corresponding elements
964     *  <i>array1</i> and <i>array2</i> are all less than or equal to
965     *  <i>maxError</i>. If both arrays are empty, return true.
966     *  If <i>maxError</i> is negative, return false.
967     *
968     *  @param array1 The first array.
969     *  @param array2 The second array.
970     *  @param maxError The threshold for the magnitude of the difference.
971     *  @return True if all the distances between corresponding elements
972     *  <i>array1</i> and <i>array2</i> are all less than or equal to
973     *  <i>maxError</i>.
974     *  @exception IllegalArgumentException If the arrays are not of the same
975     *   length.
976     */
977    public static final boolean within(Complex[] array1, Complex[] array2,
978            double maxError) {
979        int length = _commonLength(array1, array2, "ComplexArrayMath.within");
980
981        for (int i = 0; i < length; i++) {
982            if (!array1[i].isCloseTo(array2[i], maxError)) {
983                return false;
984            }
985        }
986
987        return true;
988    }
989
990    /** Return true if all the distances between corresponding elements
991     *  <i>array1</i> and <i>array2</i> are all less than or equal to
992     *  the corresponding elements in <i>maxError</i>. If both arrays
993     *  are empty, return true. If any element in <i>maxError</i> is negative,
994     *  return false.
995     *
996     *  @param array1 The first array.
997     *  @param array2 The second array.
998     *  @param maxError The array of thresholds for the magnitudes of
999     *   the difference.
1000     *  @return True if all the distances between corresponding elements
1001     *  <i>array1</i> and <i>array2</i> are all less than or equal to
1002     *  the corresponding elements in <i>maxError</i>.
1003     *  @exception IllegalArgumentException If the arrays are not of the same
1004     *   length.
1005     */
1006    public static final boolean within(Complex[] array1, Complex[] array2,
1007            double[] maxError) {
1008        int length = _commonLength(array1, array2, "ComplexArrayMath.within");
1009
1010        for (int i = 0; i < length; i++) {
1011            if (!array1[i].isCloseTo(array2[i], maxError[i])) {
1012                return false;
1013            }
1014        }
1015
1016        return true;
1017    }
1018
1019    /** Return true if all the distances between corresponding elements
1020     *  <i>array1</i> and <i>array2</i> are all less than or equal to
1021     *  the magnitudes of the corresponding elements in <i>maxError</i>.
1022     *  If both arrays are empty, return true.
1023     *
1024     *  <p>Note that there is no notion of negative distance with
1025     *  complex numbers, so unlike the within() methods for other
1026     *  types, this method will not return false if an element of the
1027     *  maxError matrix is negative.
1028     *
1029     *  @param array1 The first array.
1030     *  @param array2 The second array.
1031     *  @param maxError An array of complex numbers whose magnitude
1032     *  for each element is taken to be the distance threshold.
1033     *  @return True if all the distances between corresponding elements
1034     *  <i>array1</i> and <i>array2</i> are all less than or equal to
1035     *  the magnitudes of the corresponding elements in <i>maxError</i>.
1036     *  @exception IllegalArgumentException If the arrays are not of the same
1037     *   length.
1038     */
1039    public static final boolean within(Complex[] array1, Complex[] array2,
1040            Complex[] maxError) {
1041        int length = maxError.length;
1042        double[] doubleError = new double[length];
1043
1044        for (int i = 0; i < length; i++) {
1045            doubleError[i] = maxError[i].magnitude();
1046        }
1047
1048        return within(array1, array2, doubleError);
1049    }
1050
1051    ///////////////////////////////////////////////////////////////////
1052    //    protected methods
1053
1054    /** Throw an exception if the two arrays are not of the same length,
1055     *  or if either array is null. An exception is NOT thrown if both
1056     *  arrays are of length 0. If no exception is thrown, return the common
1057     *  length of the arrays.
1058     *  @param array1 The first array of doubles.
1059     *  @param array2 The second array of doubles.
1060     *  @param methodName A String representing the method name of the caller,
1061     *  without parentheses.
1062     *  @return The common length of both arrays.
1063     */
1064    protected static final int _commonLength(final Complex[] array1,
1065            final Complex[] array2, final String methodName) {
1066        if (array1 == null) {
1067            throw new IllegalArgumentException("ptolemy.math." + methodName
1068                    + "() : first input array is null.");
1069        }
1070
1071        if (array2 == null) {
1072            throw new IllegalArgumentException("ptolemy.math." + methodName
1073                    + "() : second input array is null.");
1074        }
1075
1076        if (array1.length != array2.length) {
1077            throw new IllegalArgumentException("ptolemy.math." + methodName
1078                    + "() : input arrays must have the same length, "
1079                    + "but the first array has length " + array1.length
1080                    + " and the second array has length " + array2.length
1081                    + '.');
1082        }
1083
1084        return array1.length;
1085    }
1086}