View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.analysis.differentiation;
18  
19  import java.util.ArrayList;
20  import java.util.Arrays;
21  import java.util.List;
22  import java.util.concurrent.atomic.AtomicReference;
23  
24  import org.apache.commons.math3.exception.DimensionMismatchException;
25  import org.apache.commons.math3.exception.MathArithmeticException;
26  import org.apache.commons.math3.exception.MathInternalError;
27  import org.apache.commons.math3.exception.NotPositiveException;
28  import org.apache.commons.math3.exception.NumberIsTooLargeException;
29  import org.apache.commons.math3.util.CombinatoricsUtils;
30  import org.apache.commons.math3.util.FastMath;
31  import org.apache.commons.math3.util.MathArrays;
32  
33  /** Class holding "compiled" computation rules for derivative structures.
34   * <p>This class implements the computation rules described in Dan Kalman's paper <a
35   * href="http://www1.american.edu/cas/mathstat/People/kalman/pdffiles/mmgautodiff.pdf">Doubly
36   * Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75,
37   * no. 3, June 2002. However, in order to avoid performances bottlenecks, the recursive
38   * rules are "compiled" once in an unfold form. This class does this recursion unrolling
39   * and stores the computation rules as simple loops with pre-computed indirection arrays.</p>
40   * <p>
41   * This class maps all derivative computation into single dimension arrays that hold the
42   * value and partial derivatives. The class does not hold these arrays, which remains under
43   * the responsibility of the caller. For each combination of number of free parameters and
44   * derivation order, only one compiler is necessary, and this compiler will be used to
45   * perform computations on all arrays provided to it, which can represent hundreds or
46   * thousands of different parameters kept together with all theur partial derivatives.
47   * </p>
48   * <p>
49   * The arrays on which compilers operate contain only the partial derivatives together
50   * with the 0<sup>th</sup> derivative, i.e. the value. The partial derivatives are stored in
51   * a compiler-specific order, which can be retrieved using methods {@link
52   * #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} and {@link
53   * #getPartialDerivativeOrders(int)}. The value is guaranteed to be stored as the first element
54   * (i.e. the {@link #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} method returns
55   * 0 when called with 0 for all derivation orders and {@link #getPartialDerivativeOrders(int)
56   * getPartialDerivativeOrders} returns an array filled with 0 when called with 0 as the index).
57   * </p>
58   * <p>
59   * Note that the ordering changes with number of parameters and derivation order. For example
60   * given 2 parameters x and y, df/dy is stored at index 2 when derivation order is set to 1 (in
61   * this case the array has three elements: f, df/dx and df/dy). If derivation order is set to
62   * 2, then df/dy will be stored at index 3 (in this case the array has six elements: f, df/dx,
63   * df/dxdx, df/dy, df/dxdy and df/dydy).
64   * </p>
65   * <p>
66   * Given this structure, users can perform some simple operations like adding, subtracting
67   * or multiplying constants and negating the elements by themselves, knowing if they want to
68   * mutate their array or create a new array. These simple operations are not provided by
69   * the compiler. The compiler provides only the more complex operations between several arrays.
70   * </p>
71   * <p>This class is mainly used as the engine for scalar variable {@link DerivativeStructure}.
72   * It can also be used directly to hold several variables in arrays for more complex data
73   * structures. User can for example store a vector of n variables depending on three x, y
74   * and z free parameters in one array as follows:</p> <pre>
75   *   // parameter 0 is x, parameter 1 is y, parameter 2 is z
76   *   int parameters = 3;
77   *   DSCompiler compiler = DSCompiler.getCompiler(parameters, order);
78   *   int size = compiler.getSize();
79   *
80   *   // pack all elements in a single array
81   *   double[] array = new double[n * size];
82   *   for (int i = 0; i &lt; n; ++i) {
83   *
84   *     // we know value is guaranteed to be the first element
85   *     array[i * size] = v[i];
86   *
87   *     // we don't know where first derivatives are stored, so we ask the compiler
88   *     array[i * size + compiler.getPartialDerivativeIndex(1, 0, 0) = dvOnDx[i][0];
89   *     array[i * size + compiler.getPartialDerivativeIndex(0, 1, 0) = dvOnDy[i][0];
90   *     array[i * size + compiler.getPartialDerivativeIndex(0, 0, 1) = dvOnDz[i][0];
91   *
92   *     // we let all higher order derivatives set to 0
93   *
94   *   }
95   * </pre>
96   * <p>Then in another function, user can perform some operations on all elements stored
97   * in the single array, such as a simple product of all variables:</p> <pre>
98   *   // compute the product of all elements
99   *   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  */
125 public 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 }