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