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.math3.analysis.differentiation;
018
019import java.util.ArrayList;
020import java.util.Arrays;
021import java.util.List;
022import java.util.concurrent.atomic.AtomicReference;
023
024import org.apache.commons.math3.exception.DimensionMismatchException;
025import org.apache.commons.math3.exception.MathArithmeticException;
026import org.apache.commons.math3.exception.MathInternalError;
027import org.apache.commons.math3.exception.NotPositiveException;
028import org.apache.commons.math3.exception.NumberIsTooLargeException;
029import org.apache.commons.math3.util.CombinatoricsUtils;
030import org.apache.commons.math3.util.FastMath;
031import org.apache.commons.math3.util.MathArrays;
032
033/** Class holding "compiled" computation rules for derivative structures.
034 * <p>This class implements the computation rules described in Dan Kalman's paper <a
035 * href="http://www1.american.edu/cas/mathstat/People/kalman/pdffiles/mmgautodiff.pdf">Doubly
036 * Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75,
037 * no. 3, June 2002. However, in order to avoid performances bottlenecks, the recursive
038 * rules are "compiled" once in an unfold form. This class does this recursion unrolling
039 * and stores the computation rules as simple loops with pre-computed indirection arrays.</p>
040 * <p>
041 * This class maps all derivative computation into single dimension arrays that hold the
042 * value and partial derivatives. The class does not hold these arrays, which remains under
043 * the responsibility of the caller. For each combination of number of free parameters and
044 * derivation order, only one compiler is necessary, and this compiler will be used to
045 * perform computations on all arrays provided to it, which can represent hundreds or
046 * thousands of different parameters kept together with all theur partial derivatives.
047 * </p>
048 * <p>
049 * The arrays on which compilers operate contain only the partial derivatives together
050 * with the 0<sup>th</sup> derivative, i.e. the value. The partial derivatives are stored in
051 * a compiler-specific order, which can be retrieved using methods {@link
052 * #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} and {@link
053 * #getPartialDerivativeOrders(int)}. The value is guaranteed to be stored as the first element
054 * (i.e. the {@link #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} method returns
055 * 0 when called with 0 for all derivation orders and {@link #getPartialDerivativeOrders(int)
056 * getPartialDerivativeOrders} returns an array filled with 0 when called with 0 as the index).
057 * </p>
058 * <p>
059 * Note that the ordering changes with number of parameters and derivation order. For example
060 * given 2 parameters x and y, df/dy is stored at index 2 when derivation order is set to 1 (in
061 * this case the array has three elements: f, df/dx and df/dy). If derivation order is set to
062 * 2, then df/dy will be stored at index 3 (in this case the array has six elements: f, df/dx,
063 * df/dxdx, df/dy, df/dxdy and df/dydy).
064 * </p>
065 * <p>
066 * Given this structure, users can perform some simple operations like adding, subtracting
067 * or multiplying constants and negating the elements by themselves, knowing if they want to
068 * mutate their array or create a new array. These simple operations are not provided by
069 * the compiler. The compiler provides only the more complex operations between several arrays.
070 * </p>
071 * <p>This class is mainly used as the engine for scalar variable {@link DerivativeStructure}.
072 * It can also be used directly to hold several variables in arrays for more complex data
073 * structures. User can for example store a vector of n variables depending on three x, y
074 * and z free parameters in one array as follows:</p> <pre>
075 *   // parameter 0 is x, parameter 1 is y, parameter 2 is z
076 *   int parameters = 3;
077 *   DSCompiler compiler = DSCompiler.getCompiler(parameters, order);
078 *   int size = compiler.getSize();
079 *
080 *   // pack all elements in a single array
081 *   double[] array = new double[n * size];
082 *   for (int i = 0; i &lt; n; ++i) {
083 *
084 *     // we know value is guaranteed to be the first element
085 *     array[i * size] = v[i];
086 *
087 *     // we don't know where first derivatives are stored, so we ask the compiler
088 *     array[i * size + compiler.getPartialDerivativeIndex(1, 0, 0) = dvOnDx[i][0];
089 *     array[i * size + compiler.getPartialDerivativeIndex(0, 1, 0) = dvOnDy[i][0];
090 *     array[i * size + compiler.getPartialDerivativeIndex(0, 0, 1) = dvOnDz[i][0];
091 *
092 *     // we let all higher order derivatives set to 0
093 *
094 *   }
095 * </pre>
096 * <p>Then in another function, user can perform some operations on all elements stored
097 * in the single array, such as a simple product of all variables:</p> <pre>
098 *   // compute the product of all elements
099 *   double[] product = new double[size];
100 *   prod[0] = 1.0;
101 *   for (int i = 0; i &lt; n; ++i) {
102 *     double[] tmp = product.clone();
103 *     compiler.multiply(tmp, 0, array, i * size, product, 0);
104 *   }
105 *
106 *   // value
107 *   double p = product[0];
108 *
109 *   // first derivatives
110 *   double dPdX = product[compiler.getPartialDerivativeIndex(1, 0, 0)];
111 *   double dPdY = product[compiler.getPartialDerivativeIndex(0, 1, 0)];
112 *   double dPdZ = product[compiler.getPartialDerivativeIndex(0, 0, 1)];
113 *
114 *   // cross derivatives (assuming order was at least 2)
115 *   double dPdXdX = product[compiler.getPartialDerivativeIndex(2, 0, 0)];
116 *   double dPdXdY = product[compiler.getPartialDerivativeIndex(1, 1, 0)];
117 *   double dPdXdZ = product[compiler.getPartialDerivativeIndex(1, 0, 1)];
118 *   double dPdYdY = product[compiler.getPartialDerivativeIndex(0, 2, 0)];
119 *   double dPdYdZ = product[compiler.getPartialDerivativeIndex(0, 1, 1)];
120 *   double dPdZdZ = product[compiler.getPartialDerivativeIndex(0, 0, 2)];
121 * </pre>
122 * @see DerivativeStructure
123 * @since 3.1
124 */
125public class DSCompiler {
126
127    /** Array of all compilers created so far. */
128    private static AtomicReference<DSCompiler[][]> compilers =
129            new AtomicReference<DSCompiler[][]>(null);
130
131    /** Number of free parameters. */
132    private final int parameters;
133
134    /** Derivation order. */
135    private final int order;
136
137    /** Number of partial derivatives (including the single 0 order derivative element). */
138    private final int[][] sizes;
139
140    /** Indirection array for partial derivatives. */
141    private final int[][] derivativesIndirection;
142
143    /** Indirection array of the lower derivative elements. */
144    private final int[] lowerIndirection;
145
146    /** Indirection arrays for multiplication. */
147    private final int[][][] multIndirection;
148
149    /** Indirection arrays for function composition. */
150    private final int[][][] compIndirection;
151
152    /** Private constructor, reserved for the factory method {@link #getCompiler(int, int)}.
153     * @param parameters number of free parameters
154     * @param order derivation order
155     * @param valueCompiler compiler for the value part
156     * @param derivativeCompiler compiler for the derivative part
157     * @throws NumberIsTooLargeException if order is too large
158     */
159    private DSCompiler(final int parameters, final int order,
160                       final DSCompiler valueCompiler, final DSCompiler derivativeCompiler)
161        throws NumberIsTooLargeException {
162
163        this.parameters = parameters;
164        this.order      = order;
165        this.sizes      = compileSizes(parameters, order, valueCompiler);
166        this.derivativesIndirection =
167                compileDerivativesIndirection(parameters, order,
168                                              valueCompiler, derivativeCompiler);
169        this.lowerIndirection =
170                compileLowerIndirection(parameters, order,
171                                        valueCompiler, derivativeCompiler);
172        this.multIndirection =
173                compileMultiplicationIndirection(parameters, order,
174                                                 valueCompiler, derivativeCompiler, lowerIndirection);
175        this.compIndirection =
176                compileCompositionIndirection(parameters, order,
177                                              valueCompiler, derivativeCompiler,
178                                              sizes, derivativesIndirection);
179
180    }
181
182    /** Get the compiler for number of free parameters and order.
183     * @param parameters number of free parameters
184     * @param order derivation order
185     * @return cached rules set
186     * @throws NumberIsTooLargeException if order is too large
187     */
188    public static DSCompiler getCompiler(int parameters, int order)
189        throws NumberIsTooLargeException {
190
191        // get the cached compilers
192        final DSCompiler[][] cache = compilers.get();
193        if (cache != null && cache.length > parameters &&
194            cache[parameters].length > order && cache[parameters][order] != null) {
195            // the compiler has already been created
196            return cache[parameters][order];
197        }
198
199        // we need to create more compilers
200        final int maxParameters = FastMath.max(parameters, cache == null ? 0 : cache.length);
201        final int maxOrder      = FastMath.max(order,     cache == null ? 0 : cache[0].length);
202        final DSCompiler[][] newCache = new DSCompiler[maxParameters + 1][maxOrder + 1];
203
204        if (cache != null) {
205            // preserve the already created compilers
206            for (int i = 0; i < cache.length; ++i) {
207                System.arraycopy(cache[i], 0, newCache[i], 0, cache[i].length);
208            }
209        }
210
211        // create the array in increasing diagonal order
212        for (int diag = 0; diag <= parameters + order; ++diag) {
213            for (int o = FastMath.max(0, diag - parameters); o <= FastMath.min(order, diag); ++o) {
214                final int p = diag - o;
215                if (newCache[p][o] == null) {
216                    final DSCompiler valueCompiler      = (p == 0) ? null : newCache[p - 1][o];
217                    final DSCompiler derivativeCompiler = (o == 0) ? null : newCache[p][o - 1];
218                    newCache[p][o] = new DSCompiler(p, o, valueCompiler, derivativeCompiler);
219                }
220            }
221        }
222
223        // atomically reset the cached compilers array
224        compilers.compareAndSet(cache, newCache);
225
226        return newCache[parameters][order];
227
228    }
229
230    /** Compile the sizes array.
231     * @param parameters number of free parameters
232     * @param order derivation order
233     * @param valueCompiler compiler for the value part
234     * @return sizes array
235     */
236    private static int[][] compileSizes(final int parameters, final int order,
237                                        final DSCompiler valueCompiler) {
238
239        final int[][] sizes = new int[parameters + 1][order + 1];
240        if (parameters == 0) {
241            Arrays.fill(sizes[0], 1);
242        } else {
243            System.arraycopy(valueCompiler.sizes, 0, sizes, 0, parameters);
244            sizes[parameters][0] = 1;
245            for (int i = 0; i < order; ++i) {
246                sizes[parameters][i + 1] = sizes[parameters][i] + sizes[parameters - 1][i + 1];
247            }
248        }
249
250        return sizes;
251
252    }
253
254    /** Compile the derivatives indirection array.
255     * @param parameters number of free parameters
256     * @param order derivation order
257     * @param valueCompiler compiler for the value part
258     * @param derivativeCompiler compiler for the derivative part
259     * @return derivatives indirection array
260     */
261    private static int[][] compileDerivativesIndirection(final int parameters, final int order,
262                                                      final DSCompiler valueCompiler,
263                                                      final DSCompiler derivativeCompiler) {
264
265        if (parameters == 0 || order == 0) {
266            return new int[1][parameters];
267        }
268
269        final int vSize = valueCompiler.derivativesIndirection.length;
270        final int dSize = derivativeCompiler.derivativesIndirection.length;
271        final int[][] derivativesIndirection = new int[vSize + dSize][parameters];
272
273        // set up the indices for the value part
274        for (int i = 0; i < vSize; ++i) {
275            // copy the first indices, the last one remaining set to 0
276            System.arraycopy(valueCompiler.derivativesIndirection[i], 0,
277                             derivativesIndirection[i], 0,
278                             parameters - 1);
279        }
280
281        // set up the indices for the derivative part
282        for (int i = 0; i < dSize; ++i) {
283
284            // copy the indices
285            System.arraycopy(derivativeCompiler.derivativesIndirection[i], 0,
286                             derivativesIndirection[vSize + i], 0,
287                             parameters);
288
289            // increment the derivation order for the last parameter
290            derivativesIndirection[vSize + i][parameters - 1]++;
291
292        }
293
294        return derivativesIndirection;
295
296    }
297
298    /** Compile the lower derivatives indirection array.
299     * <p>
300     * This indirection array contains the indices of all elements
301     * except derivatives for last derivation order.
302     * </p>
303     * @param parameters number of free parameters
304     * @param order derivation order
305     * @param valueCompiler compiler for the value part
306     * @param derivativeCompiler compiler for the derivative part
307     * @return lower derivatives indirection array
308     */
309    private static int[] compileLowerIndirection(final int parameters, final int order,
310                                              final DSCompiler valueCompiler,
311                                              final DSCompiler derivativeCompiler) {
312
313        if (parameters == 0 || order <= 1) {
314            return new int[] { 0 };
315        }
316
317        // this is an implementation of definition 6 in Dan Kalman's paper.
318        final int vSize = valueCompiler.lowerIndirection.length;
319        final int dSize = derivativeCompiler.lowerIndirection.length;
320        final int[] lowerIndirection = new int[vSize + dSize];
321        System.arraycopy(valueCompiler.lowerIndirection, 0, lowerIndirection, 0, vSize);
322        for (int i = 0; i < dSize; ++i) {
323            lowerIndirection[vSize + i] = valueCompiler.getSize() + derivativeCompiler.lowerIndirection[i];
324        }
325
326        return lowerIndirection;
327
328    }
329
330    /** Compile the multiplication indirection array.
331     * <p>
332     * This indirection array contains the indices of all pairs of elements
333     * involved when computing a multiplication. This allows a straightforward
334     * loop-based multiplication (see {@link #multiply(double[], int, double[], int, double[], int)}).
335     * </p>
336     * @param parameters number of free parameters
337     * @param order derivation order
338     * @param valueCompiler compiler for the value part
339     * @param derivativeCompiler compiler for the derivative part
340     * @param lowerIndirection lower derivatives indirection array
341     * @return multiplication indirection array
342     */
343    private static int[][][] compileMultiplicationIndirection(final int parameters, final int order,
344                                                           final DSCompiler valueCompiler,
345                                                           final DSCompiler derivativeCompiler,
346                                                           final int[] lowerIndirection) {
347
348        if ((parameters == 0) || (order == 0)) {
349            return new int[][][] { { { 1, 0, 0 } } };
350        }
351
352        // this is an implementation of definition 3 in Dan Kalman's paper.
353        final int vSize = valueCompiler.multIndirection.length;
354        final int dSize = derivativeCompiler.multIndirection.length;
355        final int[][][] multIndirection = new int[vSize + dSize][][];
356
357        System.arraycopy(valueCompiler.multIndirection, 0, multIndirection, 0, vSize);
358
359        for (int i = 0; i < dSize; ++i) {
360            final int[][] dRow = derivativeCompiler.multIndirection[i];
361            List<int[]> row = new ArrayList<int[]>(dRow.length * 2);
362            for (int j = 0; j < dRow.length; ++j) {
363                row.add(new int[] { dRow[j][0], lowerIndirection[dRow[j][1]], vSize + dRow[j][2] });
364                row.add(new int[] { dRow[j][0], vSize + dRow[j][1], lowerIndirection[dRow[j][2]] });
365            }
366
367            // combine terms with similar derivation orders
368            final List<int[]> combined = new ArrayList<int[]>(row.size());
369            for (int j = 0; j < row.size(); ++j) {
370                final int[] termJ = row.get(j);
371                if (termJ[0] > 0) {
372                    for (int k = j + 1; k < row.size(); ++k) {
373                        final int[] termK = row.get(k);
374                        if (termJ[1] == termK[1] && termJ[2] == termK[2]) {
375                            // combine termJ and termK
376                            termJ[0] += termK[0];
377                            // make sure we will skip termK later on in the outer loop
378                            termK[0] = 0;
379                        }
380                    }
381                    combined.add(termJ);
382                }
383            }
384
385            multIndirection[vSize + i] = combined.toArray(new int[combined.size()][]);
386
387        }
388
389        return multIndirection;
390
391    }
392
393    /** Compile the function composition indirection array.
394     * <p>
395     * This indirection array contains the indices of all sets of elements
396     * involved when computing a composition. This allows a straightforward
397     * loop-based composition (see {@link #compose(double[], int, double[], double[], int)}).
398     * </p>
399     * @param parameters number of free parameters
400     * @param order derivation order
401     * @param valueCompiler compiler for the value part
402     * @param derivativeCompiler compiler for the derivative part
403     * @param sizes sizes array
404     * @param derivativesIndirection derivatives indirection array
405     * @return multiplication indirection array
406     * @throws NumberIsTooLargeException if order is too large
407     */
408    private static int[][][] compileCompositionIndirection(final int parameters, final int order,
409                                                           final DSCompiler valueCompiler,
410                                                           final DSCompiler derivativeCompiler,
411                                                           final int[][] sizes,
412                                                           final int[][] derivativesIndirection)
413       throws NumberIsTooLargeException {
414
415        if ((parameters == 0) || (order == 0)) {
416            return new int[][][] { { { 1, 0 } } };
417        }
418
419        final int vSize = valueCompiler.compIndirection.length;
420        final int dSize = derivativeCompiler.compIndirection.length;
421        final int[][][] compIndirection = new int[vSize + dSize][][];
422
423        // the composition rules from the value part can be reused as is
424        System.arraycopy(valueCompiler.compIndirection, 0, compIndirection, 0, vSize);
425
426        // the composition rules for the derivative part are deduced by
427        // differentiation the rules from the underlying compiler once
428        // with respect to the parameter this compiler handles and the
429        // underlying one did not handle
430        for (int i = 0; i < dSize; ++i) {
431            List<int[]> row = new ArrayList<int[]>();
432            for (int[] term : derivativeCompiler.compIndirection[i]) {
433
434                // handle term p * f_k(g(x)) * g_l1(x) * g_l2(x) * ... * g_lp(x)
435
436                // derive the first factor in the term: f_k with respect to new parameter
437                int[] derivedTermF = new int[term.length + 1];
438                derivedTermF[0] = term[0];     // p
439                derivedTermF[1] = term[1] + 1; // f_(k+1)
440                int[] orders = new int[parameters];
441                orders[parameters - 1] = 1;
442                derivedTermF[term.length] = getPartialDerivativeIndex(parameters, order, sizes, orders);  // g_1
443                for (int j = 2; j < term.length; ++j) {
444                    // convert the indices as the mapping for the current order
445                    // is different from the mapping with one less order
446                    derivedTermF[j] = convertIndex(term[j], parameters,
447                                                   derivativeCompiler.derivativesIndirection,
448                                                   parameters, order, sizes);
449                }
450                Arrays.sort(derivedTermF, 2, derivedTermF.length);
451                row.add(derivedTermF);
452
453                // derive the various g_l
454                for (int l = 2; l < term.length; ++l) {
455                    int[] derivedTermG = new int[term.length];
456                    derivedTermG[0] = term[0];
457                    derivedTermG[1] = term[1];
458                    for (int j = 2; j < term.length; ++j) {
459                        // convert the indices as the mapping for the current order
460                        // is different from the mapping with one less order
461                        derivedTermG[j] = convertIndex(term[j], parameters,
462                                                       derivativeCompiler.derivativesIndirection,
463                                                       parameters, order, sizes);
464                        if (j == l) {
465                            // derive this term
466                            System.arraycopy(derivativesIndirection[derivedTermG[j]], 0, orders, 0, parameters);
467                            orders[parameters - 1]++;
468                            derivedTermG[j] = getPartialDerivativeIndex(parameters, order, sizes, orders);
469                        }
470                    }
471                    Arrays.sort(derivedTermG, 2, derivedTermG.length);
472                    row.add(derivedTermG);
473                }
474
475            }
476
477            // combine terms with similar derivation orders
478            final List<int[]> combined = new ArrayList<int[]>(row.size());
479            for (int j = 0; j < row.size(); ++j) {
480                final int[] termJ = row.get(j);
481                if (termJ[0] > 0) {
482                    for (int k = j + 1; k < row.size(); ++k) {
483                        final int[] termK = row.get(k);
484                        boolean equals = termJ.length == termK.length;
485                        for (int l = 1; equals && l < termJ.length; ++l) {
486                            equals &= termJ[l] == termK[l];
487                        }
488                        if (equals) {
489                            // combine termJ and termK
490                            termJ[0] += termK[0];
491                            // make sure we will skip termK later on in the outer loop
492                            termK[0] = 0;
493                        }
494                    }
495                    combined.add(termJ);
496                }
497            }
498
499            compIndirection[vSize + i] = combined.toArray(new int[combined.size()][]);
500
501        }
502
503        return compIndirection;
504
505    }
506
507    /** Get the index of a partial derivative in the array.
508     * <p>
509     * If all orders are set to 0, then the 0<sup>th</sup> order derivative
510     * is returned, which is the value of the function.
511     * </p>
512     * <p>The indices of derivatives are between 0 and {@link #getSize() getSize()} - 1.
513     * Their specific order is fixed for a given compiler, but otherwise not
514     * publicly specified. There are however some simple cases which have guaranteed
515     * indices:
516     * </p>
517     * <ul>
518     *   <li>the index of 0<sup>th</sup> order derivative is always 0</li>
519     *   <li>if there is only 1 {@link #getFreeParameters() free parameter}, then the
520     *   derivatives are sorted in increasing derivation order (i.e. f at index 0, df/dp
521     *   at index 1, d<sup>2</sup>f/dp<sup>2</sup> at index 2 ...
522     *   d<sup>k</sup>f/dp<sup>k</sup> at index k),</li>
523     *   <li>if the {@link #getOrder() derivation order} is 1, then the derivatives
524     *   are sorted in increasing free parameter order (i.e. f at index 0, df/dx<sub>1</sub>
525     *   at index 1, df/dx<sub>2</sub> at index 2 ... df/dx<sub>k</sub> at index k),</li>
526     *   <li>all other cases are not publicly specified</li>
527     * </ul>
528     * <p>
529     * This method is the inverse of method {@link #getPartialDerivativeOrders(int)}
530     * </p>
531     * @param orders derivation orders with respect to each parameter
532     * @return index of the partial derivative
533     * @exception DimensionMismatchException if the numbers of parameters does not
534     * match the instance
535     * @exception NumberIsTooLargeException if sum of derivation orders is larger
536     * than the instance limits
537     * @see #getPartialDerivativeOrders(int)
538     */
539    public int getPartialDerivativeIndex(final int ... orders)
540            throws DimensionMismatchException, NumberIsTooLargeException {
541
542        // safety check
543        if (orders.length != getFreeParameters()) {
544            throw new DimensionMismatchException(orders.length, getFreeParameters());
545        }
546
547        return getPartialDerivativeIndex(parameters, order, sizes, orders);
548
549    }
550
551    /** Get the index of a partial derivative in an array.
552     * @param parameters number of free parameters
553     * @param order derivation order
554     * @param sizes sizes array
555     * @param orders derivation orders with respect to each parameter
556     * (the lenght of this array must match the number of parameters)
557     * @return index of the partial derivative
558     * @exception NumberIsTooLargeException if sum of derivation orders is larger
559     * than the instance limits
560     */
561    private static int getPartialDerivativeIndex(final int parameters, final int order,
562                                                 final int[][] sizes, final int ... orders)
563        throws NumberIsTooLargeException {
564
565        // the value is obtained by diving into the recursive Dan Kalman's structure
566        // this is theorem 2 of his paper, with recursion replaced by iteration
567        int index     = 0;
568        int m         = order;
569        int ordersSum = 0;
570        for (int i = parameters - 1; i >= 0; --i) {
571
572            // derivative order for current free parameter
573            int derivativeOrder = orders[i];
574
575            // safety check
576            ordersSum += derivativeOrder;
577            if (ordersSum > order) {
578                throw new NumberIsTooLargeException(ordersSum, order, true);
579            }
580
581            while (derivativeOrder-- > 0) {
582                // as long as we differentiate according to current free parameter,
583                // we have to skip the value part and dive into the derivative part
584                // so we add the size of the value part to the base index
585                index += sizes[i][m--];
586            }
587
588        }
589
590        return index;
591
592    }
593
594    /** Convert an index from one (parameters, order) structure to another.
595     * @param index index of a partial derivative in source derivative structure
596     * @param srcP number of free parameters in source derivative structure
597     * @param srcDerivativesIndirection derivatives indirection array for the source
598     * derivative structure
599     * @param destP number of free parameters in destination derivative structure
600     * @param destO derivation order in destination derivative structure
601     * @param destSizes sizes array for the destination derivative structure
602     * @return index of the partial derivative with the <em>same</em> characteristics
603     * in destination derivative structure
604     * @throws NumberIsTooLargeException if order is too large
605     */
606    private static int convertIndex(final int index,
607                                    final int srcP, final int[][] srcDerivativesIndirection,
608                                    final int destP, final int destO, final int[][] destSizes)
609        throws NumberIsTooLargeException {
610        int[] orders = new int[destP];
611        System.arraycopy(srcDerivativesIndirection[index], 0, orders, 0, FastMath.min(srcP, destP));
612        return getPartialDerivativeIndex(destP, destO, destSizes, orders);
613    }
614
615    /** Get the derivation orders for a specific index in the array.
616     * <p>
617     * This method is the inverse of {@link #getPartialDerivativeIndex(int...)}.
618     * </p>
619     * @param index of the partial derivative
620     * @return orders derivation orders with respect to each parameter
621     * @see #getPartialDerivativeIndex(int...)
622     */
623    public int[] getPartialDerivativeOrders(final int index) {
624        return derivativesIndirection[index];
625    }
626
627    /** Get the number of free parameters.
628     * @return number of free parameters
629     */
630    public int getFreeParameters() {
631        return parameters;
632    }
633
634    /** Get the derivation order.
635     * @return derivation order
636     */
637    public int getOrder() {
638        return order;
639    }
640
641    /** Get the array size required for holding partial derivatives data.
642     * <p>
643     * This number includes the single 0 order derivative element, which is
644     * guaranteed to be stored in the first element of the array.
645     * </p>
646     * @return array size required for holding partial derivatives data
647     */
648    public int getSize() {
649        return sizes[parameters][order];
650    }
651
652    /** Compute linear combination.
653     * The derivative structure built will be a1 * ds1 + a2 * ds2
654     * @param a1 first scale factor
655     * @param c1 first base (unscaled) component
656     * @param offset1 offset of first operand in its array
657     * @param a2 second scale factor
658     * @param c2 second base (unscaled) component
659     * @param offset2 offset of second operand in its array
660     * @param result array where result must be stored (it may be
661     * one of the input arrays)
662     * @param resultOffset offset of the result in its array
663     */
664    public void linearCombination(final double a1, final double[] c1, final int offset1,
665                                  final double a2, final double[] c2, final int offset2,
666                                  final double[] result, final int resultOffset) {
667        for (int i = 0; i < getSize(); ++i) {
668            result[resultOffset + i] =
669                    MathArrays.linearCombination(a1, c1[offset1 + i], a2, c2[offset2 + i]);
670        }
671    }
672
673    /** Compute linear combination.
674     * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
675     * @param a1 first scale factor
676     * @param c1 first base (unscaled) component
677     * @param offset1 offset of first operand in its array
678     * @param a2 second scale factor
679     * @param c2 second base (unscaled) component
680     * @param offset2 offset of second operand in its array
681     * @param a3 third scale factor
682     * @param c3 third base (unscaled) component
683     * @param offset3 offset of third operand in its array
684     * @param result array where result must be stored (it may be
685     * one of the input arrays)
686     * @param resultOffset offset of the result in its array
687     */
688    public void linearCombination(final double a1, final double[] c1, final int offset1,
689                                  final double a2, final double[] c2, final int offset2,
690                                  final double a3, final double[] c3, final int offset3,
691                                  final double[] result, final int resultOffset) {
692        for (int i = 0; i < getSize(); ++i) {
693            result[resultOffset + i] =
694                    MathArrays.linearCombination(a1, c1[offset1 + i],
695                                                 a2, c2[offset2 + i],
696                                                 a3, c3[offset3 + i]);
697        }
698    }
699
700    /** Compute linear combination.
701     * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
702     * @param a1 first scale factor
703     * @param c1 first base (unscaled) component
704     * @param offset1 offset of first operand in its array
705     * @param a2 second scale factor
706     * @param c2 second base (unscaled) component
707     * @param offset2 offset of second operand in its array
708     * @param a3 third scale factor
709     * @param c3 third base (unscaled) component
710     * @param offset3 offset of third operand in its array
711     * @param a4 fourth scale factor
712     * @param c4 fourth base (unscaled) component
713     * @param offset4 offset of fourth operand in its array
714     * @param result array where result must be stored (it may be
715     * one of the input arrays)
716     * @param resultOffset offset of the result in its array
717     */
718    public void linearCombination(final double a1, final double[] c1, final int offset1,
719                                  final double a2, final double[] c2, final int offset2,
720                                  final double a3, final double[] c3, final int offset3,
721                                  final double a4, final double[] c4, final int offset4,
722                                  final double[] result, final int resultOffset) {
723        for (int i = 0; i < getSize(); ++i) {
724            result[resultOffset + i] =
725                    MathArrays.linearCombination(a1, c1[offset1 + i],
726                                                 a2, c2[offset2 + i],
727                                                 a3, c3[offset3 + i],
728                                                 a4, c4[offset4 + i]);
729        }
730    }
731
732    /** Perform addition of two derivative structures.
733     * @param lhs array holding left hand side of addition
734     * @param lhsOffset offset of the left hand side in its array
735     * @param rhs array right hand side of addition
736     * @param rhsOffset offset of the right hand side in its array
737     * @param result array where result must be stored (it may be
738     * one of the input arrays)
739     * @param resultOffset offset of the result in its array
740     */
741    public void add(final double[] lhs, final int lhsOffset,
742                    final double[] rhs, final int rhsOffset,
743                    final double[] result, final int resultOffset) {
744        for (int i = 0; i < getSize(); ++i) {
745            result[resultOffset + i] = lhs[lhsOffset + i] + rhs[rhsOffset + i];
746        }
747    }
748    /** Perform subtraction of two derivative structures.
749     * @param lhs array holding left hand side of subtraction
750     * @param lhsOffset offset of the left hand side in its array
751     * @param rhs array right hand side of subtraction
752     * @param rhsOffset offset of the right hand side in its array
753     * @param result array where result must be stored (it may be
754     * one of the input arrays)
755     * @param resultOffset offset of the result in its array
756     */
757    public void subtract(final double[] lhs, final int lhsOffset,
758                         final double[] rhs, final int rhsOffset,
759                         final double[] result, final int resultOffset) {
760        for (int i = 0; i < getSize(); ++i) {
761            result[resultOffset + i] = lhs[lhsOffset + i] - rhs[rhsOffset + i];
762        }
763    }
764
765    /** Perform multiplication of two derivative structures.
766     * @param lhs array holding left hand side of multiplication
767     * @param lhsOffset offset of the left hand side in its array
768     * @param rhs array right hand side of multiplication
769     * @param rhsOffset offset of the right hand side in its array
770     * @param result array where result must be stored (for
771     * multiplication the result array <em>cannot</em> be one of
772     * the input arrays)
773     * @param resultOffset offset of the result in its array
774     */
775    public void multiply(final double[] lhs, final int lhsOffset,
776                         final double[] rhs, final int rhsOffset,
777                         final double[] result, final int resultOffset) {
778        for (int i = 0; i < multIndirection.length; ++i) {
779            final int[][] mappingI = multIndirection[i];
780            double r = 0;
781            for (int j = 0; j < mappingI.length; ++j) {
782                r += mappingI[j][0] *
783                     lhs[lhsOffset + mappingI[j][1]] *
784                     rhs[rhsOffset + mappingI[j][2]];
785            }
786            result[resultOffset + i] = r;
787        }
788    }
789
790    /** Perform division of two derivative structures.
791     * @param lhs array holding left hand side of division
792     * @param lhsOffset offset of the left hand side in its array
793     * @param rhs array right hand side of division
794     * @param rhsOffset offset of the right hand side in its array
795     * @param result array where result must be stored (for
796     * division the result array <em>cannot</em> be one of
797     * the input arrays)
798     * @param resultOffset offset of the result in its array
799     */
800    public void divide(final double[] lhs, final int lhsOffset,
801                       final double[] rhs, final int rhsOffset,
802                       final double[] result, final int resultOffset) {
803        final double[] reciprocal = new double[getSize()];
804        pow(rhs, lhsOffset, -1, reciprocal, 0);
805        multiply(lhs, lhsOffset, reciprocal, 0, result, resultOffset);
806    }
807
808    /** Perform remainder of two derivative structures.
809     * @param lhs array holding left hand side of remainder
810     * @param lhsOffset offset of the left hand side in its array
811     * @param rhs array right hand side of remainder
812     * @param rhsOffset offset of the right hand side in its array
813     * @param result array where result must be stored (it may be
814     * one of the input arrays)
815     * @param resultOffset offset of the result in its array
816     */
817    public void remainder(final double[] lhs, final int lhsOffset,
818                          final double[] rhs, final int rhsOffset,
819                          final double[] result, final int resultOffset) {
820
821        // compute k such that lhs % rhs = lhs - k rhs
822        final double rem = FastMath.IEEEremainder(lhs[lhsOffset], rhs[rhsOffset]);
823        final double k   = FastMath.rint((lhs[lhsOffset] - rem) / rhs[rhsOffset]);
824
825        // set up value
826        result[resultOffset] = rem;
827
828        // set up partial derivatives
829        for (int i = 1; i < getSize(); ++i) {
830            result[resultOffset + i] = lhs[lhsOffset + i] - k * rhs[rhsOffset + i];
831        }
832
833    }
834
835    /** Compute power of a double to a derivative structure.
836     * @param a number to exponentiate
837     * @param operand array holding the power
838     * @param operandOffset offset of the power in its array
839     * @param result array where result must be stored (for
840     * power the result array <em>cannot</em> be the input
841     * array)
842     * @param resultOffset offset of the result in its array
843     * @since 3.3
844     */
845    public void pow(final double a,
846                    final double[] operand, final int operandOffset,
847                    final double[] result, final int resultOffset) {
848
849        // create the function value and derivatives
850        // [a^x, ln(a) a^x, ln(a)^2 a^x,, ln(a)^3 a^x, ... ]
851        final double[] function = new double[1 + order];
852        if (a == 0) {
853            if (operand[operandOffset] == 0) {
854                function[0] = 1;
855                double infinity = Double.POSITIVE_INFINITY;
856                for (int i = 1; i < function.length; ++i) {
857                    infinity = -infinity;
858                    function[i] = infinity;
859                }
860            } else if (operand[operandOffset] < 0) {
861                Arrays.fill(function, Double.NaN);
862            }
863        } else {
864            function[0] = FastMath.pow(a, operand[operandOffset]);
865            final double lnA = FastMath.log(a);
866            for (int i = 1; i < function.length; ++i) {
867                function[i] = lnA * function[i - 1];
868            }
869        }
870
871
872        // apply function composition
873        compose(operand, operandOffset, function, result, resultOffset);
874
875    }
876
877    /** Compute power of a derivative structure.
878     * @param operand array holding the operand
879     * @param operandOffset offset of the operand in its array
880     * @param p power to apply
881     * @param result array where result must be stored (for
882     * power the result array <em>cannot</em> be the input
883     * array)
884     * @param resultOffset offset of the result in its array
885     */
886    public void pow(final double[] operand, final int operandOffset, final double p,
887                    final double[] result, final int resultOffset) {
888
889        // create the function value and derivatives
890        // [x^p, px^(p-1), p(p-1)x^(p-2), ... ]
891        double[] function = new double[1 + order];
892        double xk = FastMath.pow(operand[operandOffset], p - order);
893        for (int i = order; i > 0; --i) {
894            function[i] = xk;
895            xk *= operand[operandOffset];
896        }
897        function[0] = xk;
898        double coefficient = p;
899        for (int i = 1; i <= order; ++i) {
900            function[i] *= coefficient;
901            coefficient *= p - i;
902        }
903
904        // apply function composition
905        compose(operand, operandOffset, function, result, resultOffset);
906
907    }
908
909    /** Compute integer power of a derivative structure.
910     * @param operand array holding the operand
911     * @param operandOffset offset of the operand in its array
912     * @param n power to apply
913     * @param result array where result must be stored (for
914     * power the result array <em>cannot</em> be the input
915     * array)
916     * @param resultOffset offset of the result in its array
917     */
918    public void pow(final double[] operand, final int operandOffset, final int n,
919                    final double[] result, final int resultOffset) {
920
921        if (n == 0) {
922            // special case, x^0 = 1 for all x
923            result[resultOffset] = 1.0;
924            Arrays.fill(result, resultOffset + 1, resultOffset + getSize(), 0);
925            return;
926        }
927
928        // create the power function value and derivatives
929        // [x^n, nx^(n-1), n(n-1)x^(n-2), ... ]
930        double[] function = new double[1 + order];
931
932        if (n > 0) {
933            // strictly positive power
934            final int maxOrder = FastMath.min(order, n);
935            double xk = FastMath.pow(operand[operandOffset], n - maxOrder);
936            for (int i = maxOrder; i > 0; --i) {
937                function[i] = xk;
938                xk *= operand[operandOffset];
939            }
940            function[0] = xk;
941        } else {
942            // strictly negative power
943            final double inv = 1.0 / operand[operandOffset];
944            double xk = FastMath.pow(inv, -n);
945            for (int i = 0; i <= order; ++i) {
946                function[i] = xk;
947                xk *= inv;
948            }
949        }
950
951        double coefficient = n;
952        for (int i = 1; i <= order; ++i) {
953            function[i] *= coefficient;
954            coefficient *= n - i;
955        }
956
957        // apply function composition
958        compose(operand, operandOffset, function, result, resultOffset);
959
960    }
961
962    /** Compute power of a derivative structure.
963     * @param x array holding the base
964     * @param xOffset offset of the base in its array
965     * @param y array holding the exponent
966     * @param yOffset offset of the exponent in its array
967     * @param result array where result must be stored (for
968     * power the result array <em>cannot</em> be the input
969     * array)
970     * @param resultOffset offset of the result in its array
971     */
972    public void pow(final double[] x, final int xOffset,
973                    final double[] y, final int yOffset,
974                    final double[] result, final int resultOffset) {
975        final double[] logX = new double[getSize()];
976        log(x, xOffset, logX, 0);
977        final double[] yLogX = new double[getSize()];
978        multiply(logX, 0, y, yOffset, yLogX, 0);
979        exp(yLogX, 0, result, resultOffset);
980    }
981
982    /** Compute n<sup>th</sup> root of a derivative structure.
983     * @param operand array holding the operand
984     * @param operandOffset offset of the operand in its array
985     * @param n order of the root
986     * @param result array where result must be stored (for
987     * n<sup>th</sup> root the result array <em>cannot</em> be the input
988     * array)
989     * @param resultOffset offset of the result in its array
990     */
991    public void rootN(final double[] operand, final int operandOffset, final int n,
992                      final double[] result, final int resultOffset) {
993
994        // create the function value and derivatives
995        // [x^(1/n), (1/n)x^((1/n)-1), (1-n)/n^2x^((1/n)-2), ... ]
996        double[] function = new double[1 + order];
997        double xk;
998        if (n == 2) {
999            function[0] = FastMath.sqrt(operand[operandOffset]);
1000            xk          = 0.5 / function[0];
1001        } else if (n == 3) {
1002            function[0] = FastMath.cbrt(operand[operandOffset]);
1003            xk          = 1.0 / (3.0 * function[0] * function[0]);
1004        } else {
1005            function[0] = FastMath.pow(operand[operandOffset], 1.0 / n);
1006            xk          = 1.0 / (n * FastMath.pow(function[0], n - 1));
1007        }
1008        final double nReciprocal = 1.0 / n;
1009        final double xReciprocal = 1.0 / operand[operandOffset];
1010        for (int i = 1; i <= order; ++i) {
1011            function[i] = xk;
1012            xk *= xReciprocal * (nReciprocal - i);
1013        }
1014
1015        // apply function composition
1016        compose(operand, operandOffset, function, result, resultOffset);
1017
1018    }
1019
1020    /** Compute exponential of a derivative structure.
1021     * @param operand array holding the operand
1022     * @param operandOffset offset of the operand in its array
1023     * @param result array where result must be stored (for
1024     * exponential the result array <em>cannot</em> be the input
1025     * array)
1026     * @param resultOffset offset of the result in its array
1027     */
1028    public void exp(final double[] operand, final int operandOffset,
1029                    final double[] result, final int resultOffset) {
1030
1031        // create the function value and derivatives
1032        double[] function = new double[1 + order];
1033        Arrays.fill(function, FastMath.exp(operand[operandOffset]));
1034
1035        // apply function composition
1036        compose(operand, operandOffset, function, result, resultOffset);
1037
1038    }
1039
1040    /** Compute exp(x) - 1 of a derivative structure.
1041     * @param operand array holding the operand
1042     * @param operandOffset offset of the operand in its array
1043     * @param result array where result must be stored (for
1044     * exponential the result array <em>cannot</em> be the input
1045     * array)
1046     * @param resultOffset offset of the result in its array
1047     */
1048    public void expm1(final double[] operand, final int operandOffset,
1049                      final double[] result, final int resultOffset) {
1050
1051        // create the function value and derivatives
1052        double[] function = new double[1 + order];
1053        function[0] = FastMath.expm1(operand[operandOffset]);
1054        Arrays.fill(function, 1, 1 + order, FastMath.exp(operand[operandOffset]));
1055
1056        // apply function composition
1057        compose(operand, operandOffset, function, result, resultOffset);
1058
1059    }
1060
1061    /** Compute natural logarithm of a derivative structure.
1062     * @param operand array holding the operand
1063     * @param operandOffset offset of the operand in its array
1064     * @param result array where result must be stored (for
1065     * logarithm the result array <em>cannot</em> be the input
1066     * array)
1067     * @param resultOffset offset of the result in its array
1068     */
1069    public void log(final double[] operand, final int operandOffset,
1070                    final double[] result, final int resultOffset) {
1071
1072        // create the function value and derivatives
1073        double[] function = new double[1 + order];
1074        function[0] = FastMath.log(operand[operandOffset]);
1075        if (order > 0) {
1076            double inv = 1.0 / operand[operandOffset];
1077            double xk  = inv;
1078            for (int i = 1; i <= order; ++i) {
1079                function[i] = xk;
1080                xk *= -i * inv;
1081            }
1082        }
1083
1084        // apply function composition
1085        compose(operand, operandOffset, function, result, resultOffset);
1086
1087    }
1088
1089    /** Computes shifted logarithm of a derivative structure.
1090     * @param operand array holding the operand
1091     * @param operandOffset offset of the operand in its array
1092     * @param result array where result must be stored (for
1093     * shifted logarithm the result array <em>cannot</em> be the input array)
1094     * @param resultOffset offset of the result in its array
1095     */
1096    public void log1p(final double[] operand, final int operandOffset,
1097                      final double[] result, final int resultOffset) {
1098
1099        // create the function value and derivatives
1100        double[] function = new double[1 + order];
1101        function[0] = FastMath.log1p(operand[operandOffset]);
1102        if (order > 0) {
1103            double inv = 1.0 / (1.0 + operand[operandOffset]);
1104            double xk  = inv;
1105            for (int i = 1; i <= order; ++i) {
1106                function[i] = xk;
1107                xk *= -i * inv;
1108            }
1109        }
1110
1111        // apply function composition
1112        compose(operand, operandOffset, function, result, resultOffset);
1113
1114    }
1115
1116    /** Computes base 10 logarithm of a derivative structure.
1117     * @param operand array holding the operand
1118     * @param operandOffset offset of the operand in its array
1119     * @param result array where result must be stored (for
1120     * base 10 logarithm the result array <em>cannot</em> be the input array)
1121     * @param resultOffset offset of the result in its array
1122     */
1123    public void log10(final double[] operand, final int operandOffset,
1124                      final double[] result, final int resultOffset) {
1125
1126        // create the function value and derivatives
1127        double[] function = new double[1 + order];
1128        function[0] = FastMath.log10(operand[operandOffset]);
1129        if (order > 0) {
1130            double inv = 1.0 / operand[operandOffset];
1131            double xk  = inv / FastMath.log(10.0);
1132            for (int i = 1; i <= order; ++i) {
1133                function[i] = xk;
1134                xk *= -i * inv;
1135            }
1136        }
1137
1138        // apply function composition
1139        compose(operand, operandOffset, function, result, resultOffset);
1140
1141    }
1142
1143    /** Compute cosine of a derivative structure.
1144     * @param operand array holding the operand
1145     * @param operandOffset offset of the operand in its array
1146     * @param result array where result must be stored (for
1147     * cosine the result array <em>cannot</em> be the input
1148     * array)
1149     * @param resultOffset offset of the result in its array
1150     */
1151    public void cos(final double[] operand, final int operandOffset,
1152                    final double[] result, final int resultOffset) {
1153
1154        // create the function value and derivatives
1155        double[] function = new double[1 + order];
1156        function[0] = FastMath.cos(operand[operandOffset]);
1157        if (order > 0) {
1158            function[1] = -FastMath.sin(operand[operandOffset]);
1159            for (int i = 2; i <= order; ++i) {
1160                function[i] = -function[i - 2];
1161            }
1162        }
1163
1164        // apply function composition
1165        compose(operand, operandOffset, function, result, resultOffset);
1166
1167    }
1168
1169    /** Compute sine of a derivative structure.
1170     * @param operand array holding the operand
1171     * @param operandOffset offset of the operand in its array
1172     * @param result array where result must be stored (for
1173     * sine the result array <em>cannot</em> be the input
1174     * array)
1175     * @param resultOffset offset of the result in its array
1176     */
1177    public void sin(final double[] operand, final int operandOffset,
1178                    final double[] result, final int resultOffset) {
1179
1180        // create the function value and derivatives
1181        double[] function = new double[1 + order];
1182        function[0] = FastMath.sin(operand[operandOffset]);
1183        if (order > 0) {
1184            function[1] = FastMath.cos(operand[operandOffset]);
1185            for (int i = 2; i <= order; ++i) {
1186                function[i] = -function[i - 2];
1187            }
1188        }
1189
1190        // apply function composition
1191        compose(operand, operandOffset, function, result, resultOffset);
1192
1193    }
1194
1195    /** Compute tangent of a derivative structure.
1196     * @param operand array holding the operand
1197     * @param operandOffset offset of the operand in its array
1198     * @param result array where result must be stored (for
1199     * tangent the result array <em>cannot</em> be the input
1200     * array)
1201     * @param resultOffset offset of the result in its array
1202     */
1203    public void tan(final double[] operand, final int operandOffset,
1204                    final double[] result, final int resultOffset) {
1205
1206        // create the function value and derivatives
1207        final double[] function = new double[1 + order];
1208        final double t = FastMath.tan(operand[operandOffset]);
1209        function[0] = t;
1210
1211        if (order > 0) {
1212
1213            // the nth order derivative of tan has the form:
1214            // dn(tan(x)/dxn = P_n(tan(x))
1215            // where P_n(t) is a degree n+1 polynomial with same parity as n+1
1216            // P_0(t) = t, P_1(t) = 1 + t^2, P_2(t) = 2 t (1 + t^2) ...
1217            // the general recurrence relation for P_n is:
1218            // P_n(x) = (1+t^2) P_(n-1)'(t)
1219            // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1220            final double[] p = new double[order + 2];
1221            p[1] = 1;
1222            final double t2 = t * t;
1223            for (int n = 1; n <= order; ++n) {
1224
1225                // update and evaluate polynomial P_n(t)
1226                double v = 0;
1227                p[n + 1] = n * p[n];
1228                for (int k = n + 1; k >= 0; k -= 2) {
1229                    v = v * t2 + p[k];
1230                    if (k > 2) {
1231                        p[k - 2] = (k - 1) * p[k - 1] + (k - 3) * p[k - 3];
1232                    } else if (k == 2) {
1233                        p[0] = p[1];
1234                    }
1235                }
1236                if ((n & 0x1) == 0) {
1237                    v *= t;
1238                }
1239
1240                function[n] = v;
1241
1242            }
1243        }
1244
1245        // apply function composition
1246        compose(operand, operandOffset, function, result, resultOffset);
1247
1248    }
1249
1250    /** Compute arc cosine of a derivative structure.
1251     * @param operand array holding the operand
1252     * @param operandOffset offset of the operand in its array
1253     * @param result array where result must be stored (for
1254     * arc cosine the result array <em>cannot</em> be the input
1255     * array)
1256     * @param resultOffset offset of the result in its array
1257     */
1258    public void acos(final double[] operand, final int operandOffset,
1259                    final double[] result, final int resultOffset) {
1260
1261        // create the function value and derivatives
1262        double[] function = new double[1 + order];
1263        final double x = operand[operandOffset];
1264        function[0] = FastMath.acos(x);
1265        if (order > 0) {
1266            // the nth order derivative of acos has the form:
1267            // dn(acos(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
1268            // where P_n(x) is a degree n-1 polynomial with same parity as n-1
1269            // P_1(x) = -1, P_2(x) = -x, P_3(x) = -2x^2 - 1 ...
1270            // the general recurrence relation for P_n is:
1271            // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
1272            // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1273            final double[] p = new double[order];
1274            p[0] = -1;
1275            final double x2    = x * x;
1276            final double f     = 1.0 / (1 - x2);
1277            double coeff = FastMath.sqrt(f);
1278            function[1] = coeff * p[0];
1279            for (int n = 2; n <= order; ++n) {
1280
1281                // update and evaluate polynomial P_n(x)
1282                double v = 0;
1283                p[n - 1] = (n - 1) * p[n - 2];
1284                for (int k = n - 1; k >= 0; k -= 2) {
1285                    v = v * x2 + p[k];
1286                    if (k > 2) {
1287                        p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
1288                    } else if (k == 2) {
1289                        p[0] = p[1];
1290                    }
1291                }
1292                if ((n & 0x1) == 0) {
1293                    v *= x;
1294                }
1295
1296                coeff *= f;
1297                function[n] = coeff * v;
1298
1299            }
1300        }
1301
1302        // apply function composition
1303        compose(operand, operandOffset, function, result, resultOffset);
1304
1305    }
1306
1307    /** Compute arc sine of a derivative structure.
1308     * @param operand array holding the operand
1309     * @param operandOffset offset of the operand in its array
1310     * @param result array where result must be stored (for
1311     * arc sine the result array <em>cannot</em> be the input
1312     * array)
1313     * @param resultOffset offset of the result in its array
1314     */
1315    public void asin(final double[] operand, final int operandOffset,
1316                    final double[] result, final int resultOffset) {
1317
1318        // create the function value and derivatives
1319        double[] function = new double[1 + order];
1320        final double x = operand[operandOffset];
1321        function[0] = FastMath.asin(x);
1322        if (order > 0) {
1323            // the nth order derivative of asin has the form:
1324            // dn(asin(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
1325            // where P_n(x) is a degree n-1 polynomial with same parity as n-1
1326            // P_1(x) = 1, P_2(x) = x, P_3(x) = 2x^2 + 1 ...
1327            // the general recurrence relation for P_n is:
1328            // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
1329            // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1330            final double[] p = new double[order];
1331            p[0] = 1;
1332            final double x2    = x * x;
1333            final double f     = 1.0 / (1 - x2);
1334            double coeff = FastMath.sqrt(f);
1335            function[1] = coeff * p[0];
1336            for (int n = 2; n <= order; ++n) {
1337
1338                // update and evaluate polynomial P_n(x)
1339                double v = 0;
1340                p[n - 1] = (n - 1) * p[n - 2];
1341                for (int k = n - 1; k >= 0; k -= 2) {
1342                    v = v * x2 + p[k];
1343                    if (k > 2) {
1344                        p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
1345                    } else if (k == 2) {
1346                        p[0] = p[1];
1347                    }
1348                }
1349                if ((n & 0x1) == 0) {
1350                    v *= x;
1351                }
1352
1353                coeff *= f;
1354                function[n] = coeff * v;
1355
1356            }
1357        }
1358
1359        // apply function composition
1360        compose(operand, operandOffset, function, result, resultOffset);
1361
1362    }
1363
1364    /** Compute arc tangent of a derivative structure.
1365     * @param operand array holding the operand
1366     * @param operandOffset offset of the operand in its array
1367     * @param result array where result must be stored (for
1368     * arc tangent the result array <em>cannot</em> be the input
1369     * array)
1370     * @param resultOffset offset of the result in its array
1371     */
1372    public void atan(final double[] operand, final int operandOffset,
1373                     final double[] result, final int resultOffset) {
1374
1375        // create the function value and derivatives
1376        double[] function = new double[1 + order];
1377        final double x = operand[operandOffset];
1378        function[0] = FastMath.atan(x);
1379        if (order > 0) {
1380            // the nth order derivative of atan has the form:
1381            // dn(atan(x)/dxn = Q_n(x) / (1 + x^2)^n
1382            // where Q_n(x) is a degree n-1 polynomial with same parity as n-1
1383            // Q_1(x) = 1, Q_2(x) = -2x, Q_3(x) = 6x^2 - 2 ...
1384            // the general recurrence relation for Q_n is:
1385            // Q_n(x) = (1+x^2) Q_(n-1)'(x) - 2(n-1) x Q_(n-1)(x)
1386            // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array
1387            final double[] q = new double[order];
1388            q[0] = 1;
1389            final double x2    = x * x;
1390            final double f     = 1.0 / (1 + x2);
1391            double coeff = f;
1392            function[1] = coeff * q[0];
1393            for (int n = 2; n <= order; ++n) {
1394
1395                // update and evaluate polynomial Q_n(x)
1396                double v = 0;
1397                q[n - 1] = -n * q[n - 2];
1398                for (int k = n - 1; k >= 0; k -= 2) {
1399                    v = v * x2 + q[k];
1400                    if (k > 2) {
1401                        q[k - 2] = (k - 1) * q[k - 1] + (k - 1 - 2 * n) * q[k - 3];
1402                    } else if (k == 2) {
1403                        q[0] = q[1];
1404                    }
1405                }
1406                if ((n & 0x1) == 0) {
1407                    v *= x;
1408                }
1409
1410                coeff *= f;
1411                function[n] = coeff * v;
1412
1413            }
1414        }
1415
1416        // apply function composition
1417        compose(operand, operandOffset, function, result, resultOffset);
1418
1419    }
1420
1421    /** Compute two arguments arc tangent of a derivative structure.
1422     * @param y array holding the first operand
1423     * @param yOffset offset of the first operand in its array
1424     * @param x array holding the second operand
1425     * @param xOffset offset of the second operand in its array
1426     * @param result array where result must be stored (for
1427     * two arguments arc tangent the result array <em>cannot</em>
1428     * be the input array)
1429     * @param resultOffset offset of the result in its array
1430     */
1431    public void atan2(final double[] y, final int yOffset,
1432                      final double[] x, final int xOffset,
1433                      final double[] result, final int resultOffset) {
1434
1435        // compute r = sqrt(x^2+y^2)
1436        double[] tmp1 = new double[getSize()];
1437        multiply(x, xOffset, x, xOffset, tmp1, 0);      // x^2
1438        double[] tmp2 = new double[getSize()];
1439        multiply(y, yOffset, y, yOffset, tmp2, 0);      // y^2
1440        add(tmp1, 0, tmp2, 0, tmp2, 0);                 // x^2 + y^2
1441        rootN(tmp2, 0, 2, tmp1, 0);                     // r = sqrt(x^2 + y^2)
1442
1443        if (x[xOffset] >= 0) {
1444
1445            // compute atan2(y, x) = 2 atan(y / (r + x))
1446            add(tmp1, 0, x, xOffset, tmp2, 0);          // r + x
1447            divide(y, yOffset, tmp2, 0, tmp1, 0);       // y /(r + x)
1448            atan(tmp1, 0, tmp2, 0);                     // atan(y / (r + x))
1449            for (int i = 0; i < tmp2.length; ++i) {
1450                result[resultOffset + i] = 2 * tmp2[i]; // 2 * atan(y / (r + x))
1451            }
1452
1453        } else {
1454
1455            // compute atan2(y, x) = +/- pi - 2 atan(y / (r - x))
1456            subtract(tmp1, 0, x, xOffset, tmp2, 0);     // r - x
1457            divide(y, yOffset, tmp2, 0, tmp1, 0);       // y /(r - x)
1458            atan(tmp1, 0, tmp2, 0);                     // atan(y / (r - x))
1459            result[resultOffset] =
1460                    ((tmp2[0] <= 0) ? -FastMath.PI : FastMath.PI) - 2 * tmp2[0]; // +/-pi - 2 * atan(y / (r - x))
1461            for (int i = 1; i < tmp2.length; ++i) {
1462                result[resultOffset + i] = -2 * tmp2[i]; // +/-pi - 2 * atan(y / (r - x))
1463            }
1464
1465        }
1466
1467        // fix value to take special cases (+0/+0, +0/-0, -0/+0, -0/-0, +/-infinity) correctly
1468        result[resultOffset] = FastMath.atan2(y[yOffset], x[xOffset]);
1469
1470    }
1471
1472    /** Compute hyperbolic cosine of a derivative structure.
1473     * @param operand array holding the operand
1474     * @param operandOffset offset of the operand in its array
1475     * @param result array where result must be stored (for
1476     * hyperbolic cosine the result array <em>cannot</em> be the input
1477     * array)
1478     * @param resultOffset offset of the result in its array
1479     */
1480    public void cosh(final double[] operand, final int operandOffset,
1481                     final double[] result, final int resultOffset) {
1482
1483        // create the function value and derivatives
1484        double[] function = new double[1 + order];
1485        function[0] = FastMath.cosh(operand[operandOffset]);
1486        if (order > 0) {
1487            function[1] = FastMath.sinh(operand[operandOffset]);
1488            for (int i = 2; i <= order; ++i) {
1489                function[i] = function[i - 2];
1490            }
1491        }
1492
1493        // apply function composition
1494        compose(operand, operandOffset, function, result, resultOffset);
1495
1496    }
1497
1498    /** Compute hyperbolic sine of a derivative structure.
1499     * @param operand array holding the operand
1500     * @param operandOffset offset of the operand in its array
1501     * @param result array where result must be stored (for
1502     * hyperbolic sine the result array <em>cannot</em> be the input
1503     * array)
1504     * @param resultOffset offset of the result in its array
1505     */
1506    public void sinh(final double[] operand, final int operandOffset,
1507                     final double[] result, final int resultOffset) {
1508
1509        // create the function value and derivatives
1510        double[] function = new double[1 + order];
1511        function[0] = FastMath.sinh(operand[operandOffset]);
1512        if (order > 0) {
1513            function[1] = FastMath.cosh(operand[operandOffset]);
1514            for (int i = 2; i <= order; ++i) {
1515                function[i] = function[i - 2];
1516            }
1517        }
1518
1519        // apply function composition
1520        compose(operand, operandOffset, function, result, resultOffset);
1521
1522    }
1523
1524    /** Compute hyperbolic tangent of a derivative structure.
1525     * @param operand array holding the operand
1526     * @param operandOffset offset of the operand in its array
1527     * @param result array where result must be stored (for
1528     * hyperbolic tangent the result array <em>cannot</em> be the input
1529     * array)
1530     * @param resultOffset offset of the result in its array
1531     */
1532    public void tanh(final double[] operand, final int operandOffset,
1533                     final double[] result, final int resultOffset) {
1534
1535        // create the function value and derivatives
1536        final double[] function = new double[1 + order];
1537        final double t = FastMath.tanh(operand[operandOffset]);
1538        function[0] = t;
1539
1540        if (order > 0) {
1541
1542            // the nth order derivative of tanh has the form:
1543            // dn(tanh(x)/dxn = P_n(tanh(x))
1544            // where P_n(t) is a degree n+1 polynomial with same parity as n+1
1545            // P_0(t) = t, P_1(t) = 1 - t^2, P_2(t) = -2 t (1 - t^2) ...
1546            // the general recurrence relation for P_n is:
1547            // P_n(x) = (1-t^2) P_(n-1)'(t)
1548            // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1549            final double[] p = new double[order + 2];
1550            p[1] = 1;
1551            final double t2 = t * t;
1552            for (int n = 1; n <= order; ++n) {
1553
1554                // update and evaluate polynomial P_n(t)
1555                double v = 0;
1556                p[n + 1] = -n * p[n];
1557                for (int k = n + 1; k >= 0; k -= 2) {
1558                    v = v * t2 + p[k];
1559                    if (k > 2) {
1560                        p[k - 2] = (k - 1) * p[k - 1] - (k - 3) * p[k - 3];
1561                    } else if (k == 2) {
1562                        p[0] = p[1];
1563                    }
1564                }
1565                if ((n & 0x1) == 0) {
1566                    v *= t;
1567                }
1568
1569                function[n] = v;
1570
1571            }
1572        }
1573
1574        // apply function composition
1575        compose(operand, operandOffset, function, result, resultOffset);
1576
1577    }
1578
1579    /** Compute inverse hyperbolic cosine of a derivative structure.
1580     * @param operand array holding the operand
1581     * @param operandOffset offset of the operand in its array
1582     * @param result array where result must be stored (for
1583     * inverse hyperbolic cosine the result array <em>cannot</em> be the input
1584     * array)
1585     * @param resultOffset offset of the result in its array
1586     */
1587    public void acosh(final double[] operand, final int operandOffset,
1588                     final double[] result, final int resultOffset) {
1589
1590        // create the function value and derivatives
1591        double[] function = new double[1 + order];
1592        final double x = operand[operandOffset];
1593        function[0] = FastMath.acosh(x);
1594        if (order > 0) {
1595            // the nth order derivative of acosh has the form:
1596            // dn(acosh(x)/dxn = P_n(x) / [x^2 - 1]^((2n-1)/2)
1597            // where P_n(x) is a degree n-1 polynomial with same parity as n-1
1598            // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 + 1 ...
1599            // the general recurrence relation for P_n is:
1600            // P_n(x) = (x^2-1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x)
1601            // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1602            final double[] p = new double[order];
1603            p[0] = 1;
1604            final double x2  = x * x;
1605            final double f   = 1.0 / (x2 - 1);
1606            double coeff = FastMath.sqrt(f);
1607            function[1] = coeff * p[0];
1608            for (int n = 2; n <= order; ++n) {
1609
1610                // update and evaluate polynomial P_n(x)
1611                double v = 0;
1612                p[n - 1] = (1 - n) * p[n - 2];
1613                for (int k = n - 1; k >= 0; k -= 2) {
1614                    v = v * x2 + p[k];
1615                    if (k > 2) {
1616                        p[k - 2] = (1 - k) * p[k - 1] + (k - 2 * n) * p[k - 3];
1617                    } else if (k == 2) {
1618                        p[0] = -p[1];
1619                    }
1620                }
1621                if ((n & 0x1) == 0) {
1622                    v *= x;
1623                }
1624
1625                coeff *= f;
1626                function[n] = coeff * v;
1627
1628            }
1629        }
1630
1631        // apply function composition
1632        compose(operand, operandOffset, function, result, resultOffset);
1633
1634    }
1635
1636    /** Compute inverse hyperbolic sine of a derivative structure.
1637     * @param operand array holding the operand
1638     * @param operandOffset offset of the operand in its array
1639     * @param result array where result must be stored (for
1640     * inverse hyperbolic sine the result array <em>cannot</em> be the input
1641     * array)
1642     * @param resultOffset offset of the result in its array
1643     */
1644    public void asinh(final double[] operand, final int operandOffset,
1645                     final double[] result, final int resultOffset) {
1646
1647        // create the function value and derivatives
1648        double[] function = new double[1 + order];
1649        final double x = operand[operandOffset];
1650        function[0] = FastMath.asinh(x);
1651        if (order > 0) {
1652            // the nth order derivative of asinh has the form:
1653            // dn(asinh(x)/dxn = P_n(x) / [x^2 + 1]^((2n-1)/2)
1654            // where P_n(x) is a degree n-1 polynomial with same parity as n-1
1655            // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 - 1 ...
1656            // the general recurrence relation for P_n is:
1657            // P_n(x) = (x^2+1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x)
1658            // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1659            final double[] p = new double[order];
1660            p[0] = 1;
1661            final double x2    = x * x;
1662            final double f     = 1.0 / (1 + x2);
1663            double coeff = FastMath.sqrt(f);
1664            function[1] = coeff * p[0];
1665            for (int n = 2; n <= order; ++n) {
1666
1667                // update and evaluate polynomial P_n(x)
1668                double v = 0;
1669                p[n - 1] = (1 - n) * p[n - 2];
1670                for (int k = n - 1; k >= 0; k -= 2) {
1671                    v = v * x2 + p[k];
1672                    if (k > 2) {
1673                        p[k - 2] = (k - 1) * p[k - 1] + (k - 2 * n) * p[k - 3];
1674                    } else if (k == 2) {
1675                        p[0] = p[1];
1676                    }
1677                }
1678                if ((n & 0x1) == 0) {
1679                    v *= x;
1680                }
1681
1682                coeff *= f;
1683                function[n] = coeff * v;
1684
1685            }
1686        }
1687
1688        // apply function composition
1689        compose(operand, operandOffset, function, result, resultOffset);
1690
1691    }
1692
1693    /** Compute inverse hyperbolic tangent of a derivative structure.
1694     * @param operand array holding the operand
1695     * @param operandOffset offset of the operand in its array
1696     * @param result array where result must be stored (for
1697     * inverse hyperbolic tangent the result array <em>cannot</em> be the input
1698     * array)
1699     * @param resultOffset offset of the result in its array
1700     */
1701    public void atanh(final double[] operand, final int operandOffset,
1702                      final double[] result, final int resultOffset) {
1703
1704        // create the function value and derivatives
1705        double[] function = new double[1 + order];
1706        final double x = operand[operandOffset];
1707        function[0] = FastMath.atanh(x);
1708        if (order > 0) {
1709            // the nth order derivative of atanh has the form:
1710            // dn(atanh(x)/dxn = Q_n(x) / (1 - x^2)^n
1711            // where Q_n(x) is a degree n-1 polynomial with same parity as n-1
1712            // Q_1(x) = 1, Q_2(x) = 2x, Q_3(x) = 6x^2 + 2 ...
1713            // the general recurrence relation for Q_n is:
1714            // Q_n(x) = (1-x^2) Q_(n-1)'(x) + 2(n-1) x Q_(n-1)(x)
1715            // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array
1716            final double[] q = new double[order];
1717            q[0] = 1;
1718            final double x2 = x * x;
1719            final double f  = 1.0 / (1 - x2);
1720            double coeff = f;
1721            function[1] = coeff * q[0];
1722            for (int n = 2; n <= order; ++n) {
1723
1724                // update and evaluate polynomial Q_n(x)
1725                double v = 0;
1726                q[n - 1] = n * q[n - 2];
1727                for (int k = n - 1; k >= 0; k -= 2) {
1728                    v = v * x2 + q[k];
1729                    if (k > 2) {
1730                        q[k - 2] = (k - 1) * q[k - 1] + (2 * n - k + 1) * q[k - 3];
1731                    } else if (k == 2) {
1732                        q[0] = q[1];
1733                    }
1734                }
1735                if ((n & 0x1) == 0) {
1736                    v *= x;
1737                }
1738
1739                coeff *= f;
1740                function[n] = coeff * v;
1741
1742            }
1743        }
1744
1745        // apply function composition
1746        compose(operand, operandOffset, function, result, resultOffset);
1747
1748    }
1749
1750    /** Compute composition of a derivative structure by a function.
1751     * @param operand array holding the operand
1752     * @param operandOffset offset of the operand in its array
1753     * @param f array of value and derivatives of the function at
1754     * the current point (i.e. at {@code operand[operandOffset]}).
1755     * @param result array where result must be stored (for
1756     * composition the result array <em>cannot</em> be the input
1757     * array)
1758     * @param resultOffset offset of the result in its array
1759     */
1760    public void compose(final double[] operand, final int operandOffset, final double[] f,
1761                        final double[] result, final int resultOffset) {
1762        for (int i = 0; i < compIndirection.length; ++i) {
1763            final int[][] mappingI = compIndirection[i];
1764            double r = 0;
1765            for (int j = 0; j < mappingI.length; ++j) {
1766                final int[] mappingIJ = mappingI[j];
1767                double product = mappingIJ[0] * f[mappingIJ[1]];
1768                for (int k = 2; k < mappingIJ.length; ++k) {
1769                    product *= operand[operandOffset + mappingIJ[k]];
1770                }
1771                r += product;
1772            }
1773            result[resultOffset + i] = r;
1774        }
1775    }
1776
1777    /** Evaluate Taylor expansion of a derivative structure.
1778     * @param ds array holding the derivative structure
1779     * @param dsOffset offset of the derivative structure in its array
1780     * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
1781     * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
1782     * @throws MathArithmeticException if factorials becomes too large
1783     */
1784    public double taylor(final double[] ds, final int dsOffset, final double ... delta)
1785       throws MathArithmeticException {
1786        double value = 0;
1787        for (int i = getSize() - 1; i >= 0; --i) {
1788            final int[] orders = getPartialDerivativeOrders(i);
1789            double term = ds[dsOffset + i];
1790            for (int k = 0; k < orders.length; ++k) {
1791                if (orders[k] > 0) {
1792                    try {
1793                        term *= FastMath.pow(delta[k], orders[k]) /
1794                        CombinatoricsUtils.factorial(orders[k]);
1795                    } catch (NotPositiveException e) {
1796                        // this cannot happen
1797                        throw new MathInternalError(e);
1798                    }
1799                }
1800            }
1801            value += term;
1802        }
1803        return value;
1804    }
1805
1806    /** Check rules set compatibility.
1807     * @param compiler other compiler to check against instance
1808     * @exception DimensionMismatchException if number of free parameters or orders are inconsistent
1809     */
1810    public void checkCompatibility(final DSCompiler compiler)
1811            throws DimensionMismatchException {
1812        if (parameters != compiler.parameters) {
1813            throw new DimensionMismatchException(parameters, compiler.parameters);
1814        }
1815        if (order != compiler.order) {
1816            throw new DimensionMismatchException(order, compiler.order);
1817        }
1818    }
1819
1820}