001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.analysis.differentiation;
018
019import java.io.Serializable;
020
021import org.apache.commons.numbers.arrays.LinearCombination;
022import org.apache.commons.math4.Field;
023import org.apache.commons.math4.FieldElement;
024import org.apache.commons.math4.RealFieldElement;
025import org.apache.commons.math4.exception.DimensionMismatchException;
026import org.apache.commons.math4.exception.MathArithmeticException;
027import org.apache.commons.math4.exception.NumberIsTooLargeException;
028import org.apache.commons.math4.util.FastMath;
029import org.apache.commons.math4.util.MathArrays;
030import org.apache.commons.math4.util.MathUtils;
031
032/** Class representing both the value and the differentials of a function.
033 * <p>This class is the workhorse of the differentiation package.</p>
034 * <p>This class is an implementation of the extension to Rall's
035 * numbers described in Dan Kalman's paper <a
036 * href="http://www1.american.edu/cas/mathstat/People/kalman/pdffiles/mmgautodiff.pdf">Doubly
037 * Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75,
038 * no. 3, June 2002. Rall's numbers are an extension to the real numbers used
039 * throughout mathematical expressions; they hold the derivative together with the
040 * value of a function. Dan Kalman's derivative structures hold all partial derivatives
041 * up to any specified order, with respect to any number of free parameters. Rall's
042 * numbers therefore can be seen as derivative structures for order one derivative and
043 * one free parameter, and real numbers can be seen as derivative structures with zero
044 * order derivative and no free parameters.</p>
045 * <p>{@link DerivativeStructure} instances can be used directly thanks to
046 * the arithmetic operators to the mathematical functions provided as
047 * methods by this class (+, -, *, /, %, sin, cos ...).</p>
048 * <p>Implementing complex expressions by hand using these classes is
049 * a tedious and error-prone task but has the advantage of having no limitation
050 * on the derivation order despite not requiring users to compute the derivatives by
051 * themselves. Implementing complex expression can also be done by developing computation
052 * code using standard primitive double values and to use {@link
053 * UnivariateFunctionDifferentiator differentiators} to create the {@link
054 * DerivativeStructure}-based instances. This method is simpler but may be limited in
055 * the accuracy and derivation orders and may be computationally intensive (this is
056 * typically the case for {@link FiniteDifferencesDifferentiator finite differences
057 * differentiator}.</p>
058 * <p>Instances of this class are guaranteed to be immutable.</p>
059 * @see DSCompiler
060 * @since 3.1
061 */
062public class DerivativeStructure implements RealFieldElement<DerivativeStructure>, Serializable {
063
064    /** Serializable UID. */
065    private static final long serialVersionUID = 20120730L;
066
067    /** Compiler for the current dimensions. */
068    private transient DSCompiler compiler;
069
070    /** Combined array holding all values. */
071    private final double[] data;
072
073    /** Build an instance with all values and derivatives set to 0.
074     * @param compiler compiler to use for computation
075     */
076    private DerivativeStructure(final DSCompiler compiler) {
077        this.compiler = compiler;
078        this.data     = new double[compiler.getSize()];
079    }
080
081    /** Build an instance with all values and derivatives set to 0.
082     * @param parameters number of free parameters
083     * @param order derivation order
084     * @throws NumberIsTooLargeException if order is too large
085     */
086    public DerivativeStructure(final int parameters, final int order)
087        throws NumberIsTooLargeException {
088        this(DSCompiler.getCompiler(parameters, order));
089    }
090
091    /** Build an instance representing a constant value.
092     * @param parameters number of free parameters
093     * @param order derivation order
094     * @param value value of the constant
095     * @throws NumberIsTooLargeException if order is too large
096     * @see #DerivativeStructure(int, int, int, double)
097     */
098    public DerivativeStructure(final int parameters, final int order, final double value)
099        throws NumberIsTooLargeException {
100        this(parameters, order);
101        this.data[0] = value;
102    }
103
104    /** Build an instance representing a variable.
105     * <p>Instances built using this constructor are considered
106     * to be the free variables with respect to which differentials
107     * are computed. As such, their differential with respect to
108     * themselves is +1.</p>
109     * @param parameters number of free parameters
110     * @param order derivation order
111     * @param index index of the variable (from 0 to {@code parameters - 1})
112     * @param value value of the variable
113     * @exception NumberIsTooLargeException if {@code index &ge; parameters}.
114     * @see #DerivativeStructure(int, int, double)
115     */
116    public DerivativeStructure(final int parameters, final int order,
117                               final int index, final double value)
118        throws NumberIsTooLargeException {
119        this(parameters, order, value);
120
121        if (index >= parameters) {
122            throw new NumberIsTooLargeException(index, parameters, false);
123        }
124
125        if (order > 0) {
126            // the derivative of the variable with respect to itself is 1.
127            data[DSCompiler.getCompiler(index, order).getSize()] = 1.0;
128        }
129
130    }
131
132    /** Linear combination constructor.
133     * The derivative structure built will be a1 * ds1 + a2 * ds2
134     * @param a1 first scale factor
135     * @param ds1 first base (unscaled) derivative structure
136     * @param a2 second scale factor
137     * @param ds2 second base (unscaled) derivative structure
138     * @exception DimensionMismatchException if number of free parameters or orders are inconsistent
139     */
140    public DerivativeStructure(final double a1, final DerivativeStructure ds1,
141                               final double a2, final DerivativeStructure ds2)
142        throws DimensionMismatchException {
143        this(ds1.compiler);
144        compiler.checkCompatibility(ds2.compiler);
145        compiler.linearCombination(a1, ds1.data, 0, a2, ds2.data, 0, data, 0);
146    }
147
148    /** Linear combination constructor.
149     * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3
150     * @param a1 first scale factor
151     * @param ds1 first base (unscaled) derivative structure
152     * @param a2 second scale factor
153     * @param ds2 second base (unscaled) derivative structure
154     * @param a3 third scale factor
155     * @param ds3 third base (unscaled) derivative structure
156     * @exception DimensionMismatchException if number of free parameters or orders are inconsistent
157     */
158    public DerivativeStructure(final double a1, final DerivativeStructure ds1,
159                               final double a2, final DerivativeStructure ds2,
160                               final double a3, final DerivativeStructure ds3)
161        throws DimensionMismatchException {
162        this(ds1.compiler);
163        compiler.checkCompatibility(ds2.compiler);
164        compiler.checkCompatibility(ds3.compiler);
165        compiler.linearCombination(a1, ds1.data, 0, a2, ds2.data, 0, a3, ds3.data, 0, data, 0);
166    }
167
168    /** Linear combination constructor.
169     * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
170     * @param a1 first scale factor
171     * @param ds1 first base (unscaled) derivative structure
172     * @param a2 second scale factor
173     * @param ds2 second base (unscaled) derivative structure
174     * @param a3 third scale factor
175     * @param ds3 third base (unscaled) derivative structure
176     * @param a4 fourth scale factor
177     * @param ds4 fourth base (unscaled) derivative structure
178     * @exception DimensionMismatchException if number of free parameters or orders are inconsistent
179     */
180    public DerivativeStructure(final double a1, final DerivativeStructure ds1,
181                               final double a2, final DerivativeStructure ds2,
182                               final double a3, final DerivativeStructure ds3,
183                               final double a4, final DerivativeStructure ds4)
184        throws DimensionMismatchException {
185        this(ds1.compiler);
186        compiler.checkCompatibility(ds2.compiler);
187        compiler.checkCompatibility(ds3.compiler);
188        compiler.checkCompatibility(ds4.compiler);
189        compiler.linearCombination(a1, ds1.data, 0, a2, ds2.data, 0,
190                                   a3, ds3.data, 0, a4, ds4.data, 0,
191                                   data, 0);
192    }
193
194    /** Build an instance from all its derivatives.
195     * @param parameters number of free parameters
196     * @param order derivation order
197     * @param derivatives derivatives sorted according to
198     * {@link DSCompiler#getPartialDerivativeIndex(int...)}
199     * @exception DimensionMismatchException if derivatives array does not match the
200     * {@link DSCompiler#getSize() size} expected by the compiler
201     * @throws NumberIsTooLargeException if order is too large
202     * @see #getAllDerivatives()
203     */
204    public DerivativeStructure(final int parameters, final int order, final double ... derivatives)
205        throws DimensionMismatchException, NumberIsTooLargeException {
206        this(parameters, order);
207        if (derivatives.length != data.length) {
208            throw new DimensionMismatchException(derivatives.length, data.length);
209        }
210        System.arraycopy(derivatives, 0, data, 0, data.length);
211    }
212
213    /** Copy constructor.
214     * @param ds instance to copy
215     */
216    private DerivativeStructure(final DerivativeStructure ds) {
217        this.compiler = ds.compiler;
218        this.data     = ds.data.clone();
219    }
220
221    /** Get the number of free parameters.
222     * @return number of free parameters
223     */
224    public int getFreeParameters() {
225        return compiler.getFreeParameters();
226    }
227
228    /** Get the derivation order.
229     * @return derivation order
230     */
231    public int getOrder() {
232        return compiler.getOrder();
233    }
234
235    /** Create a constant compatible with instance order and number of parameters.
236     * <p>
237     * This method is a convenience factory method, it simply calls
238     * {@code new DerivativeStructure(getFreeParameters(), getOrder(), c)}
239     * </p>
240     * @param c value of the constant
241     * @return a constant compatible with instance order and number of parameters
242     * @see #DerivativeStructure(int, int, double)
243     * @since 3.3
244     */
245    public DerivativeStructure createConstant(final double c) {
246        return new DerivativeStructure(getFreeParameters(), getOrder(), c);
247    }
248
249    /** {@inheritDoc}
250     * @since 3.2
251     */
252    @Override
253    public double getReal() {
254        return data[0];
255    }
256
257    /** Get the value part of the derivative structure.
258     * @return value part of the derivative structure
259     * @see #getPartialDerivative(int...)
260     */
261    public double getValue() {
262        return data[0];
263    }
264
265    /** Get a partial derivative.
266     * @param orders derivation orders with respect to each variable (if all orders are 0,
267     * the value is returned)
268     * @return partial derivative
269     * @see #getValue()
270     * @exception DimensionMismatchException if the numbers of variables does not
271     * match the instance
272     * @exception NumberIsTooLargeException if sum of derivation orders is larger
273     * than the instance limits
274     */
275    public double getPartialDerivative(final int ... orders)
276        throws DimensionMismatchException, NumberIsTooLargeException {
277        return data[compiler.getPartialDerivativeIndex(orders)];
278    }
279
280    /** Get all partial derivatives.
281     * @return a fresh copy of partial derivatives, in an array sorted according to
282     * {@link DSCompiler#getPartialDerivativeIndex(int...)}
283     */
284    public double[] getAllDerivatives() {
285        return data.clone();
286    }
287
288    /** {@inheritDoc}
289     * @since 3.2
290     */
291    @Override
292    public DerivativeStructure add(final double a) {
293        final DerivativeStructure ds = new DerivativeStructure(this);
294        ds.data[0] += a;
295        return ds;
296    }
297
298    /** {@inheritDoc}
299     * @exception DimensionMismatchException if number of free parameters
300     * or orders do not match
301     */
302    @Override
303    public DerivativeStructure add(final DerivativeStructure a)
304        throws DimensionMismatchException {
305        compiler.checkCompatibility(a.compiler);
306        final DerivativeStructure ds = new DerivativeStructure(this);
307        compiler.add(data, 0, a.data, 0, ds.data, 0);
308        return ds;
309    }
310
311    /** {@inheritDoc}
312     * @since 3.2
313     */
314    @Override
315    public DerivativeStructure subtract(final double a) {
316        return add(-a);
317    }
318
319    /** {@inheritDoc}
320     * @exception DimensionMismatchException if number of free parameters
321     * or orders do not match
322     */
323    @Override
324    public DerivativeStructure subtract(final DerivativeStructure a)
325        throws DimensionMismatchException {
326        compiler.checkCompatibility(a.compiler);
327        final DerivativeStructure ds = new DerivativeStructure(this);
328        compiler.subtract(data, 0, a.data, 0, ds.data, 0);
329        return ds;
330    }
331
332    /** {@inheritDoc} */
333    @Override
334    public DerivativeStructure multiply(final int n) {
335        return multiply((double) n);
336    }
337
338    /** {@inheritDoc}
339     * @since 3.2
340     */
341    @Override
342    public DerivativeStructure multiply(final double a) {
343        final DerivativeStructure ds = new DerivativeStructure(this);
344        for (int i = 0; i < ds.data.length; ++i) {
345            ds.data[i] *= a;
346        }
347        return ds;
348    }
349
350    /** {@inheritDoc}
351     * @exception DimensionMismatchException if number of free parameters
352     * or orders do not match
353     */
354    @Override
355    public DerivativeStructure multiply(final DerivativeStructure a)
356        throws DimensionMismatchException {
357        compiler.checkCompatibility(a.compiler);
358        final DerivativeStructure result = new DerivativeStructure(compiler);
359        compiler.multiply(data, 0, a.data, 0, result.data, 0);
360        return result;
361    }
362
363    /** {@inheritDoc}
364     * @since 3.2
365     */
366    @Override
367    public DerivativeStructure divide(final double a) {
368        final DerivativeStructure ds = new DerivativeStructure(this);
369        for (int i = 0; i < ds.data.length; ++i) {
370            ds.data[i] /= a;
371        }
372        return ds;
373    }
374
375    /** {@inheritDoc}
376     * @exception DimensionMismatchException if number of free parameters
377     * or orders do not match
378     */
379    @Override
380    public DerivativeStructure divide(final DerivativeStructure a)
381        throws DimensionMismatchException {
382        compiler.checkCompatibility(a.compiler);
383        final DerivativeStructure result = new DerivativeStructure(compiler);
384        compiler.divide(data, 0, a.data, 0, result.data, 0);
385        return result;
386    }
387
388    /** {@inheritDoc} */
389    @Override
390    public DerivativeStructure remainder(final double a) {
391        final DerivativeStructure ds = new DerivativeStructure(this);
392        ds.data[0] = FastMath.IEEEremainder(ds.data[0], a);
393        return ds;
394    }
395
396    /** {@inheritDoc}
397     * @exception DimensionMismatchException if number of free parameters
398     * or orders do not match
399     * @since 3.2
400     */
401    @Override
402    public DerivativeStructure remainder(final DerivativeStructure a)
403        throws DimensionMismatchException {
404        compiler.checkCompatibility(a.compiler);
405        final DerivativeStructure result = new DerivativeStructure(compiler);
406        compiler.remainder(data, 0, a.data, 0, result.data, 0);
407        return result;
408    }
409
410    /** {@inheritDoc} */
411    @Override
412    public DerivativeStructure negate() {
413        final DerivativeStructure ds = new DerivativeStructure(compiler);
414        for (int i = 0; i < ds.data.length; ++i) {
415            ds.data[i] = -data[i];
416        }
417        return ds;
418    }
419
420    /** {@inheritDoc}
421     * @since 3.2
422     */
423    @Override
424    public DerivativeStructure abs() {
425        if (Double.doubleToLongBits(data[0]) < 0) {
426            // we use the bits representation to also handle -0.0
427            return negate();
428        } else {
429            return this;
430        }
431    }
432
433    /** {@inheritDoc}
434     * @since 3.2
435     */
436    @Override
437    public DerivativeStructure ceil() {
438        return new DerivativeStructure(compiler.getFreeParameters(),
439                                       compiler.getOrder(),
440                                       FastMath.ceil(data[0]));
441    }
442
443    /** {@inheritDoc}
444     * @since 3.2
445     */
446    @Override
447    public DerivativeStructure floor() {
448        return new DerivativeStructure(compiler.getFreeParameters(),
449                                       compiler.getOrder(),
450                                       FastMath.floor(data[0]));
451    }
452
453    /** {@inheritDoc}
454     * @since 3.2
455     */
456    @Override
457    public DerivativeStructure rint() {
458        return new DerivativeStructure(compiler.getFreeParameters(),
459                                       compiler.getOrder(),
460                                       FastMath.rint(data[0]));
461    }
462
463    /** {@inheritDoc} */
464    @Override
465    public long round() {
466        return FastMath.round(data[0]);
467    }
468
469    /** {@inheritDoc}
470     * @since 3.2
471     */
472    @Override
473    public DerivativeStructure signum() {
474        return new DerivativeStructure(compiler.getFreeParameters(),
475                                       compiler.getOrder(),
476                                       FastMath.signum(data[0]));
477    }
478
479    /** {@inheritDoc}
480     * @since 3.2
481     */
482    @Override
483    public DerivativeStructure copySign(final DerivativeStructure sign){
484        long m = Double.doubleToLongBits(data[0]);
485        long s = Double.doubleToLongBits(sign.data[0]);
486        if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
487            return this;
488        }
489        return negate(); // flip sign
490    }
491
492    /** {@inheritDoc}
493     * @since 3.2
494     */
495    @Override
496    public DerivativeStructure copySign(final double sign) {
497        long m = Double.doubleToLongBits(data[0]);
498        long s = Double.doubleToLongBits(sign);
499        if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
500            return this;
501        }
502        return negate(); // flip sign
503    }
504
505    /**
506     * Return the exponent of the instance value, removing the bias.
507     * <p>
508     * For double numbers of the form 2<sup>x</sup>, the unbiased
509     * exponent is exactly x.
510     * </p>
511     * @return exponent for instance in IEEE754 representation, without bias
512     */
513    public int getExponent() {
514        return FastMath.getExponent(data[0]);
515    }
516
517    /** {@inheritDoc}
518     * @since 3.2
519     */
520    @Override
521    public DerivativeStructure scalb(final int n) {
522        final DerivativeStructure ds = new DerivativeStructure(compiler);
523        for (int i = 0; i < ds.data.length; ++i) {
524            ds.data[i] = FastMath.scalb(data[i], n);
525        }
526        return ds;
527    }
528
529    /** {@inheritDoc}
530     * @exception DimensionMismatchException if number of free parameters
531     * or orders do not match
532     * @since 3.2
533     */
534    @Override
535    public DerivativeStructure hypot(final DerivativeStructure y)
536        throws DimensionMismatchException {
537
538        compiler.checkCompatibility(y.compiler);
539
540        if (Double.isInfinite(data[0]) || Double.isInfinite(y.data[0])) {
541            return new DerivativeStructure(compiler.getFreeParameters(),
542                                           compiler.getFreeParameters(),
543                                           Double.POSITIVE_INFINITY);
544        } else if (Double.isNaN(data[0]) || Double.isNaN(y.data[0])) {
545            return new DerivativeStructure(compiler.getFreeParameters(),
546                                           compiler.getFreeParameters(),
547                                           Double.NaN);
548        } else {
549
550            final int expX = getExponent();
551            final int expY = y.getExponent();
552            if (expX > expY + 27) {
553                // y is neglectible with respect to x
554                return abs();
555            } else if (expY > expX + 27) {
556                // x is neglectible with respect to y
557                return y.abs();
558            } else {
559
560                // find an intermediate scale to avoid both overflow and underflow
561                final int middleExp = (expX + expY) / 2;
562
563                // scale parameters without losing precision
564                final DerivativeStructure scaledX = scalb(-middleExp);
565                final DerivativeStructure scaledY = y.scalb(-middleExp);
566
567                // compute scaled hypotenuse
568                final DerivativeStructure scaledH =
569                        scaledX.multiply(scaledX).add(scaledY.multiply(scaledY)).sqrt();
570
571                // remove scaling
572                return scaledH.scalb(middleExp);
573
574            }
575
576        }
577    }
578
579    /**
580     * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
581     * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
582     * avoiding intermediate overflow or underflow.
583     *
584     * <ul>
585     * <li> If either argument is infinite, then the result is positive infinity.</li>
586     * <li> else, if either argument is NaN then the result is NaN.</li>
587     * </ul>
588     *
589     * @param x a value
590     * @param y a value
591     * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
592     * @exception DimensionMismatchException if number of free parameters
593     * or orders do not match
594     * @since 3.2
595     */
596    public static DerivativeStructure hypot(final DerivativeStructure x, final DerivativeStructure y)
597        throws DimensionMismatchException {
598        return x.hypot(y);
599    }
600
601    /** Compute composition of the instance by a univariate function.
602     * @param f array of value and derivatives of the function at
603     * the current point (i.e. [f({@link #getValue()}),
604     * f'({@link #getValue()}), f''({@link #getValue()})...]).
605     * @return f(this)
606     * @exception DimensionMismatchException if the number of derivatives
607     * in the array is not equal to {@link #getOrder() order} + 1
608     */
609    public DerivativeStructure compose(final double ... f)
610        throws DimensionMismatchException {
611        if (f.length != getOrder() + 1) {
612            throw new DimensionMismatchException(f.length, getOrder() + 1);
613        }
614        final DerivativeStructure result = new DerivativeStructure(compiler);
615        compiler.compose(data, 0, f, result.data, 0);
616        return result;
617    }
618
619    /** {@inheritDoc} */
620    @Override
621    public DerivativeStructure reciprocal() {
622        final DerivativeStructure result = new DerivativeStructure(compiler);
623        compiler.pow(data, 0, -1, result.data, 0);
624        return result;
625    }
626
627    /** {@inheritDoc}
628     * @since 3.2
629     */
630    @Override
631    public DerivativeStructure sqrt() {
632        return rootN(2);
633    }
634
635    /** {@inheritDoc}
636     * @since 3.2
637     */
638    @Override
639    public DerivativeStructure cbrt() {
640        return rootN(3);
641    }
642
643    /** {@inheritDoc}
644     * @since 3.2
645     */
646    @Override
647    public DerivativeStructure rootN(final int n) {
648        final DerivativeStructure result = new DerivativeStructure(compiler);
649        compiler.rootN(data, 0, n, result.data, 0);
650        return result;
651    }
652
653    /** {@inheritDoc} */
654    @Override
655    public Field<DerivativeStructure> getField() {
656        return new Field<DerivativeStructure>() {
657
658            /** {@inheritDoc} */
659            @Override
660            public DerivativeStructure getZero() {
661                return new DerivativeStructure(compiler.getFreeParameters(), compiler.getOrder(), 0.0);
662            }
663
664            /** {@inheritDoc} */
665            @Override
666            public DerivativeStructure getOne() {
667                return new DerivativeStructure(compiler.getFreeParameters(), compiler.getOrder(), 1.0);
668            }
669
670            /** {@inheritDoc} */
671            @Override
672            public Class<? extends FieldElement<DerivativeStructure>> getRuntimeClass() {
673                return DerivativeStructure.class;
674            }
675
676        };
677    }
678
679    /** Compute a<sup>x</sup> where a is a double and x a {@link DerivativeStructure}
680     * @param a number to exponentiate
681     * @param x power to apply
682     * @return a<sup>x</sup>
683     * @since 3.3
684     */
685    public static DerivativeStructure pow(final double a, final DerivativeStructure x) {
686        final DerivativeStructure result = new DerivativeStructure(x.compiler);
687        x.compiler.pow(a, x.data, 0, result.data, 0);
688        return result;
689    }
690
691    /** {@inheritDoc}
692     * @since 3.2
693     */
694    @Override
695    public DerivativeStructure pow(final double p) {
696        final DerivativeStructure result = new DerivativeStructure(compiler);
697        compiler.pow(data, 0, p, result.data, 0);
698        return result;
699    }
700
701    /** {@inheritDoc}
702     * @since 3.2
703     */
704    @Override
705    public DerivativeStructure pow(final int n) {
706        final DerivativeStructure result = new DerivativeStructure(compiler);
707        compiler.pow(data, 0, n, result.data, 0);
708        return result;
709    }
710
711    /** {@inheritDoc}
712     * @exception DimensionMismatchException if number of free parameters
713     * or orders do not match
714     * @since 3.2
715     */
716    @Override
717    public DerivativeStructure pow(final DerivativeStructure e)
718        throws DimensionMismatchException {
719        compiler.checkCompatibility(e.compiler);
720        final DerivativeStructure result = new DerivativeStructure(compiler);
721        compiler.pow(data, 0, e.data, 0, result.data, 0);
722        return result;
723    }
724
725    /** {@inheritDoc}
726     * @since 3.2
727     */
728    @Override
729    public DerivativeStructure exp() {
730        final DerivativeStructure result = new DerivativeStructure(compiler);
731        compiler.exp(data, 0, result.data, 0);
732        return result;
733    }
734
735    /** {@inheritDoc}
736     * @since 3.2
737     */
738    @Override
739    public DerivativeStructure expm1() {
740        final DerivativeStructure result = new DerivativeStructure(compiler);
741        compiler.expm1(data, 0, result.data, 0);
742        return result;
743    }
744
745    /** {@inheritDoc}
746     * @since 3.2
747     */
748    @Override
749    public DerivativeStructure log() {
750        final DerivativeStructure result = new DerivativeStructure(compiler);
751        compiler.log(data, 0, result.data, 0);
752        return result;
753    }
754
755    /** {@inheritDoc}
756     * @since 3.2
757     */
758    @Override
759    public DerivativeStructure log1p() {
760        final DerivativeStructure result = new DerivativeStructure(compiler);
761        compiler.log1p(data, 0, result.data, 0);
762        return result;
763    }
764
765    /** Base 10 logarithm.
766     * @return base 10 logarithm of the instance
767     */
768    @Override
769    public DerivativeStructure log10() {
770        final DerivativeStructure result = new DerivativeStructure(compiler);
771        compiler.log10(data, 0, result.data, 0);
772        return result;
773    }
774
775    /** {@inheritDoc}
776     * @since 3.2
777     */
778    @Override
779    public DerivativeStructure cos() {
780        final DerivativeStructure result = new DerivativeStructure(compiler);
781        compiler.cos(data, 0, result.data, 0);
782        return result;
783    }
784
785    /** {@inheritDoc}
786     * @since 3.2
787     */
788    @Override
789    public DerivativeStructure sin() {
790        final DerivativeStructure result = new DerivativeStructure(compiler);
791        compiler.sin(data, 0, result.data, 0);
792        return result;
793    }
794
795    /** {@inheritDoc}
796     * @since 3.2
797     */
798    @Override
799    public DerivativeStructure tan() {
800        final DerivativeStructure result = new DerivativeStructure(compiler);
801        compiler.tan(data, 0, result.data, 0);
802        return result;
803    }
804
805    /** {@inheritDoc}
806     * @since 3.2
807     */
808    @Override
809    public DerivativeStructure acos() {
810        final DerivativeStructure result = new DerivativeStructure(compiler);
811        compiler.acos(data, 0, result.data, 0);
812        return result;
813    }
814
815    /** {@inheritDoc}
816     * @since 3.2
817     */
818    @Override
819    public DerivativeStructure asin() {
820        final DerivativeStructure result = new DerivativeStructure(compiler);
821        compiler.asin(data, 0, result.data, 0);
822        return result;
823    }
824
825    /** {@inheritDoc}
826     * @since 3.2
827     */
828    @Override
829    public DerivativeStructure atan() {
830        final DerivativeStructure result = new DerivativeStructure(compiler);
831        compiler.atan(data, 0, result.data, 0);
832        return result;
833    }
834
835    /** {@inheritDoc}
836     * @since 3.2
837     */
838    @Override
839    public DerivativeStructure atan2(final DerivativeStructure x)
840        throws DimensionMismatchException {
841        compiler.checkCompatibility(x.compiler);
842        final DerivativeStructure result = new DerivativeStructure(compiler);
843        compiler.atan2(data, 0, x.data, 0, result.data, 0);
844        return result;
845    }
846
847    /** Two arguments arc tangent operation.
848     * @param y first argument of the arc tangent
849     * @param x second argument of the arc tangent
850     * @return atan2(y, x)
851     * @exception DimensionMismatchException if number of free parameters
852     * or orders do not match
853     * @since 3.2
854     */
855    public static DerivativeStructure atan2(final DerivativeStructure y, final DerivativeStructure x)
856        throws DimensionMismatchException {
857        return y.atan2(x);
858    }
859
860    /** {@inheritDoc}
861     * @since 3.2
862     */
863    @Override
864    public DerivativeStructure cosh() {
865        final DerivativeStructure result = new DerivativeStructure(compiler);
866        compiler.cosh(data, 0, result.data, 0);
867        return result;
868    }
869
870    /** {@inheritDoc}
871     * @since 3.2
872     */
873    @Override
874    public DerivativeStructure sinh() {
875        final DerivativeStructure result = new DerivativeStructure(compiler);
876        compiler.sinh(data, 0, result.data, 0);
877        return result;
878    }
879
880    /** {@inheritDoc}
881     * @since 3.2
882     */
883    @Override
884    public DerivativeStructure tanh() {
885        final DerivativeStructure result = new DerivativeStructure(compiler);
886        compiler.tanh(data, 0, result.data, 0);
887        return result;
888    }
889
890    /** {@inheritDoc}
891     * @since 3.2
892     */
893    @Override
894    public DerivativeStructure acosh() {
895        final DerivativeStructure result = new DerivativeStructure(compiler);
896        compiler.acosh(data, 0, result.data, 0);
897        return result;
898    }
899
900    /** {@inheritDoc}
901     * @since 3.2
902     */
903    @Override
904    public DerivativeStructure asinh() {
905        final DerivativeStructure result = new DerivativeStructure(compiler);
906        compiler.asinh(data, 0, result.data, 0);
907        return result;
908    }
909
910    /** {@inheritDoc}
911     * @since 3.2
912     */
913    @Override
914    public DerivativeStructure atanh() {
915        final DerivativeStructure result = new DerivativeStructure(compiler);
916        compiler.atanh(data, 0, result.data, 0);
917        return result;
918    }
919
920    /** Convert radians to degrees, with error of less than 0.5 ULP
921     *  @return instance converted into degrees
922     */
923    public DerivativeStructure toDegrees() {
924        final DerivativeStructure ds = new DerivativeStructure(compiler);
925        for (int i = 0; i < ds.data.length; ++i) {
926            ds.data[i] = FastMath.toDegrees(data[i]);
927        }
928        return ds;
929    }
930
931    /** Convert degrees to radians, with error of less than 0.5 ULP
932     *  @return instance converted into radians
933     */
934    public DerivativeStructure toRadians() {
935        final DerivativeStructure ds = new DerivativeStructure(compiler);
936        for (int i = 0; i < ds.data.length; ++i) {
937            ds.data[i] = FastMath.toRadians(data[i]);
938        }
939        return ds;
940    }
941
942    /** Evaluate Taylor expansion a derivative structure.
943     * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
944     * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
945     * @throws MathArithmeticException if factorials becomes too large
946     */
947    public double taylor(final double ... delta) throws MathArithmeticException {
948        return compiler.taylor(data, 0, delta);
949    }
950
951    /** {@inheritDoc}
952     * @exception DimensionMismatchException if number of free parameters
953     * or orders do not match
954     * @since 3.2
955     */
956    @Override
957    public DerivativeStructure linearCombination(final DerivativeStructure[] a, final DerivativeStructure[] b)
958        throws DimensionMismatchException {
959
960        // compute an accurate value, taking care of cancellations
961        final double[] aDouble = new double[a.length];
962        for (int i = 0; i < a.length; ++i) {
963            aDouble[i] = a[i].getValue();
964        }
965        final double[] bDouble = new double[b.length];
966        for (int i = 0; i < b.length; ++i) {
967            bDouble[i] = b[i].getValue();
968        }
969        final double accurateValue = LinearCombination.value(aDouble, bDouble);
970
971        // compute a simple value, with all partial derivatives
972        DerivativeStructure simpleValue = a[0].getField().getZero();
973        for (int i = 0; i < a.length; ++i) {
974            simpleValue = simpleValue.add(a[i].multiply(b[i]));
975        }
976
977        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
978        final double[] all = simpleValue.getAllDerivatives();
979        all[0] = accurateValue;
980        return new DerivativeStructure(simpleValue.getFreeParameters(), simpleValue.getOrder(), all);
981
982    }
983
984    /** {@inheritDoc}
985     * @exception DimensionMismatchException if number of free parameters
986     * or orders do not match
987     * @since 3.2
988     */
989    @Override
990    public DerivativeStructure linearCombination(final double[] a, final DerivativeStructure[] b)
991        throws DimensionMismatchException {
992
993        // compute an accurate value, taking care of cancellations
994        final double[] bDouble = new double[b.length];
995        for (int i = 0; i < b.length; ++i) {
996            bDouble[i] = b[i].getValue();
997        }
998        final double accurateValue = LinearCombination.value(a, bDouble);
999
1000        // compute a simple value, with all partial derivatives
1001        DerivativeStructure simpleValue = b[0].getField().getZero();
1002        for (int i = 0; i < a.length; ++i) {
1003            simpleValue = simpleValue.add(b[i].multiply(a[i]));
1004        }
1005
1006        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1007        final double[] all = simpleValue.getAllDerivatives();
1008        all[0] = accurateValue;
1009        return new DerivativeStructure(simpleValue.getFreeParameters(), simpleValue.getOrder(), all);
1010
1011    }
1012
1013    /** {@inheritDoc}
1014     * @exception DimensionMismatchException if number of free parameters
1015     * or orders do not match
1016     * @since 3.2
1017     */
1018    @Override
1019    public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
1020                                                 final DerivativeStructure a2, final DerivativeStructure b2)
1021        throws DimensionMismatchException {
1022
1023        // compute an accurate value, taking care of cancellations
1024        final double accurateValue = LinearCombination.value(a1.getValue(), b1.getValue(),
1025                                                                  a2.getValue(), b2.getValue());
1026
1027        // compute a simple value, with all partial derivatives
1028        final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2));
1029
1030        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1031        final double[] all = simpleValue.getAllDerivatives();
1032        all[0] = accurateValue;
1033        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1034
1035    }
1036
1037    /** {@inheritDoc}
1038     * @exception DimensionMismatchException if number of free parameters
1039     * or orders do not match
1040     * @since 3.2
1041     */
1042    @Override
1043    public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
1044                                                 final double a2, final DerivativeStructure b2)
1045        throws DimensionMismatchException {
1046
1047        // compute an accurate value, taking care of cancellations
1048        final double accurateValue = LinearCombination.value(a1, b1.getValue(),
1049                                                                  a2, b2.getValue());
1050
1051        // compute a simple value, with all partial derivatives
1052        final DerivativeStructure simpleValue = b1.multiply(a1).add(b2.multiply(a2));
1053
1054        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1055        final double[] all = simpleValue.getAllDerivatives();
1056        all[0] = accurateValue;
1057        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1058
1059    }
1060
1061    /** {@inheritDoc}
1062     * @exception DimensionMismatchException if number of free parameters
1063     * or orders do not match
1064     * @since 3.2
1065     */
1066    @Override
1067    public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
1068                                                 final DerivativeStructure a2, final DerivativeStructure b2,
1069                                                 final DerivativeStructure a3, final DerivativeStructure b3)
1070        throws DimensionMismatchException {
1071
1072        // compute an accurate value, taking care of cancellations
1073        final double accurateValue = LinearCombination.value(a1.getValue(), b1.getValue(),
1074                                                                  a2.getValue(), b2.getValue(),
1075                                                                  a3.getValue(), b3.getValue());
1076
1077        // compute a simple value, with all partial derivatives
1078        final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3));
1079
1080        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1081        final double[] all = simpleValue.getAllDerivatives();
1082        all[0] = accurateValue;
1083        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1084
1085    }
1086
1087    /** {@inheritDoc}
1088     * @exception DimensionMismatchException if number of free parameters
1089     * or orders do not match
1090     * @since 3.2
1091     */
1092    @Override
1093    public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
1094                                                 final double a2, final DerivativeStructure b2,
1095                                                 final double a3, final DerivativeStructure b3)
1096        throws DimensionMismatchException {
1097
1098        // compute an accurate value, taking care of cancellations
1099        final double accurateValue = LinearCombination.value(a1, b1.getValue(),
1100                                                                  a2, b2.getValue(),
1101                                                                  a3, b3.getValue());
1102
1103        // compute a simple value, with all partial derivatives
1104        final DerivativeStructure simpleValue = b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3));
1105
1106        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1107        final double[] all = simpleValue.getAllDerivatives();
1108        all[0] = accurateValue;
1109        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1110
1111    }
1112
1113    /** {@inheritDoc}
1114     * @exception DimensionMismatchException if number of free parameters
1115     * or orders do not match
1116     * @since 3.2
1117     */
1118    @Override
1119    public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
1120                                                 final DerivativeStructure a2, final DerivativeStructure b2,
1121                                                 final DerivativeStructure a3, final DerivativeStructure b3,
1122                                                 final DerivativeStructure a4, final DerivativeStructure b4)
1123        throws DimensionMismatchException {
1124
1125        // compute an accurate value, taking care of cancellations
1126        final double accurateValue = LinearCombination.value(a1.getValue(), b1.getValue(),
1127                                                                  a2.getValue(), b2.getValue(),
1128                                                                  a3.getValue(), b3.getValue(),
1129                                                                  a4.getValue(), b4.getValue());
1130
1131        // compute a simple value, with all partial derivatives
1132        final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3)).add(a4.multiply(b4));
1133
1134        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1135        final double[] all = simpleValue.getAllDerivatives();
1136        all[0] = accurateValue;
1137        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1138
1139    }
1140
1141    /** {@inheritDoc}
1142     * @exception DimensionMismatchException if number of free parameters
1143     * or orders do not match
1144     * @since 3.2
1145     */
1146    @Override
1147    public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
1148                                                 final double a2, final DerivativeStructure b2,
1149                                                 final double a3, final DerivativeStructure b3,
1150                                                 final double a4, final DerivativeStructure b4)
1151        throws DimensionMismatchException {
1152
1153        // compute an accurate value, taking care of cancellations
1154        final double accurateValue = LinearCombination.value(a1, b1.getValue(),
1155                                                                  a2, b2.getValue(),
1156                                                                  a3, b3.getValue(),
1157                                                                  a4, b4.getValue());
1158
1159        // compute a simple value, with all partial derivatives
1160        final DerivativeStructure simpleValue = b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3)).add(b4.multiply(a4));
1161
1162        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1163        final double[] all = simpleValue.getAllDerivatives();
1164        all[0] = accurateValue;
1165        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1166
1167    }
1168
1169    /**
1170     * Test for the equality of two derivative structures.
1171     * <p>
1172     * Derivative structures are considered equal if they have the same number
1173     * of free parameters, the same derivation order, and the same derivatives.
1174     * </p>
1175     * @param other Object to test for equality to this
1176     * @return true if two derivative structures are equal
1177     * @since 3.2
1178     */
1179    @Override
1180    public boolean equals(Object other) {
1181
1182        if (this == other) {
1183            return true;
1184        }
1185
1186        if (other instanceof DerivativeStructure) {
1187            final DerivativeStructure rhs = (DerivativeStructure)other;
1188            return (getFreeParameters() == rhs.getFreeParameters()) &&
1189                   (getOrder() == rhs.getOrder()) &&
1190                   MathArrays.equals(data, rhs.data);
1191        }
1192
1193        return false;
1194
1195    }
1196
1197    /**
1198     * Get a hashCode for the derivative structure.
1199     * @return a hash code value for this object
1200     * @since 3.2
1201     */
1202    @Override
1203    public int hashCode() {
1204        return 227 + 229 * getFreeParameters() + 233 * getOrder() + 239 * MathUtils.hash(data);
1205    }
1206
1207    /**
1208     * Replace the instance with a data transfer object for serialization.
1209     * @return data transfer object that will be serialized
1210     */
1211    private Object writeReplace() {
1212        return new DataTransferObject(compiler.getFreeParameters(), compiler.getOrder(), data);
1213    }
1214
1215    /** Internal class used only for serialization. */
1216    private static class DataTransferObject implements Serializable {
1217
1218        /** Serializable UID. */
1219        private static final long serialVersionUID = 20120730L;
1220
1221        /** Number of variables.
1222         * @serial
1223         */
1224        private final int variables;
1225
1226        /** Derivation order.
1227         * @serial
1228         */
1229        private final int order;
1230
1231        /** Partial derivatives.
1232         * @serial
1233         */
1234        private final double[] data;
1235
1236        /** Simple constructor.
1237         * @param variables number of variables
1238         * @param order derivation order
1239         * @param data partial derivatives
1240         */
1241        DataTransferObject(final int variables, final int order, final double[] data) {
1242            this.variables = variables;
1243            this.order     = order;
1244            this.data      = data;
1245        }
1246
1247        /** Replace the deserialized data transfer object with a {@link DerivativeStructure}.
1248         * @return replacement {@link DerivativeStructure}
1249         */
1250        private Object readResolve() {
1251            return new DerivativeStructure(variables, order, data);
1252        }
1253
1254    }
1255
1256}