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    *      https://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  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.differentiation;
23  
24  import java.lang.reflect.Array;
25  import java.util.ArrayList;
26  import java.util.Arrays;
27  import java.util.List;
28  import java.util.concurrent.atomic.AtomicReference;
29  
30  import org.hipparchus.CalculusFieldElement;
31  import org.hipparchus.Field;
32  import org.hipparchus.exception.LocalizedCoreFormats;
33  import org.hipparchus.exception.MathIllegalArgumentException;
34  import org.hipparchus.exception.MathRuntimeException;
35  import org.hipparchus.util.CombinatoricsUtils;
36  import org.hipparchus.util.FastMath;
37  import org.hipparchus.util.FieldSinCos;
38  import org.hipparchus.util.FieldSinhCosh;
39  import org.hipparchus.util.MathArrays;
40  import org.hipparchus.util.MathUtils;
41  import org.hipparchus.util.SinCos;
42  import org.hipparchus.util.SinhCosh;
43  
44  /** Class holding "compiled" computation rules for derivative structures.
45   * <p>This class implements the computation rules described in Dan Kalman's paper <a
46   * href="http://www1.american.edu/cas/mathstat/People/kalman/pdffiles/mmgautodiff.pdf">Doubly
47   * Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75,
48   * no. 3, June 2002. However, in order to avoid performances bottlenecks, the recursive
49   * rules are "compiled" once in an unfold form. This class does this recursion unrolling
50   * and stores the computation rules as simple loops with pre-computed indirection arrays.</p>
51   * <p>
52   * This class maps all derivative computation into single dimension arrays that hold the
53   * value and partial derivatives. The class does not hold these arrays, which remains under
54   * the responsibility of the caller. For each combination of number of free parameters and
55   * derivation order, only one compiler is necessary, and this compiler will be used to
56   * perform computations on all arrays provided to it, which can represent hundreds or
57   * thousands of different parameters kept together with all their partial derivatives.
58   * </p>
59   * <p>
60   * The arrays on which compilers operate contain only the partial derivatives together
61   * with the 0<sup>th</sup> derivative, i.e. the value. The partial derivatives are stored in
62   * a compiler-specific order, which can be retrieved using methods {@link
63   * #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} and {@link
64   * #getPartialDerivativeOrders(int)}. The value is guaranteed to be stored as the first element
65   * (i.e. the {@link #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} method returns
66   * 0 when called with 0 for all derivation orders and {@link #getPartialDerivativeOrders(int)
67   * getPartialDerivativeOrders} returns an array filled with 0 when called with 0 as the index).
68   * </p>
69   * <p>
70   * Note that the ordering changes with number of parameters and derivation order. For example
71   * given 2 parameters x and y, df/dy is stored at index 2 when derivation order is set to 1 (in
72   * this case the array has three elements: f, df/dx and df/dy). If derivation order is set to
73   * 2, then df/dy will be stored at index 3 (in this case the array has six elements: f, df/dx,
74   * d²f/dxdx, df/dy, d²f/dxdy and d²f/dydy).
75   * </p>
76   * <p>
77   * Given this structure, users can perform some simple operations like adding, subtracting
78   * or multiplying constants and negating the elements by themselves, knowing if they want to
79   * mutate their array or create a new array. These simple operations are not provided by
80   * the compiler. The compiler provides only the more complex operations between several arrays.
81   * </p>
82   * <p>This class is mainly used as the engine for scalar variable {@link DerivativeStructure}.
83   * It can also be used directly to hold several variables in arrays for more complex data
84   * structures. User can for example store a vector of n variables depending on three x, y
85   * and z free parameters in one array as follows:</p> <pre>
86   *   // parameter 0 is x, parameter 1 is y, parameter 2 is z
87   *   int parameters = 3;
88   *   DSCompiler compiler = DSCompiler.getCompiler(parameters, order);
89   *   int size = compiler.getSize();
90   *
91   *   // pack all elements in a single array
92   *   double[] array = new double[n * size];
93   *   for (int i = 0; i &lt; n; ++i) {
94   *
95   *     // we know value is guaranteed to be the first element
96   *     array[i * size] = v[i];
97   *
98   *     // we don't know where first derivatives are stored, so we ask the compiler
99   *     array[i * size + compiler.getPartialDerivativeIndex(1, 0, 0) = dvOnDx[i][0];
100  *     array[i * size + compiler.getPartialDerivativeIndex(0, 1, 0) = dvOnDy[i][0];
101  *     array[i * size + compiler.getPartialDerivativeIndex(0, 0, 1) = dvOnDz[i][0];
102  *
103  *     // we let all higher order derivatives set to 0
104  *
105  *   }
106  * </pre>
107  * <p>Then in another function, user can perform some operations on all elements stored
108  * in the single array, such as a simple product of all variables:</p> <pre>
109  *   // compute the product of all elements
110  *   double[] product = new double[size];
111  *   prod[0] = 1.0;
112  *   for (int i = 0; i &lt; n; ++i) {
113  *     double[] tmp = product.clone();
114  *     compiler.multiply(tmp, 0, array, i * size, product, 0);
115  *   }
116  *
117  *   // value
118  *   double p = product[0];
119  *
120  *   // first derivatives
121  *   double dPdX = product[compiler.getPartialDerivativeIndex(1, 0, 0)];
122  *   double dPdY = product[compiler.getPartialDerivativeIndex(0, 1, 0)];
123  *   double dPdZ = product[compiler.getPartialDerivativeIndex(0, 0, 1)];
124  *
125  *   // cross derivatives (assuming order was at least 2)
126  *   double dPdXdX = product[compiler.getPartialDerivativeIndex(2, 0, 0)];
127  *   double dPdXdY = product[compiler.getPartialDerivativeIndex(1, 1, 0)];
128  *   double dPdXdZ = product[compiler.getPartialDerivativeIndex(1, 0, 1)];
129  *   double dPdYdY = product[compiler.getPartialDerivativeIndex(0, 2, 0)];
130  *   double dPdYdZ = product[compiler.getPartialDerivativeIndex(0, 1, 1)];
131  *   double dPdZdZ = product[compiler.getPartialDerivativeIndex(0, 0, 2)];
132  * </pre>
133  * @see DerivativeStructure
134  * @see FieldDerivativeStructure
135  */
136 public class DSCompiler {
137 
138     /** Array of all compilers created so far. */
139     private static AtomicReference<DSCompiler[][]> compilers = new AtomicReference<>(null);
140 
141     /** Number of free parameters. */
142     private final int parameters;
143 
144     /** Derivation order. */
145     private final int order;
146 
147     /** Number of partial derivatives (including the single 0 order derivative element). */
148     private final int[][] sizes;
149 
150     /** Orders array for partial derivatives. */
151     private final int[][] derivativesOrders;
152 
153     /** Sum of orders array for partial derivatives. */
154     private final int[] derivativesOrdersSum;
155 
156     /** Indirection array of the lower derivative elements. */
157     private final int[] lowerIndirection;
158 
159     /** Indirection arrays for multiplication. */
160     private final MultiplicationMapper[][] multIndirection;
161 
162     /** Indirection arrays for univariate function composition. */
163     private final UnivariateCompositionMapper[][] compIndirection;
164 
165     /** Indirection arrays for multivariate function rebasing. */
166     private final List<MultivariateCompositionMapper[][]> rebaseIndirection;
167 
168     /** Private constructor, reserved for the factory method {@link #getCompiler(int, int)}.
169      * @param parameters number of free parameters
170      * @param order derivation order
171      * @param valueCompiler compiler for the value part
172      * @param derivativeCompiler compiler for the derivative part
173      * @throws MathIllegalArgumentException if order is too large
174      */
175     private DSCompiler(final int parameters, final int order,
176                        final DSCompiler valueCompiler, final DSCompiler derivativeCompiler)
177         throws MathIllegalArgumentException {
178 
179         this.parameters           = parameters;
180         this.order                = order;
181         this.sizes                = compileSizes(parameters, order, valueCompiler);
182         this.derivativesOrders    = compileDerivativesOrders(parameters, order,
183                                                              valueCompiler, derivativeCompiler);
184         this.derivativesOrdersSum = compileDerivativesOrdersSum(derivativesOrders);
185         this.lowerIndirection     = compileLowerIndirection(parameters, order,
186                                                             valueCompiler, derivativeCompiler);
187         this.multIndirection      = compileMultiplicationIndirection(parameters, order,
188                                                                      valueCompiler, derivativeCompiler,
189                                                                      lowerIndirection);
190         this.compIndirection      = compileCompositionIndirection(parameters, order,
191                                                                   valueCompiler, derivativeCompiler,
192                                                                   sizes, derivativesOrders);
193 
194         this.rebaseIndirection = new ArrayList<>();
195     }
196 
197     /** Get the compiler for number of free parameters and order.
198      * @param parameters number of free parameters
199      * @param order derivation order
200      * @return cached rules set
201      * @throws MathIllegalArgumentException if order is too large
202      */
203     public static DSCompiler getCompiler(int parameters, int order)
204         throws MathIllegalArgumentException {
205 
206         // get the cached compilers
207         final DSCompiler[][] cache = compilers.get();
208         if (cache != null && cache.length > parameters &&
209             cache[parameters].length > order && cache[parameters][order] != null) {
210             // the compiler has already been created
211             return cache[parameters][order];
212         }
213 
214         // we need to create more compilers
215         final int maxParameters = FastMath.max(parameters, cache == null ? 0 : cache.length);
216         final int maxOrder      = FastMath.max(order,     cache == null ? 0 : cache[0].length);
217         final DSCompiler[][] newCache = new DSCompiler[maxParameters + 1][maxOrder + 1];
218 
219         if (cache != null) {
220             // preserve the already created compilers
221             for (int i = 0; i < cache.length; ++i) {
222                 System.arraycopy(cache[i], 0, newCache[i], 0, cache[i].length);
223             }
224         }
225 
226         // create the array in increasing diagonal order
227         for (int diag = 0; diag <= parameters + order; ++diag) {
228             for (int o = FastMath.max(0, diag - parameters); o <= FastMath.min(order, diag); ++o) {
229                 final int p = diag - o;
230                 if (newCache[p][o] == null) {
231                     final DSCompiler valueCompiler      = (p == 0) ? null : newCache[p - 1][o];
232                     final DSCompiler derivativeCompiler = (o == 0) ? null : newCache[p][o - 1];
233                     newCache[p][o] = new DSCompiler(p, o, valueCompiler, derivativeCompiler);
234                 }
235             }
236         }
237 
238         // atomically reset the cached compilers array
239         compilers.compareAndSet(cache, newCache);
240 
241         return newCache[parameters][order];
242 
243     }
244 
245     /** Compile the sizes array.
246      * @param parameters number of free parameters
247      * @param order derivation order
248      * @param valueCompiler compiler for the value part
249      * @return sizes array
250      */
251     private static int[][] compileSizes(final int parameters, final int order,
252                                         final DSCompiler valueCompiler) {
253 
254         final int[][] sizes = new int[parameters + 1][order + 1];
255         if (parameters == 0) {
256             Arrays.fill(sizes[0], 1);
257         } else {
258             System.arraycopy(valueCompiler.sizes, 0, sizes, 0, parameters);
259             sizes[parameters][0] = 1;
260             for (int i = 0; i < order; ++i) {
261                 sizes[parameters][i + 1] = sizes[parameters][i] + sizes[parameters - 1][i + 1];
262             }
263         }
264 
265         return sizes;
266 
267     }
268 
269     /** Compile the derivatives orders array.
270      * @param parameters number of free parameters
271      * @param order derivation order
272      * @param valueCompiler compiler for the value part
273      * @param derivativeCompiler compiler for the derivative part
274      * @return derivatives orders array
275      */
276     private static int[][] compileDerivativesOrders(final int parameters, final int order,
277                                                     final DSCompiler valueCompiler,
278                                                     final DSCompiler derivativeCompiler) {
279 
280         if (parameters == 0 || order == 0) {
281             return new int[1][parameters];
282         }
283 
284         final int vSize = valueCompiler.derivativesOrders.length;
285         final int dSize = derivativeCompiler.derivativesOrders.length;
286         final int[][] derivativesOrders = new int[vSize + dSize][parameters];
287 
288         // set up the indices for the value part
289         for (int i = 0; i < vSize; ++i) {
290             // copy the first indices, the last one remaining set to 0
291             System.arraycopy(valueCompiler.derivativesOrders[i], 0,
292                              derivativesOrders[i], 0,
293                              parameters - 1);
294         }
295 
296         // set up the indices for the derivative part
297         for (int i = 0; i < dSize; ++i) {
298 
299             // copy the indices
300             System.arraycopy(derivativeCompiler.derivativesOrders[i], 0,
301                              derivativesOrders[vSize + i], 0,
302                              parameters);
303 
304             // increment the derivation order for the last parameter
305             derivativesOrders[vSize + i][parameters - 1]++;
306 
307         }
308 
309         return derivativesOrders;
310 
311     }
312 
313     /** Compile the sum of orders array for partial derivatives.
314      * @param derivativesOrders orders array for partial derivatives
315      * @return sum of orders array for partial derivatives
316      */
317     private static int[] compileDerivativesOrdersSum(final int[][] derivativesOrders) {
318 
319         final int[] derivativesOrdersSum = new int[derivativesOrders.length];
320 
321         // locate the partial derivatives at order 1
322         for (int i = 0; i < derivativesOrdersSum.length; ++i) {
323             for (final int o : derivativesOrders[i]) {
324                 derivativesOrdersSum[i] += o;
325             }
326         }
327 
328         return derivativesOrdersSum;
329 
330     }
331 
332     /** Compile the lower derivatives indirection array.
333      * <p>
334      * This indirection array contains the indices of all elements
335      * except derivatives for last derivation order.
336      * </p>
337      * @param parameters number of free parameters
338      * @param order derivation order
339      * @param valueCompiler compiler for the value part
340      * @param derivativeCompiler compiler for the derivative part
341      * @return lower derivatives indirection array
342      */
343     private static int[] compileLowerIndirection(final int parameters, final int order,
344                                                  final DSCompiler valueCompiler,
345                                                  final DSCompiler derivativeCompiler) {
346 
347         if (parameters == 0 || order <= 1) {
348             return new int[] { 0 };
349         }
350 
351         // this is an implementation of definition 6 in Dan Kalman's paper.
352         final int vSize = valueCompiler.lowerIndirection.length;
353         final int dSize = derivativeCompiler.lowerIndirection.length;
354         final int[] lowerIndirection = new int[vSize + dSize];
355         System.arraycopy(valueCompiler.lowerIndirection, 0, lowerIndirection, 0, vSize);
356         for (int i = 0; i < dSize; ++i) {
357             lowerIndirection[vSize + i] = valueCompiler.getSize() + derivativeCompiler.lowerIndirection[i];
358         }
359 
360         return lowerIndirection;
361 
362     }
363 
364     /** Compile the multiplication indirection array.
365      * <p>
366      * This indirection array contains the indices of all pairs of elements
367      * involved when computing a multiplication. This allows a straightforward
368      * loop-based multiplication (see {@link #multiply(double[], int, double[], int, double[], int)}).
369      * </p>
370      * @param parameters number of free parameters
371      * @param order derivation order
372      * @param valueCompiler compiler for the value part
373      * @param derivativeCompiler compiler for the derivative part
374      * @param lowerIndirection lower derivatives indirection array
375      * @return multiplication indirection array
376      */
377     private static MultiplicationMapper[][] compileMultiplicationIndirection(final int parameters, final int order,
378                                                                              final DSCompiler valueCompiler,
379                                                                              final DSCompiler derivativeCompiler,
380                                                                              final int[] lowerIndirection) {
381 
382         if (parameters == 0 || order == 0) {
383             return new MultiplicationMapper[][] { { new MultiplicationMapper(1, 0, 0) } };
384         }
385 
386         // this is an implementation of definition 3 in Dan Kalman's paper.
387         final int vSize = valueCompiler.multIndirection.length;
388         final int dSize = derivativeCompiler.multIndirection.length;
389         final MultiplicationMapper[][] multIndirection = new MultiplicationMapper[vSize + dSize][];
390 
391         System.arraycopy(valueCompiler.multIndirection, 0, multIndirection, 0, vSize);
392 
393         for (int i = 0; i < dSize; ++i) {
394             final MultiplicationMapper[] dRow = derivativeCompiler.multIndirection[i];
395             final List<MultiplicationMapper> row = new ArrayList<>(dRow.length * 2);
396             for (MultiplicationMapper dj : dRow) {
397                 row.add(new MultiplicationMapper(dj.getCoeff(), lowerIndirection[dj.lhsIndex], vSize + dj.rhsIndex));
398                 row.add(new MultiplicationMapper(dj.getCoeff(), vSize + dj.lhsIndex, lowerIndirection[dj.rhsIndex]));
399             }
400             multIndirection[vSize + i] = combineSimilarTerms(row);
401         }
402 
403         return multIndirection;
404 
405     }
406 
407     /** Compile the function composition indirection array.
408      * <p>
409      * This indirection array contains the indices of all sets of elements
410      * involved when computing a composition. This allows a straightforward
411      * loop-based composition (see {@link #compose(double[], int, double[], double[], int)}).
412      * </p>
413      * @param parameters number of free parameters
414      * @param order derivation order
415      * @param valueCompiler compiler for the value part
416      * @param derivativeCompiler compiler for the derivative part
417      * @param sizes sizes array
418      * @param derivativesIndirection derivatives indirection array
419      * @return multiplication indirection array
420      * @throws MathIllegalArgumentException if order is too large
421      */
422     private static UnivariateCompositionMapper[][] compileCompositionIndirection(final int parameters, final int order,
423                                                                                  final DSCompiler valueCompiler,
424                                                                                  final DSCompiler derivativeCompiler,
425                                                                                  final int[][] sizes,
426                                                                                  final int[][] derivativesIndirection)
427        throws MathIllegalArgumentException {
428 
429         if (parameters == 0 || order == 0) {
430             return new UnivariateCompositionMapper[][] { { new UnivariateCompositionMapper(1, 0, new int[0]) } };
431         }
432 
433         final int vSize = valueCompiler.compIndirection.length;
434         final int dSize = derivativeCompiler.compIndirection.length;
435         final UnivariateCompositionMapper[][] compIndirection = new UnivariateCompositionMapper[vSize + dSize][];
436 
437         // the composition rules from the value part can be reused as is
438         System.arraycopy(valueCompiler.compIndirection, 0, compIndirection, 0, vSize);
439 
440         // the composition rules for the derivative part are deduced by
441         // differentiating the rules from the underlying compiler once
442         // with respect to the parameter this compiler handles and the
443         // underlying one did not handle
444         for (int i = 0; i < dSize; ++i) {
445             List<UnivariateCompositionMapper> row = new ArrayList<>();
446             for (UnivariateCompositionMapper term : derivativeCompiler.compIndirection[i]) {
447 
448                 // handle term p * f_k(g(x)) * g_l1(x) * g_l2(x) * ... * g_lp(x)
449 
450                 // derive the first factor in the term: f_k with respect to new parameter
451                 UnivariateCompositionMapper derivedTermF = new UnivariateCompositionMapper(term.getCoeff(),  // p
452                                                                                            term.fIndex + 1,  // f_(k+1)
453                                                                                            new int[term.dsIndices.length + 1]);
454                 int[] orders = new int[parameters];
455                 orders[parameters - 1] = 1;
456                 derivedTermF.dsIndices[term.dsIndices.length] = getPartialDerivativeIndex(parameters, order, sizes, orders);  // g_1
457                 for (int j = 0; j < term.dsIndices.length; ++j) {
458                     // convert the indices as the mapping for the current order
459                     // is different from the mapping with one less order
460                     derivedTermF.dsIndices[j] = convertIndex(term.dsIndices[j], parameters,
461                                                            derivativeCompiler.derivativesOrders,
462                                                            parameters, order, sizes);
463                 }
464                 derivedTermF.sort();
465                 row.add(derivedTermF);
466 
467                 // derive the various g_l
468                 for (int l = 0; l < term.dsIndices.length; ++l) {
469                     UnivariateCompositionMapper derivedTermG = new UnivariateCompositionMapper(term.getCoeff(),
470                                                                                                term.fIndex,
471                                                                                                new int[term.dsIndices.length]);
472                     for (int j = 0; j < term.dsIndices.length; ++j) {
473                         // convert the indices as the mapping for the current order
474                         // is different from the mapping with one less order
475                         derivedTermG.dsIndices[j] = convertIndex(term.dsIndices[j], parameters,
476                                                                derivativeCompiler.derivativesOrders,
477                                                                parameters, order, sizes);
478                         if (j == l) {
479                             // derive this term
480                             System.arraycopy(derivativesIndirection[derivedTermG.dsIndices[j]], 0, orders, 0, parameters);
481                             orders[parameters - 1]++;
482                             derivedTermG.dsIndices[j] = getPartialDerivativeIndex(parameters, order, sizes, orders);
483                         }
484                     }
485                     derivedTermG.sort();
486                     row.add(derivedTermG);
487                 }
488 
489             }
490 
491             // combine terms with similar derivation orders
492             compIndirection[vSize + i] = combineSimilarTerms(row);
493 
494         }
495 
496         return compIndirection;
497 
498     }
499 
500     /** Get rebaser, creating it if needed.
501      * @param baseCompiler compiler associated with the low level parameter functions
502      * @return rebaser for the number of base variables specified
503      * @since 2.2
504      */
505     private MultivariateCompositionMapper[][] getRebaser(final DSCompiler baseCompiler) {
506         synchronized (rebaseIndirection) {
507 
508             final int m = baseCompiler.getFreeParameters();
509             while (rebaseIndirection.size() <= m) {
510                 rebaseIndirection.add(null);
511             }
512 
513             if (rebaseIndirection.get(m) == null) {
514                 // we need to create the rebaser from instance to m base variables
515 
516                 if (order == 0) {
517                     // at order 0, the rebaser just copies the function value
518                     final MultivariateCompositionMapper[][] rebaser  = {
519                         { new MultivariateCompositionMapper(1, 0, new int[0]) }
520                     };
521                     rebaseIndirection.set(m, rebaser);
522                     return rebaser;
523                 }
524 
525                 // at order n > 0, the rebaser starts from the rebaser at order n-1
526                 // so the first rows of the rebaser (corresponding to orders 0 to n-1)
527                 // are just copies of the lower rebaser rows with indices adjusted,
528                 // the last row corresponding to order n is a term ∂ⁿf/∂qⱼ⋯∂qₖ∂qₗ
529                 // which can be written ∂(∂fⁿ⁻¹/∂qⱼ⋯∂qₖ)/∂qₗ, selecting any arbitrary
530                 // qₗ with non-zero derivation order as the base for recursion
531                 // the lower level rebaser provides ∂fⁿ⁻¹/∂qⱼ⋯∂qₖ as a
532                 // sum of products: Σᵢ ∂fⁿ⁻¹/∂pᵤ⋯∂pᵥ ∂pᵤ/∂qⱼ⋯∂qₖ ⋯ ∂pᵥ/∂qⱼ⋯∂qₖ
533                 // so we have to differentiate this sum of products
534                 //   - the term ∂fⁿ⁻¹/∂pᵤ⋯∂pᵥ depends on the p intermediate variables,
535                 //     not on the q base variables, so we use the composition formula
536                 //     ∂g/∂qₗ = Σᵢ ∂g/∂pᵢ ∂pᵢ/∂qₗ
537                 //   - the terms ∂pᵤ/∂qⱼ⋯∂qₖ are directly the intermediate variables p and we
538                 //     know their derivatives with respect to the base variables q
539                 final int baseSize = baseCompiler.getSize();
540                 final MultivariateCompositionMapper[][] rebaser = initializeFromLowerRebaser(baseCompiler);
541 
542                 // derivatives for last order
543                 for (int k = 1; k < baseSize; ++k) {
544                     // outer loop on rebased derivatives
545                     // at each iteration of the loop we are dealing with one derivative
546                     // like for example ∂³f/∂qⱼ∂qₖ∂qₗ, i.e. the components the rebaser produces
547                     if (rebaser[k] == null) {
548                         // the entry has not been set earlier
549                         // it is an entry of the form ∂ⁿf/∂qⱼ⋯∂qₖ∂qₗ where n is max order
550                         final List<MultivariateCompositionMapper> row = new ArrayList<>();
551 
552                         // find a variable with respect to which we have a derivative
553                         final int[] orders = baseCompiler.derivativesOrders[k].clone();
554                         int qIndex = -1;
555                         for (int j = 0; j < orders.length; ++j) {
556                             if (orders[j] > 0) {
557                                 qIndex = j;
558                                 break;
559                             }
560                         }
561 
562                         // find the entry corresponding to differentiating one order less with respect to this variable
563                         // ∂fⁿ⁻¹/∂qⱼ⋯∂qₖ
564                         orders[qIndex]--;
565                         final MultivariateCompositionMapper[] lowerRow =
566                                         rebaser[baseCompiler.getPartialDerivativeIndex(orders)];
567 
568                         // apply recursion formula
569                         for (final MultivariateCompositionMapper lowerTerm : lowerRow) {
570 
571                             for (int i = 0; i < parameters; ++i) {
572                                 // differentiate the term ∂fⁿ⁻¹/∂pᵤ⋯∂pᵥ part
573                                 row.add(differentiateFPart(lowerTerm, i, qIndex, baseCompiler));
574                             }
575 
576                             // differentiate the products ∂pᵤ/∂qⱼ⋯∂qₖ ⋯ ∂pᵥ/∂qⱼ⋯∂qₖ
577                             for (int j = 0; j < lowerTerm.productIndices.length; ++j) {
578                                 row.add(differentiateProductPart(lowerTerm, j, qIndex, baseCompiler));
579                             }
580 
581                         }
582 
583                         // simplifies and store the completed entry
584                         rebaser[k] = combineSimilarTerms(row);
585 
586                     }
587 
588                 }
589 
590                 rebaseIndirection.set(m, rebaser);
591 
592             }
593 
594             return rebaseIndirection.get(m);
595 
596         }
597     }
598 
599     /** Initialize a rebaser by copying the rules from a lower rebaser.
600      * @param baseCompiler compiler associated with the low level parameter functions
601      * @return rebaser with rules up to order - 1 copied (with indices adjusted)
602      * @since 2.2
603      */
604     private MultivariateCompositionMapper[][] initializeFromLowerRebaser(final DSCompiler baseCompiler) {
605 
606         // get the rebaser at order - 1
607         final DSCompiler lowerCompiler     = getCompiler(parameters, order - 1);
608         final DSCompiler lowerBaseCompiler = getCompiler(baseCompiler.parameters, order - 1);
609         final int        lowerBaseSize     = lowerBaseCompiler.getSize();
610         final MultivariateCompositionMapper[][] lowerRebaser = lowerCompiler.getRebaser(lowerBaseCompiler);
611 
612         // allocate array for rebaser at current order
613         final int baseSize = baseCompiler.getSize();
614         final MultivariateCompositionMapper[][] rebaser = new MultivariateCompositionMapper[baseSize][];
615 
616         // copy the rebasing rules for orders 0 to order - 1, adjusting indices
617         for (int i = 0; i < lowerRebaser.length; ++i) {
618             final int index = convertIndex(i, lowerBaseCompiler.parameters, lowerBaseCompiler.derivativesOrders,
619                                            baseCompiler.parameters, baseCompiler.order, baseCompiler.sizes);
620             rebaser[index] = new MultivariateCompositionMapper[lowerRebaser[i].length];
621             for (int j = 0; j < rebaser[index].length; ++j) {
622                 final int coeff  = lowerRebaser[i][j].getCoeff();
623                 final int dsIndex = convertIndex(lowerRebaser[i][j].dsIndex,
624                                                  lowerCompiler.parameters, lowerCompiler.derivativesOrders,
625                                                  parameters, order, sizes);
626                 final int[] productIndices = new int[lowerRebaser[i][j].productIndices.length];
627                 for (int k = 0; k < productIndices.length; ++k) {
628                     final int pIndex      = lowerRebaser[i][j].productIndices[k] / lowerBaseSize;
629                     final int baseDSIndex = lowerRebaser[i][j].productIndices[k] % lowerBaseSize;
630                     productIndices[k] = pIndex * baseSize +
631                                         convertIndex(baseDSIndex,
632                                                      lowerBaseCompiler.parameters, lowerBaseCompiler.derivativesOrders,
633                                                      baseCompiler.parameters, baseCompiler.order, baseCompiler.sizes);
634                 }
635                 rebaser[index][j] = new MultivariateCompositionMapper(coeff, dsIndex, productIndices);
636             }
637         }
638 
639         return rebaser;
640 
641     }
642 
643     /** Differentiate the ∂fⁿ⁻¹/∂pᵤ⋯∂pᵥ part of a {@link MultivariateCompositionMapper}.
644      * @param lowerTerm term to differentiate
645      * @param i index of the intermediate variable pᵢ
646      * @param qIndex index of the qₗ variable
647      * @param baseCompiler compiler associated with the low level parameter functions
648      * @return ∂fⁿ⁻¹/∂pᵤ⋯∂pᵥ
649      */
650     private MultivariateCompositionMapper differentiateFPart(final MultivariateCompositionMapper lowerTerm,
651                                                              final int i, final int qIndex,
652                                                              final DSCompiler baseCompiler) {
653 
654         // differentiate the term ∂fⁿ⁻¹/∂pᵤ⋯∂pᵥ with respect to pi
655         final int[] termOrders = derivativesOrders[lowerTerm.dsIndex].clone();
656         termOrders[i]++;
657 
658         // multiply by ∂pᵢ/∂qₗ
659         final int fDSIndex = getPartialDerivativeIndex(termOrders);
660         final int[] productIndicesF = new int[lowerTerm.productIndices.length + 1];
661         System.arraycopy(lowerTerm.productIndices, 0, productIndicesF, 0, lowerTerm.productIndices.length);
662         final int[] qOrders = new int[baseCompiler.parameters];
663         qOrders[qIndex] = 1;
664         productIndicesF[productIndicesF.length - 1] = i * baseCompiler.getSize() +
665                                                       baseCompiler.getPartialDerivativeIndex(qOrders);
666 
667         // generate the differentiated term
668         final MultivariateCompositionMapper termF =
669                         new MultivariateCompositionMapper(lowerTerm.getCoeff(), fDSIndex, productIndicesF);
670         termF.sort();
671         return termF;
672 
673     }
674 
675     /** Differentiate a product part of a {@link MultivariateCompositionMapper}.
676      * @param lowerTerm term to differentiate
677      * @param j index of the product to differentiate
678      * @param qIndex index of the qₗ variable
679      * @param baseCompiler compiler associated with the low level parameter functions
680      * @return ∂fⁿ⁻¹/∂pᵤ⋯∂pᵥ
681      */
682     private MultivariateCompositionMapper differentiateProductPart(final MultivariateCompositionMapper lowerTerm,
683                                                                    final int j, final int qIndex,
684                                                                    final DSCompiler baseCompiler) {
685 
686         // get derivation orders of ∂p/∂q
687         final int baseSize              = baseCompiler.getSize();
688         final int[] productIndicesP     = lowerTerm.productIndices.clone();
689         final int   pIndex              = productIndicesP[j] / baseSize;
690         final int   pDSIndex            = productIndicesP[j] % baseSize;
691         final int[] pOrders             = baseCompiler.getPartialDerivativeOrders(pDSIndex);
692 
693         // derive once more with respect to the selected q
694         pOrders[qIndex]++;
695         final int   pDSIndexHigherOrder = baseCompiler.getPartialDerivativeIndex(pOrders);
696         productIndicesP[j]              = pIndex * baseSize + pDSIndexHigherOrder;
697 
698         // create new term
699         final MultivariateCompositionMapper termP =
700                         new MultivariateCompositionMapper(lowerTerm.getCoeff(), lowerTerm.dsIndex, productIndicesP);
701         termP.sort();
702         return termP;
703 
704     }
705 
706     /** Get the index of a partial derivative in the array.
707      * <p>
708      * If all orders are set to 0, then the 0<sup>th</sup> order derivative
709      * is returned, which is the value of the function.
710      * </p>
711      * <p>The indices of derivatives are between 0 and {@link #getSize() getSize()} - 1.
712      * Their specific order is fixed for a given compiler, but otherwise not
713      * publicly specified. There are however some simple cases which have guaranteed
714      * indices:
715      * </p>
716      * <ul>
717      *   <li>the index of 0<sup>th</sup> order derivative is always 0</li>
718      *   <li>if there is only 1 {@link #getFreeParameters() free parameter}, then the
719      *   derivatives are sorted in increasing derivation order (i.e. f at index 0, df/dp
720      *   at index 1, d<sup>2</sup>f/dp<sup>2</sup> at index 2 ...
721      *   d<sup>k</sup>f/dp<sup>k</sup> at index k),</li>
722      *   <li>if the {@link #getOrder() derivation order} is 1, then the derivatives
723      *   are sorted in increasing free parameter order (i.e. f at index 0, df/dx<sub>1</sub>
724      *   at index 1, df/dx<sub>2</sub> at index 2 ... df/dx<sub>k</sub> at index k),</li>
725      *   <li>all other cases are not publicly specified</li>
726      * </ul>
727      * <p>
728      * This method is the inverse of method {@link #getPartialDerivativeOrders(int)}
729      * </p>
730      * @param orders derivation orders with respect to each parameter
731      * @return index of the partial derivative
732      * @exception MathIllegalArgumentException if the numbers of parameters does not
733      * match the instance
734      * @exception MathIllegalArgumentException if sum of derivation orders is larger
735      * than the instance limits
736      * @see #getPartialDerivativeOrders(int)
737      */
738     public int getPartialDerivativeIndex(final int ... orders)
739             throws MathIllegalArgumentException {
740 
741         // safety check
742         MathUtils.checkDimension(orders.length, getFreeParameters());
743         return getPartialDerivativeIndex(parameters, order, sizes, orders);
744 
745     }
746 
747     /** Get the index of a partial derivative in an array.
748      * @param parameters number of free parameters
749      * @param order derivation order
750      * @param sizes sizes array
751      * @param orders derivation orders with respect to each parameter
752      * (the length of this array must match the number of parameters)
753      * @return index of the partial derivative
754      * @exception MathIllegalArgumentException if sum of derivation orders is larger
755      * than the instance limits
756      */
757     private static int getPartialDerivativeIndex(final int parameters, final int order,
758                                                  final int[][] sizes, final int ... orders)
759         throws MathIllegalArgumentException {
760 
761         // the value is obtained by diving into the recursive Dan Kalman's structure
762         // this is theorem 2 of his paper, with recursion replaced by iteration
763         int index     = 0;
764         int m         = order;
765         int ordersSum = 0;
766         for (int i = parameters - 1; i >= 0; --i) {
767 
768             // derivative order for current free parameter
769             int derivativeOrder = orders[i];
770 
771             // safety check
772             ordersSum += derivativeOrder;
773             if (ordersSum > order) {
774                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
775                                                        ordersSum, order);
776             }
777 
778             while (derivativeOrder > 0) {
779                 --derivativeOrder;
780                 // as long as we differentiate according to current free parameter,
781                 // we have to skip the value part and dive into the derivative part
782                 // so we add the size of the value part to the base index
783                 index += sizes[i][m--];
784             }
785 
786         }
787 
788         return index;
789 
790     }
791 
792     /** Convert an index from one (parameters, order) structure to another.
793      * @param index index of a partial derivative in source derivative structure
794      * @param srcP number of free parameters in source derivative structure
795      * @param srcDerivativesOrders derivatives orders array for the source
796      * derivative structure
797      * @param destP number of free parameters in destination derivative structure
798      * @param destO derivation order in destination derivative structure
799      * @param destSizes sizes array for the destination derivative structure
800      * @return index of the partial derivative with the <em>same</em> characteristics
801      * in destination derivative structure
802      * @throws MathIllegalArgumentException if order is too large
803      */
804     private static int convertIndex(final int index,
805                                     final int srcP, final int[][] srcDerivativesOrders,
806                                     final int destP, final int destO, final int[][] destSizes)
807         throws MathIllegalArgumentException {
808         int[] orders = new int[destP];
809         System.arraycopy(srcDerivativesOrders[index], 0, orders, 0, FastMath.min(srcP, destP));
810         return getPartialDerivativeIndex(destP, destO, destSizes, orders);
811     }
812 
813     /** Get the derivation orders for a specific index in the array.
814      * <p>
815      * This method is the inverse of {@link #getPartialDerivativeIndex(int...)}.
816      * </p>
817      * @param index of the partial derivative
818      * @return derivation orders with respect to each parameter
819      * @see #getPartialDerivativeIndex(int...)
820      */
821     public int[] getPartialDerivativeOrders(final int index) {
822         return derivativesOrders[index].clone();
823     }
824 
825     /** Get the sum of derivation orders for a specific index in the array.
826      * <p>
827      * This method return the sum of the elements returned by
828      * {@link #getPartialDerivativeIndex(int...)}, using precomputed
829      * values
830      * </p>
831      * @param index of the partial derivative
832      * @return sum of derivation orders with respect to each parameter
833      * @see #getPartialDerivativeIndex(int...)
834      * @since 2.2
835      */
836     public int getPartialDerivativeOrdersSum(final int index) {
837         return derivativesOrdersSum[index];
838     }
839 
840     /** Get the number of free parameters.
841      * @return number of free parameters
842      */
843     public int getFreeParameters() {
844         return parameters;
845     }
846 
847     /** Get the derivation order.
848      * @return derivation order
849      */
850     public int getOrder() {
851         return order;
852     }
853 
854     /** Get the array size required for holding partial derivatives data.
855      * <p>
856      * This number includes the single 0 order derivative element, which is
857      * guaranteed to be stored in the first element of the array.
858      * </p>
859      * @return array size required for holding partial derivatives data
860      */
861     public int getSize() {
862         return sizes[parameters][order];
863     }
864 
865     /** Compute linear combination.
866      * The derivative structure built will be a1 * ds1 + a2 * ds2
867      * @param a1 first scale factor
868      * @param c1 first base (unscaled) component
869      * @param offset1 offset of first operand in its array
870      * @param a2 second scale factor
871      * @param c2 second base (unscaled) component
872      * @param offset2 offset of second operand in its array
873      * @param result array where result must be stored (it may be
874      * one of the input arrays)
875      * @param resultOffset offset of the result in its array
876      */
877     public void linearCombination(final double a1, final double[] c1, final int offset1,
878                                   final double a2, final double[] c2, final int offset2,
879                                   final double[] result, final int resultOffset) {
880         for (int i = 0; i < getSize(); ++i) {
881             result[resultOffset + i] =
882                     MathArrays.linearCombination(a1, c1[offset1 + i], a2, c2[offset2 + i]);
883         }
884     }
885 
886     /** Compute linear combination.
887      * The derivative structure built will be a1 * ds1 + a2 * ds2
888      * @param a1 first scale factor
889      * @param c1 first base (unscaled) component
890      * @param offset1 offset of first operand in its array
891      * @param a2 second scale factor
892      * @param c2 second base (unscaled) component
893      * @param offset2 offset of second operand in its array
894      * @param result array where result must be stored (it may be
895      * one of the input arrays)
896      * @param resultOffset offset of the result in its array
897      * @param <T> the type of the function parameters and value
898      */
899     public <T extends CalculusFieldElement<T>> void linearCombination(final T a1, final T[] c1, final int offset1,
900                                                                       final T a2, final T[] c2, final int offset2,
901                                                                       final T[] result, final int resultOffset) {
902         for (int i = 0; i < getSize(); ++i) {
903             result[resultOffset + i] =
904                     a1.linearCombination(a1, c1[offset1 + i], a2, c2[offset2 + i]);
905         }
906     }
907 
908     /** Compute linear combination.
909      * The derivative structure built will be a1 * ds1 + a2 * ds2
910      * @param a1 first scale factor
911      * @param c1 first base (unscaled) component
912      * @param offset1 offset of first operand in its array
913      * @param a2 second scale factor
914      * @param c2 second base (unscaled) component
915      * @param offset2 offset of second operand in its array
916      * @param result array where result must be stored (it may be
917      * one of the input arrays)
918      * @param resultOffset offset of the result in its array
919      * @param <T> the type of the function parameters and value
920      */
921     public <T extends CalculusFieldElement<T>> void linearCombination(final double a1, final T[] c1, final int offset1,
922                                                                       final double a2, final T[] c2, final int offset2,
923                                                                       final T[] result, final int resultOffset) {
924         for (int i = 0; i < getSize(); ++i) {
925             result[resultOffset + i] =
926                     c1[offset1].linearCombination(a1, c1[offset1 + i], a2, c2[offset2 + i]);
927         }
928     }
929 
930     /** Compute linear combination.
931      * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
932      * @param a1 first scale factor
933      * @param c1 first base (unscaled) component
934      * @param offset1 offset of first operand in its array
935      * @param a2 second scale factor
936      * @param c2 second base (unscaled) component
937      * @param offset2 offset of second operand in its array
938      * @param a3 third scale factor
939      * @param c3 third base (unscaled) component
940      * @param offset3 offset of third operand in its array
941      * @param result array where result must be stored (it may be
942      * one of the input arrays)
943      * @param resultOffset offset of the result in its array
944      */
945     public void linearCombination(final double a1, final double[] c1, final int offset1,
946                                   final double a2, final double[] c2, final int offset2,
947                                   final double a3, final double[] c3, final int offset3,
948                                   final double[] result, final int resultOffset) {
949         for (int i = 0; i < getSize(); ++i) {
950             result[resultOffset + i] =
951                     MathArrays.linearCombination(a1, c1[offset1 + i],
952                                                  a2, c2[offset2 + i],
953                                                  a3, c3[offset3 + i]);
954         }
955     }
956 
957     /** Compute linear combination.
958      * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
959      * @param a1 first scale factor
960      * @param c1 first base (unscaled) component
961      * @param offset1 offset of first operand in its array
962      * @param a2 second scale factor
963      * @param c2 second base (unscaled) component
964      * @param offset2 offset of second operand in its array
965      * @param a3 third scale factor
966      * @param c3 third base (unscaled) component
967      * @param offset3 offset of third operand in its array
968      * @param result array where result must be stored (it may be
969      * one of the input arrays)
970      * @param resultOffset offset of the result in its array
971      * @param <T> the type of the function parameters and value
972      */
973     public <T extends CalculusFieldElement<T>> void linearCombination(final T a1, final T[] c1, final int offset1,
974                                                                       final T a2, final T[] c2, final int offset2,
975                                                                       final T a3, final T[] c3, final int offset3,
976                                                                       final T[] result, final int resultOffset) {
977         for (int i = 0; i < getSize(); ++i) {
978             result[resultOffset + i] =
979                     a1.linearCombination(a1, c1[offset1 + i],
980                                          a2, c2[offset2 + i],
981                                          a3, c3[offset3 + i]);
982         }
983     }
984 
985     /** Compute linear combination.
986      * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
987      * @param a1 first scale factor
988      * @param c1 first base (unscaled) component
989      * @param offset1 offset of first operand in its array
990      * @param a2 second scale factor
991      * @param c2 second base (unscaled) component
992      * @param offset2 offset of second operand in its array
993      * @param a3 third scale factor
994      * @param c3 third base (unscaled) component
995      * @param offset3 offset of third operand in its array
996      * @param result array where result must be stored (it may be
997      * one of the input arrays)
998      * @param resultOffset offset of the result in its array
999      * @param <T> the type of the function parameters and value
1000      */
1001     public <T extends CalculusFieldElement<T>> void linearCombination(final double a1, final T[] c1, final int offset1,
1002                                                                       final double a2, final T[] c2, final int offset2,
1003                                                                       final double a3, final T[] c3, final int offset3,
1004                                                                       final T[] result, final int resultOffset) {
1005         for (int i = 0; i < getSize(); ++i) {
1006             result[resultOffset + i] =
1007                     c1[offset1].linearCombination(a1, c1[offset1 + i],
1008                                                   a2, c2[offset2 + i],
1009                                                   a3, c3[offset3 + i]);
1010         }
1011     }
1012 
1013     /** Compute linear combination.
1014      * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
1015      * @param a1 first scale factor
1016      * @param c1 first base (unscaled) component
1017      * @param offset1 offset of first operand in its array
1018      * @param a2 second scale factor
1019      * @param c2 second base (unscaled) component
1020      * @param offset2 offset of second operand in its array
1021      * @param a3 third scale factor
1022      * @param c3 third base (unscaled) component
1023      * @param offset3 offset of third operand in its array
1024      * @param a4 fourth scale factor
1025      * @param c4 fourth base (unscaled) component
1026      * @param offset4 offset of fourth operand in its array
1027      * @param result array where result must be stored (it may be
1028      * one of the input arrays)
1029      * @param resultOffset offset of the result in its array
1030      */
1031     public void linearCombination(final double a1, final double[] c1, final int offset1,
1032                                   final double a2, final double[] c2, final int offset2,
1033                                   final double a3, final double[] c3, final int offset3,
1034                                   final double a4, final double[] c4, final int offset4,
1035                                   final double[] result, final int resultOffset) {
1036         for (int i = 0; i < getSize(); ++i) {
1037             result[resultOffset + i] =
1038                     MathArrays.linearCombination(a1, c1[offset1 + i],
1039                                                  a2, c2[offset2 + i],
1040                                                  a3, c3[offset3 + i],
1041                                                  a4, c4[offset4 + i]);
1042         }
1043     }
1044 
1045     /** Compute linear combination.
1046      * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
1047      * @param a1 first scale factor
1048      * @param c1 first base (unscaled) component
1049      * @param offset1 offset of first operand in its array
1050      * @param a2 second scale factor
1051      * @param c2 second base (unscaled) component
1052      * @param offset2 offset of second operand in its array
1053      * @param a3 third scale factor
1054      * @param c3 third base (unscaled) component
1055      * @param offset3 offset of third operand in its array
1056      * @param a4 fourth scale factor
1057      * @param c4 fourth base (unscaled) component
1058      * @param offset4 offset of fourth operand in its array
1059      * @param result array where result must be stored (it may be
1060      * one of the input arrays)
1061      * @param resultOffset offset of the result in its array
1062      * @param <T> the type of the function parameters and value
1063      */
1064     public <T extends CalculusFieldElement<T>> void linearCombination(final T a1, final T[] c1, final int offset1,
1065                                                                       final T a2, final T[] c2, final int offset2,
1066                                                                       final T a3, final T[] c3, final int offset3,
1067                                                                       final T a4, final T[] c4, final int offset4,
1068                                                                       final T[] result, final int resultOffset) {
1069         for (int i = 0; i < getSize(); ++i) {
1070             result[resultOffset + i] =
1071                     a1.linearCombination(a1, c1[offset1 + i],
1072                                          a2, c2[offset2 + i],
1073                                          a3, c3[offset3 + i],
1074                                          a4, c4[offset4 + i]);
1075         }
1076     }
1077 
1078     /** Compute linear combination.
1079      * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
1080      * @param a1 first scale factor
1081      * @param c1 first base (unscaled) component
1082      * @param offset1 offset of first operand in its array
1083      * @param a2 second scale factor
1084      * @param c2 second base (unscaled) component
1085      * @param offset2 offset of second operand in its array
1086      * @param a3 third scale factor
1087      * @param c3 third base (unscaled) component
1088      * @param offset3 offset of third operand in its array
1089      * @param a4 fourth scale factor
1090      * @param c4 fourth base (unscaled) component
1091      * @param offset4 offset of fourth operand in its array
1092      * @param result array where result must be stored (it may be
1093      * one of the input arrays)
1094      * @param resultOffset offset of the result in its array
1095      * @param <T> the type of the function parameters and value
1096      */
1097     public <T extends CalculusFieldElement<T>> void linearCombination(final double a1, final T[] c1, final int offset1,
1098                                                                       final double a2, final T[] c2, final int offset2,
1099                                                                       final double a3, final T[] c3, final int offset3,
1100                                                                       final double a4, final T[] c4, final int offset4,
1101                                                                       final T[] result, final int resultOffset) {
1102         for (int i = 0; i < getSize(); ++i) {
1103             result[resultOffset + i] =
1104                             c1[offset1].linearCombination(a1, c1[offset1 + i],
1105                                                           a2, c2[offset2 + i],
1106                                                           a3, c3[offset3 + i],
1107                                                           a4, c4[offset4 + i]);
1108         }
1109     }
1110 
1111     /** Perform addition of two derivative structures.
1112      * @param lhs array holding left hand side of addition
1113      * @param lhsOffset offset of the left hand side in its array
1114      * @param rhs array right hand side of addition
1115      * @param rhsOffset offset of the right hand side in its array
1116      * @param result array where result must be stored (it may be
1117      * one of the input arrays)
1118      * @param resultOffset offset of the result in its array
1119      */
1120     public void add(final double[] lhs, final int lhsOffset,
1121                     final double[] rhs, final int rhsOffset,
1122                     final double[] result, final int resultOffset) {
1123         for (int i = 0; i < getSize(); ++i) {
1124             result[resultOffset + i] = lhs[lhsOffset + i] + rhs[rhsOffset + i];
1125         }
1126     }
1127 
1128     /** Perform addition of two derivative structures.
1129      * @param lhs array holding left hand side of addition
1130      * @param lhsOffset offset of the left hand side in its array
1131      * @param rhs array right hand side of addition
1132      * @param rhsOffset offset of the right hand side in its array
1133      * @param result array where result must be stored (it may be
1134      * one of the input arrays)
1135      * @param resultOffset offset of the result in its array
1136      * @param <T> the type of the function parameters and value
1137      */
1138     public <T extends CalculusFieldElement<T>> void add(final T[] lhs, final int lhsOffset,
1139                                                         final T[] rhs, final int rhsOffset,
1140                                                         final T[] result, final int resultOffset) {
1141         for (int i = 0; i < getSize(); ++i) {
1142             result[resultOffset + i] = lhs[lhsOffset + i].add(rhs[rhsOffset + i]);
1143         }
1144     }
1145 
1146     /** Perform subtraction of two derivative structures.
1147      * @param lhs array holding left hand side of subtraction
1148      * @param lhsOffset offset of the left hand side in its array
1149      * @param rhs array right hand side of subtraction
1150      * @param rhsOffset offset of the right hand side in its array
1151      * @param result array where result must be stored (it may be
1152      * one of the input arrays)
1153      * @param resultOffset offset of the result in its array
1154      */
1155     public void subtract(final double[] lhs, final int lhsOffset,
1156                          final double[] rhs, final int rhsOffset,
1157                          final double[] result, final int resultOffset) {
1158         for (int i = 0; i < getSize(); ++i) {
1159             result[resultOffset + i] = lhs[lhsOffset + i] - rhs[rhsOffset + i];
1160         }
1161     }
1162 
1163     /** Perform subtraction of two derivative structures.
1164      * @param lhs array holding left hand side of subtraction
1165      * @param lhsOffset offset of the left hand side in its array
1166      * @param rhs array right hand side of subtraction
1167      * @param rhsOffset offset of the right hand side in its array
1168      * @param result array where result must be stored (it may be
1169      * one of the input arrays)
1170      * @param resultOffset offset of the result in its array
1171      * @param <T> the type of the function parameters and value
1172      */
1173     public <T extends CalculusFieldElement<T>> void subtract(final T[] lhs, final int lhsOffset,
1174                                                              final T[] rhs, final int rhsOffset,
1175                                                              final T[] result, final int resultOffset) {
1176         for (int i = 0; i < getSize(); ++i) {
1177             result[resultOffset + i] = lhs[lhsOffset + i].subtract(rhs[rhsOffset + i]);
1178         }
1179     }
1180 
1181    /** Perform multiplication of two derivative structures.
1182      * @param lhs array holding left hand side of multiplication
1183      * @param lhsOffset offset of the left hand side in its array
1184      * @param rhs array right hand side of multiplication
1185      * @param rhsOffset offset of the right hand side in its array
1186      * @param result array where result must be stored (for
1187      * multiplication the result array <em>cannot</em> be one of
1188      * the input arrays)
1189      * @param resultOffset offset of the result in its array
1190      */
1191     public void multiply(final double[] lhs, final int lhsOffset,
1192                          final double[] rhs, final int rhsOffset,
1193                          final double[] result, final int resultOffset) {
1194         for (int i = 0; i < multIndirection.length; ++i) {
1195             double r = 0;
1196             for (final MultiplicationMapper mapping : multIndirection[i]) {
1197                 r += mapping.getCoeff() *
1198                      lhs[lhsOffset + mapping.lhsIndex] *
1199                      rhs[rhsOffset + mapping.rhsIndex];
1200             }
1201             result[resultOffset + i] = r;
1202         }
1203     }
1204 
1205     /** Perform multiplication of two derivative structures.
1206      * @param lhs array holding left hand side of multiplication
1207      * @param lhsOffset offset of the left hand side in its array
1208      * @param rhs array right hand side of multiplication
1209      * @param rhsOffset offset of the right hand side in its array
1210      * @param result array where result must be stored (for
1211      * multiplication the result array <em>cannot</em> be one of
1212      * the input arrays)
1213      * @param resultOffset offset of the result in its array
1214      * @param <T> the type of the function parameters and value
1215      */
1216     public <T extends CalculusFieldElement<T>> void multiply(final T[] lhs, final int lhsOffset,
1217                                                              final T[] rhs, final int rhsOffset,
1218                                                              final T[] result, final int resultOffset) {
1219         T zero = lhs[lhsOffset].getField().getZero();
1220         for (int i = 0; i < multIndirection.length; ++i) {
1221             T r = zero;
1222             for (final MultiplicationMapper mapping : multIndirection[i]) {
1223                 r = r.add(lhs[lhsOffset + mapping.lhsIndex].
1224                           multiply(rhs[rhsOffset + mapping.rhsIndex]).
1225                           multiply(mapping.getCoeff()));
1226             }
1227             result[resultOffset + i] = r;
1228         }
1229     }
1230 
1231     /** Perform division of two derivative structures. Based on the multiplication operator.
1232      * @param lhs array holding left hand side of division
1233      * @param lhsOffset offset of the left hand side in its array
1234      * @param rhs array right hand side of division
1235      * @param rhsOffset offset of the right hand side in its array
1236      * @param result array where result must be stored (for
1237      * division the result array <em>cannot</em> be one of
1238      * the input arrays)
1239      * @param resultOffset offset of the result in its array
1240      */
1241     public void divide(final double[] lhs, final int lhsOffset,
1242                        final double[] rhs, final int rhsOffset,
1243                        final double[] result, final int resultOffset) {
1244         result[resultOffset] = lhs[lhsOffset] / rhs[rhsOffset];
1245         for (int i = 1; i < multIndirection.length; ++i) {
1246             result[resultOffset + i] = lhs[lhsOffset + i];
1247             for (int j = 0; j < multIndirection[i].length - 1; j++) {
1248                 final MultiplicationMapper mapping = multIndirection[i][j];
1249                 result[resultOffset + i] -= mapping.getCoeff() *
1250                         (result[resultOffset + mapping.lhsIndex] * rhs[rhsOffset + mapping.rhsIndex]);
1251             }
1252             result[resultOffset + i] /= rhs[lhsOffset] * multIndirection[i][0].getCoeff();
1253         }
1254     }
1255 
1256     /** Perform division of two derivative structures. Based on the multiplication operator.
1257      * @param lhs array holding left hand side of division
1258      * @param lhsOffset offset of the left hand side in its array
1259      * @param rhs array right hand side of division
1260      * @param rhsOffset offset of the right hand side in its array
1261      * @param result array where result must be stored (for
1262      * division the result array <em>cannot</em> be one of
1263      * the input arrays)
1264      * @param resultOffset offset of the result in its array
1265      * @param <T> the type of the function parameters and value
1266      */
1267     public <T extends CalculusFieldElement<T>> void divide(final T[] lhs, final int lhsOffset,
1268                                                            final T[] rhs, final int rhsOffset,
1269                                                            final T[] result, final int resultOffset) {
1270         final T zero = lhs[lhsOffset].getField().getZero();
1271         result[resultOffset] = lhs[lhsOffset].divide(rhs[rhsOffset]);
1272         for (int i = 1; i < multIndirection.length; ++i) {
1273             result[resultOffset + i] = lhs[lhsOffset + i].add(zero);
1274             for (int j = 0; j < multIndirection[i].length - 1; j++) {
1275                 final MultiplicationMapper mapping = multIndirection[i][j];
1276                 result[resultOffset + i] = result[resultOffset + i].subtract(
1277                         result[resultOffset + mapping.lhsIndex].multiply(rhs[rhsOffset + mapping.rhsIndex]).
1278                                 multiply(mapping.getCoeff()));
1279             }
1280             result[resultOffset + i] = result[resultOffset + i].divide(rhs[lhsOffset].
1281                     multiply(multIndirection[i][0].getCoeff()));
1282         }
1283     }
1284 
1285     /** Compute reciprocal of derivative structure. Based on the multiplication operator.
1286      * @param operand array holding the operand
1287      * @param operandOffset offset of the operand in its array
1288      * @param result array where result must be stored
1289      * @param resultOffset offset of the result in its array
1290      */
1291     public void reciprocal(final double[] operand, final int operandOffset,
1292                            final double[] result, final int resultOffset) {
1293         result[resultOffset] = 1. / operand[operandOffset];
1294         for (int i = 1; i < multIndirection.length; ++i) {
1295             result[resultOffset + i] = 0.;
1296             for (int j = 0; j < multIndirection[i].length - 1; j++) {
1297                 final MultiplicationMapper mapping = multIndirection[i][j];
1298                 result[resultOffset + i] -= mapping.getCoeff() *
1299                         (result[resultOffset + mapping.lhsIndex] * operand[operandOffset + mapping.rhsIndex]);
1300             }
1301             result[resultOffset + i] /= operand[operandOffset] * multIndirection[i][0].getCoeff();
1302         }
1303     }
1304 
1305     /** Compute reciprocal of derivative structure. Based on the multiplication operator.
1306      * @param operand array holding the operand
1307      * @param operandOffset offset of the operand in its array
1308      * @param result array where result must be stored
1309      * @param resultOffset offset of the result in its array
1310      * @param <T> the type of the function parameters and value
1311      */
1312     public <T extends CalculusFieldElement<T>> void reciprocal(final T[] operand, final int operandOffset,
1313                                                                final T[] result, final int resultOffset) {
1314         final T zero = operand[operandOffset].getField().getZero();
1315         result[resultOffset] = operand[operandOffset].reciprocal();
1316         for (int i = 1; i < multIndirection.length; ++i) {
1317             result[resultOffset + i] = zero;
1318             for (int j = 0; j < multIndirection[i].length - 1; j++) {
1319                 final MultiplicationMapper mapping = multIndirection[i][j];
1320                 result[resultOffset + i] = result[resultOffset + i].subtract(
1321                         (result[resultOffset + mapping.lhsIndex].multiply(operand[operandOffset + mapping.rhsIndex])).
1322                                 multiply(mapping.getCoeff()));
1323             }
1324             result[resultOffset + i] = result[resultOffset + i].divide(operand[operandOffset].
1325                     multiply(multIndirection[i][0].getCoeff()));
1326         }
1327     }
1328 
1329     /** Perform remainder of two derivative structures.
1330      * @param lhs array holding left hand side of remainder
1331      * @param lhsOffset offset of the left hand side in its array
1332      * @param rhs array right hand side of remainder
1333      * @param rhsOffset offset of the right hand side in its array
1334      * @param result array where result must be stored (it may be
1335      * one of the input arrays)
1336      * @param resultOffset offset of the result in its array
1337      */
1338     public void remainder(final double[] lhs, final int lhsOffset,
1339                           final double[] rhs, final int rhsOffset,
1340                           final double[] result, final int resultOffset) {
1341 
1342         // compute k such that lhs % rhs = lhs - k rhs
1343         final double rem = FastMath.IEEEremainder(lhs[lhsOffset], rhs[rhsOffset]);
1344         final double k   = FastMath.rint((lhs[lhsOffset] - rem) / rhs[rhsOffset]);
1345 
1346         // set up value
1347         result[resultOffset] = rem;
1348 
1349         // set up partial derivatives
1350         for (int i = 1; i < getSize(); ++i) {
1351             result[resultOffset + i] = lhs[lhsOffset + i] - k * rhs[rhsOffset + i];
1352         }
1353 
1354     }
1355 
1356     /** Perform remainder of two derivative structures.
1357      * @param lhs array holding left hand side of remainder
1358      * @param lhsOffset offset of the left hand side in its array
1359      * @param rhs array right hand side of remainder
1360      * @param rhsOffset offset of the right hand side in its array
1361      * @param result array where result must be stored (it may be
1362      * one of the input arrays)
1363      * @param resultOffset offset of the result in its array
1364      * @param <T> the type of the function parameters and value
1365      */
1366     public <T extends CalculusFieldElement<T>> void remainder(final T[] lhs, final int lhsOffset,
1367                                                               final T[] rhs, final int rhsOffset,
1368                                                               final T[] result, final int resultOffset) {
1369 
1370         // compute k such that lhs % rhs = lhs - k rhs
1371         final T      rem = lhs[lhsOffset].remainder(rhs[rhsOffset]);
1372         final double k   = FastMath.rint((lhs[lhsOffset].getReal() - rem.getReal()) / rhs[rhsOffset].getReal());
1373 
1374         // set up value
1375         result[resultOffset] = rem;
1376 
1377         // set up partial derivatives
1378         for (int i = 1; i < getSize(); ++i) {
1379             result[resultOffset + i] = lhs[lhsOffset + i].subtract(rhs[rhsOffset + i].multiply(k));
1380         }
1381 
1382     }
1383 
1384     /** Compute power of a double to a derivative structure.
1385      * @param a number to exponentiate
1386      * @param operand array holding the power
1387      * @param operandOffset offset of the power in its array
1388      * @param result array where result must be stored (for
1389      * power the result array <em>cannot</em> be the input
1390      * array)
1391      * @param resultOffset offset of the result in its array
1392      */
1393     public void pow(final double a,
1394                     final double[] operand, final int operandOffset,
1395                     final double[] result, final int resultOffset) {
1396 
1397         // create the function value and derivatives
1398         // [a^x, ln(a) a^x, ln(a)^2 a^x,, ln(a)^3 a^x, ... ]
1399         final double[] function = new double[1 + order];
1400         if (a == 0) {
1401             if (operand[operandOffset] == 0) {
1402                 function[0] = 1;
1403                 double infinity = Double.POSITIVE_INFINITY;
1404                 for (int i = 1; i < function.length; ++i) {
1405                     infinity = -infinity;
1406                     function[i] = infinity;
1407                 }
1408             } else if (operand[operandOffset] < 0) {
1409                 Arrays.fill(function, Double.NaN);
1410             }
1411         } else {
1412             function[0] = FastMath.pow(a, operand[operandOffset]);
1413             final double lnA = FastMath.log(a);
1414             for (int i = 1; i < function.length; ++i) {
1415                 function[i] = lnA * function[i - 1];
1416             }
1417         }
1418 
1419 
1420         // apply function composition
1421         compose(operand, operandOffset, function, result, resultOffset);
1422 
1423     }
1424 
1425     /** Compute power of a double to a derivative structure.
1426      * @param a number to exponentiate
1427      * @param operand array holding the power
1428      * @param operandOffset offset of the power in its array
1429      * @param result array where result must be stored (for
1430      * power the result array <em>cannot</em> be the input
1431      * array)
1432      * @param resultOffset offset of the result in its array
1433      * @param <T> the type of the function parameters and value
1434      */
1435     public <T extends CalculusFieldElement<T>> void pow(final double a,
1436                                                         final T[] operand, final int operandOffset,
1437                                                         final T[] result, final int resultOffset) {
1438 
1439         final T zero = operand[operandOffset].getField().getZero();
1440 
1441         // create the function value and derivatives
1442         // [a^x, ln(a) a^x, ln(a)^2 a^x,, ln(a)^3 a^x, ... ]
1443         final T[] function = MathArrays.buildArray(operand[operandOffset].getField(), 1 + order);
1444         if (a == 0) {
1445             if (operand[operandOffset].getReal() == 0) {
1446                 function[0] = zero.add(1);
1447                 T infinity = zero.add(Double.POSITIVE_INFINITY);
1448                 for (int i = 1; i < function.length; ++i) {
1449                     infinity = infinity.negate();
1450                     function[i] = infinity;
1451                 }
1452             } else if (operand[operandOffset].getReal() < 0) {
1453                 Arrays.fill(function, zero.add(Double.NaN));
1454             }
1455         } else {
1456             function[0] = zero.add(a).pow(operand[operandOffset]);
1457             final double lnA = FastMath.log(a);
1458             for (int i = 1; i < function.length; ++i) {
1459                 function[i] = function[i - 1].multiply(lnA);
1460             }
1461         }
1462 
1463 
1464         // apply function composition
1465         compose(operand, operandOffset, function, result, resultOffset);
1466 
1467     }
1468 
1469     /** Compute power of a derivative structure.
1470      * @param operand array holding the operand
1471      * @param operandOffset offset of the operand in its array
1472      * @param p power to apply
1473      * @param result array where result must be stored (for
1474      * power the result array <em>cannot</em> be the input
1475      * array)
1476      * @param resultOffset offset of the result in its array
1477      */
1478     public void pow(final double[] operand, final int operandOffset, final double p,
1479                     final double[] result, final int resultOffset) {
1480 
1481         if (p == 0) {
1482             // special case, x^0 = 1 for all x
1483             result[resultOffset] = 1.0;
1484             Arrays.fill(result, resultOffset + 1, resultOffset + getSize(), 0);
1485             return;
1486         }
1487 
1488         if (operand[operandOffset] == 0) {
1489             // special case, 0^p = 0 for all p
1490             Arrays.fill(result, resultOffset, resultOffset + getSize(), 0);
1491             return;
1492         }
1493 
1494         // create the function value and derivatives
1495         // [x^p, px^(p-1), p(p-1)x^(p-2), ... ]
1496         double[] function = new double[1 + order];
1497         double xk = FastMath.pow(operand[operandOffset], p - order);
1498         for (int i = order; i > 0; --i) {
1499             function[i] = xk;
1500             xk *= operand[operandOffset];
1501         }
1502         function[0] = xk;
1503         double coefficient = p;
1504         for (int i = 1; i <= order; ++i) {
1505             function[i] *= coefficient;
1506             coefficient *= p - i;
1507         }
1508 
1509         // apply function composition
1510         compose(operand, operandOffset, function, result, resultOffset);
1511 
1512     }
1513 
1514     /** Compute power of a derivative structure.
1515      * @param operand array holding the operand
1516      * @param operandOffset offset of the operand in its array
1517      * @param p power to apply
1518      * @param result array where result must be stored (for
1519      * power the result array <em>cannot</em> be the input
1520      * array)
1521      * @param resultOffset offset of the result in its array
1522      * @param <T> the type of the function parameters and value
1523      */
1524     public <T extends CalculusFieldElement<T>> void pow(final T[] operand, final int operandOffset, final double p,
1525                                                         final T[] result, final int resultOffset) {
1526 
1527         final Field<T> field = operand[operandOffset].getField();
1528 
1529         if (p == 0) {
1530             // special case, x^0 = 1 for all x
1531             result[resultOffset] = field.getOne();
1532             Arrays.fill(result, resultOffset + 1, resultOffset + getSize(), field.getZero());
1533             return;
1534         }
1535 
1536         if (operand[operandOffset].getReal() == 0) {
1537             // special case, 0^p = 0 for all p
1538             Arrays.fill(result, resultOffset, resultOffset + getSize(), field.getZero());
1539             return;
1540         }
1541 
1542         // create the function value and derivatives
1543         // [x^p, px^(p-1), p(p-1)x^(p-2), ... ]
1544         T[] function = MathArrays.buildArray(field, 1 + order);
1545         T xk = operand[operandOffset].pow(p - order);
1546         for (int i = order; i > 0; --i) {
1547             function[i] = xk;
1548             xk = xk.multiply(operand[operandOffset]);
1549         }
1550         function[0] = xk;
1551         double coefficient = p;
1552         for (int i = 1; i <= order; ++i) {
1553             function[i]  = function[i].multiply(coefficient);
1554             coefficient *= p - i;
1555         }
1556 
1557         // apply function composition
1558         compose(operand, operandOffset, function, result, resultOffset);
1559 
1560     }
1561 
1562     /** Compute integer power of a derivative structure.
1563      * @param operand array holding the operand
1564      * @param operandOffset offset of the operand in its array
1565      * @param n power to apply
1566      * @param result array where result must be stored (for
1567      * power the result array <em>cannot</em> be the input
1568      * array)
1569      * @param resultOffset offset of the result in its array
1570      */
1571     public void pow(final double[] operand, final int operandOffset, final int n,
1572                     final double[] result, final int resultOffset) {
1573 
1574         if (n == 0) {
1575             // special case, x^0 = 1 for all x
1576             result[resultOffset] = 1.0;
1577             Arrays.fill(result, resultOffset + 1, resultOffset + getSize(), 0);
1578             return;
1579         }
1580 
1581         // create the power function value and derivatives
1582         // [x^n, nx^(n-1), n(n-1)x^(n-2), ... ]
1583         double[] function = new double[1 + order];
1584 
1585         if (n > 0) {
1586             // strictly positive power
1587             final int maxOrder = FastMath.min(order, n);
1588             double xk = FastMath.pow(operand[operandOffset], n - maxOrder);
1589             for (int i = maxOrder; i > 0; --i) {
1590                 function[i] = xk;
1591                 xk *= operand[operandOffset];
1592             }
1593             function[0] = xk;
1594         } else {
1595             // strictly negative power
1596             final double inv = 1.0 / operand[operandOffset];
1597             double xk = FastMath.pow(inv, -n);
1598             for (int i = 0; i <= order; ++i) {
1599                 function[i] = xk;
1600                 xk *= inv;
1601             }
1602         }
1603 
1604         double coefficient = n;
1605         for (int i = 1; i <= order; ++i) {
1606             function[i] *= coefficient;
1607             coefficient *= n - i;
1608         }
1609 
1610         // apply function composition
1611         compose(operand, operandOffset, function, result, resultOffset);
1612 
1613     }
1614 
1615     /** Compute integer power of a derivative structure.
1616      * @param operand array holding the operand
1617      * @param operandOffset offset of the operand in its array
1618      * @param n power to apply
1619      * @param result array where result must be stored (for
1620      * power the result array <em>cannot</em> be the input
1621      * array)
1622      * @param resultOffset offset of the result in its array
1623      * @param <T> the type of the function parameters and value
1624      */
1625     public <T extends CalculusFieldElement<T>> void pow(final T[] operand, final int operandOffset, final int n,
1626                                                         final T[] result, final int resultOffset) {
1627 
1628         final Field<T> field = operand[operandOffset].getField();
1629 
1630         if (n == 0) {
1631             // special case, x^0 = 1 for all x
1632             result[resultOffset] = field.getOne();
1633             Arrays.fill(result, resultOffset + 1, resultOffset + getSize(), field.getZero());
1634             return;
1635         }
1636 
1637         // create the power function value and derivatives
1638         // [x^n, nx^(n-1), n(n-1)x^(n-2), ... ]
1639         T[] function = MathArrays.buildArray(field, 1 + order);
1640 
1641         if (n > 0) {
1642             // strictly positive power
1643             final int maxOrder = FastMath.min(order, n);
1644             T xk = operand[operandOffset].pow(n - maxOrder);
1645             for (int i = maxOrder; i > 0; --i) {
1646                 function[i] = xk;
1647                 xk = xk.multiply(operand[operandOffset]);
1648             }
1649             function[0] = xk;
1650         } else {
1651             // strictly negative power
1652             final T inv = operand[operandOffset].reciprocal();
1653             T xk = inv.pow(-n);
1654             for (int i = 0; i <= order; ++i) {
1655                 function[i] = xk;
1656                 xk = xk.multiply(inv);
1657             }
1658         }
1659 
1660         double coefficient = n;
1661         for (int i = 1; i <= order; ++i) {
1662             function[i]  = function[i].multiply(coefficient);
1663             coefficient *= n - i;
1664         }
1665 
1666         // apply function composition
1667         compose(operand, operandOffset, function, result, resultOffset);
1668 
1669     }
1670 
1671     /** Compute power of a derivative structure.
1672      * @param x array holding the base
1673      * @param xOffset offset of the base in its array
1674      * @param y array holding the exponent
1675      * @param yOffset offset of the exponent in its array
1676      * @param result array where result must be stored (for
1677      * power the result array <em>cannot</em> be the input
1678      * array)
1679      * @param resultOffset offset of the result in its array
1680      */
1681     public void pow(final double[] x, final int xOffset,
1682                     final double[] y, final int yOffset,
1683                     final double[] result, final int resultOffset) {
1684         final double[] logX = new double[getSize()];
1685         log(x, xOffset, logX, 0);
1686         final double[] yLogX = new double[getSize()];
1687         multiply(logX, 0, y, yOffset, yLogX, 0);
1688         exp(yLogX, 0, result, resultOffset);
1689     }
1690 
1691     /** Compute power of a derivative structure.
1692      * @param x array holding the base
1693      * @param xOffset offset of the base in its array
1694      * @param y array holding the exponent
1695      * @param yOffset offset of the exponent in its array
1696      * @param result array where result must be stored (for
1697      * power the result array <em>cannot</em> be the input
1698      * array)
1699      * @param resultOffset offset of the result in its array
1700      * @param <T> the type of the function parameters and value
1701      */
1702     public <T extends CalculusFieldElement<T>> void pow(final T[] x, final int xOffset,
1703                                                         final T[] y, final int yOffset,
1704                                                         final T[] result, final int resultOffset) {
1705         final T[] logX = MathArrays.buildArray(x[xOffset].getField(), getSize());
1706         log(x, xOffset, logX, 0);
1707         final T[] yLogX = MathArrays.buildArray(x[xOffset].getField(), getSize());
1708         multiply(logX, 0, y, yOffset, yLogX, 0);
1709         exp(yLogX, 0, result, resultOffset);
1710     }
1711 
1712     /** Compute square root of a derivative structure. Based on the multiplication operator.
1713      * @param operand array holding the operand
1714      * @param operandOffset offset of the operand in its array
1715      * @param result array where result must be stored (for
1716      * square root the result array <em>cannot</em> be the input
1717      * array)
1718      * @param resultOffset offset of the result in its array
1719      */
1720     public void sqrt(final double[] operand, final int operandOffset,
1721                      final double[] result, final int resultOffset) {
1722         final double sqrtConstant = FastMath.sqrt(operand[operandOffset]);
1723         result[resultOffset] = sqrtConstant;
1724         for (int i = 1; i < multIndirection.length; ++i) {
1725             result[resultOffset + i] = operand[operandOffset + i];
1726             for (int j = 1; j < multIndirection[i].length - 1; j++) {
1727                 final MultiplicationMapper mapping = multIndirection[i][j];
1728                 result[resultOffset + i] -= mapping.getCoeff() *
1729                         (result[resultOffset + mapping.lhsIndex] * result[operandOffset + mapping.rhsIndex]);
1730             }
1731             result[resultOffset + i] /= sqrtConstant * (multIndirection[i][multIndirection[i].length - 1].getCoeff() +
1732                     multIndirection[i][0].getCoeff());
1733         }
1734     }
1735 
1736     /** Compute square root of a derivative structure. Based on the multiplication operator.
1737      * @param operand array holding the operand
1738      * @param operandOffset offset of the operand in its array
1739      * @param result array where result must be stored (for
1740      * square root the result array <em>cannot</em> be the input
1741      * array)
1742      * @param resultOffset offset of the result in its array
1743      * @param <T> the type of the function parameters and value
1744      */
1745     public <T extends CalculusFieldElement<T>> void sqrt(final T[] operand, final int operandOffset,
1746                                                          final T[] result, final int resultOffset) {
1747         final T zero = operand[operandOffset].getField().getZero();
1748         final T sqrtConstant = operand[operandOffset].sqrt();
1749         result[resultOffset] = sqrtConstant.add(zero);
1750         for (int i = 1; i < multIndirection.length; ++i) {
1751             result[resultOffset + i] = operand[operandOffset + i].add(zero);
1752             for (int j = 1; j < multIndirection[i].length - 1; j++) {
1753                 final MultiplicationMapper mapping = multIndirection[i][j];
1754                 result[resultOffset + i] = result[resultOffset + i].subtract(
1755                         (result[resultOffset + mapping.lhsIndex].multiply(result[operandOffset + mapping.rhsIndex])).
1756                                 multiply(mapping.getCoeff()));
1757             }
1758             result[resultOffset + i] = result[resultOffset + i].divide(sqrtConstant.multiply(
1759                     multIndirection[i][0].getCoeff() + multIndirection[i][multIndirection[i].length - 1].getCoeff()));
1760         }
1761     }
1762 
1763     /** Compute n<sup>th</sup> root of a derivative structure.
1764      * @param operand array holding the operand
1765      * @param operandOffset offset of the operand in its array
1766      * @param n order of the root
1767      * @param result array where result must be stored (for
1768      * n<sup>th</sup> root the result array <em>cannot</em> be the input
1769      * array)
1770      * @param resultOffset offset of the result in its array
1771      */
1772     public void rootN(final double[] operand, final int operandOffset, final int n,
1773                       final double[] result, final int resultOffset) {
1774 
1775         // create the function value and derivatives
1776         // [x^(1/n), (1/n)x^((1/n)-1), (1-n)/n^2x^((1/n)-2), ... ]
1777         double[] function = new double[1 + order];
1778         double xk;
1779         if (n == 2) {
1780             function[0] = FastMath.sqrt(operand[operandOffset]);
1781             xk          = 0.5 / function[0];
1782         } else if (n == 3) {
1783             function[0] = FastMath.cbrt(operand[operandOffset]);
1784             xk          = 1.0 / (3.0 * function[0] * function[0]);
1785         } else {
1786             function[0] = FastMath.pow(operand[operandOffset], 1.0 / n);
1787             xk          = 1.0 / (n * FastMath.pow(function[0], n - 1));
1788         }
1789         final double nReciprocal = 1.0 / n;
1790         final double xReciprocal = 1.0 / operand[operandOffset];
1791         for (int i = 1; i <= order; ++i) {
1792             function[i] = xk;
1793             xk *= xReciprocal * (nReciprocal - i);
1794         }
1795 
1796         // apply function composition
1797         compose(operand, operandOffset, function, result, resultOffset);
1798 
1799     }
1800 
1801     /** Compute n<sup>th</sup> root of a derivative structure.
1802      * @param operand array holding the operand
1803      * @param operandOffset offset of the operand in its array
1804      * @param n order of the root
1805      * @param result array where result must be stored (for
1806      * n<sup>th</sup> root the result array <em>cannot</em> be the input
1807      * array)
1808      * @param resultOffset offset of the result in its array
1809      * @param <T> the type of the function parameters and value
1810      */
1811     public <T extends CalculusFieldElement<T>> void rootN(final T[] operand, final int operandOffset, final int n,
1812                                                           final T[] result, final int resultOffset) {
1813 
1814         final Field<T> field = operand[operandOffset].getField();
1815 
1816         // create the function value and derivatives
1817         // [x^(1/n), (1/n)x^((1/n)-1), (1-n)/n^2x^((1/n)-2), ... ]
1818         T[] function = MathArrays.buildArray(field, 1 + order);
1819         T xk;
1820         if (n == 2) {
1821             function[0] = operand[operandOffset].sqrt();
1822             xk          = function[0].add(function[0]).reciprocal();
1823         } else if (n == 3) {
1824             function[0] = operand[operandOffset].cbrt();
1825             xk          = function[0].multiply(function[0]).multiply(3).reciprocal();
1826         } else {
1827             function[0] = operand[operandOffset].pow(1.0 / n);
1828             xk          = function[0].pow(n - 1).multiply(n).reciprocal();
1829         }
1830         final double nReciprocal = 1.0 / n;
1831         final T      xReciprocal = operand[operandOffset].reciprocal();
1832         for (int i = 1; i <= order; ++i) {
1833             function[i] = xk;
1834             xk = xk.multiply(xReciprocal.multiply(nReciprocal - i));
1835         }
1836 
1837         // apply function composition
1838         compose(operand, operandOffset, function, result, resultOffset);
1839 
1840     }
1841 
1842     /** Compute exponential of a derivative structure.
1843      * @param operand array holding the operand
1844      * @param operandOffset offset of the operand in its array
1845      * @param result array where result must be stored (for
1846      * exponential the result array <em>cannot</em> be the input
1847      * array)
1848      * @param resultOffset offset of the result in its array
1849      */
1850     public void exp(final double[] operand, final int operandOffset,
1851                     final double[] result, final int resultOffset) {
1852 
1853         // create the function value and derivatives
1854         double[] function = new double[1 + order];
1855         Arrays.fill(function, FastMath.exp(operand[operandOffset]));
1856 
1857         // apply function composition
1858         compose(operand, operandOffset, function, result, resultOffset);
1859 
1860     }
1861 
1862     /** Compute exponential of a derivative structure.
1863      * @param operand array holding the operand
1864      * @param operandOffset offset of the operand in its array
1865      * @param result array where result must be stored (for
1866      * exponential the result array <em>cannot</em> be the input
1867      * array)
1868      * @param resultOffset offset of the result in its array
1869      * @param <T> the type of the function parameters and value
1870      */
1871     public <T extends CalculusFieldElement<T>> void exp(final T[] operand, final int operandOffset,
1872                                                         final T[] result, final int resultOffset) {
1873 
1874         final Field<T> field = operand[operandOffset].getField();
1875 
1876         // create the function value and derivatives
1877         T[] function = MathArrays.buildArray(field, 1 + order);
1878         Arrays.fill(function, operand[operandOffset].exp());
1879 
1880         // apply function composition
1881         compose(operand, operandOffset, function, result, resultOffset);
1882 
1883     }
1884 
1885     /** Compute exp(x) - 1 of a derivative structure.
1886      * @param operand array holding the operand
1887      * @param operandOffset offset of the operand in its array
1888      * @param result array where result must be stored (for
1889      * exponential the result array <em>cannot</em> be the input
1890      * array)
1891      * @param resultOffset offset of the result in its array
1892      */
1893     public void expm1(final double[] operand, final int operandOffset,
1894                       final double[] result, final int resultOffset) {
1895 
1896         // create the function value and derivatives
1897         double[] function = new double[1 + order];
1898         function[0] = FastMath.expm1(operand[operandOffset]);
1899         Arrays.fill(function, 1, 1 + order, FastMath.exp(operand[operandOffset]));
1900 
1901         // apply function composition
1902         compose(operand, operandOffset, function, result, resultOffset);
1903 
1904     }
1905 
1906     /** Compute exp(x) - 1 of a derivative structure.
1907      * @param operand array holding the operand
1908      * @param operandOffset offset of the operand in its array
1909      * @param result array where result must be stored (for
1910      * exponential the result array <em>cannot</em> be the input
1911      * array)
1912      * @param resultOffset offset of the result in its array
1913      * @param <T> the type of the function parameters and value
1914      */
1915     public <T extends CalculusFieldElement<T>> void expm1(final T[] operand, final int operandOffset,
1916                                                           final T[] result, final int resultOffset) {
1917 
1918         final Field<T> field = operand[operandOffset].getField();
1919 
1920         // create the function value and derivatives
1921         T[] function = MathArrays.buildArray(field, 1 + order);
1922         function[0] = operand[operandOffset].expm1();
1923         Arrays.fill(function, 1, 1 + order, operand[operandOffset].exp());
1924 
1925         // apply function composition
1926         compose(operand, operandOffset, function, result, resultOffset);
1927 
1928     }
1929 
1930     /** Compute natural logarithm of a derivative structure.
1931      * @param operand array holding the operand
1932      * @param operandOffset offset of the operand in its array
1933      * @param result array where result must be stored (for
1934      * logarithm the result array <em>cannot</em> be the input
1935      * array)
1936      * @param resultOffset offset of the result in its array
1937      */
1938     public void log(final double[] operand, final int operandOffset,
1939                     final double[] result, final int resultOffset) {
1940 
1941         // create the function value and derivatives
1942         double[] function = new double[1 + order];
1943         function[0] = FastMath.log(operand[operandOffset]);
1944         if (order > 0) {
1945             double inv = 1.0 / operand[operandOffset];
1946             double xk  = inv;
1947             for (int i = 1; i <= order; ++i) {
1948                 function[i] = xk;
1949                 xk *= -i * inv;
1950             }
1951         }
1952 
1953         // apply function composition
1954         compose(operand, operandOffset, function, result, resultOffset);
1955 
1956     }
1957 
1958     /** Compute natural logarithm of a derivative structure.
1959      * @param operand array holding the operand
1960      * @param operandOffset offset of the operand in its array
1961      * @param result array where result must be stored (for
1962      * logarithm the result array <em>cannot</em> be the input
1963      * array)
1964      * @param resultOffset offset of the result in its array
1965      * @param <T> the type of the function parameters and value
1966      */
1967     public <T extends CalculusFieldElement<T>> void log(final T[] operand, final int operandOffset,
1968                                                         final T[] result, final int resultOffset) {
1969 
1970         final Field<T> field = operand[operandOffset].getField();
1971 
1972         // create the function value and derivatives
1973         T[] function = MathArrays.buildArray(field, 1 + order);
1974         function[0] = operand[operandOffset].log();
1975         if (order > 0) {
1976             T inv = operand[operandOffset].reciprocal();
1977             T xk  = inv;
1978             for (int i = 1; i <= order; ++i) {
1979                 function[i] = xk;
1980                 xk = xk.multiply(inv.multiply(-i));
1981             }
1982         }
1983 
1984         // apply function composition
1985         compose(operand, operandOffset, function, result, resultOffset);
1986 
1987     }
1988 
1989     /** Computes shifted logarithm of a derivative structure.
1990      * @param operand array holding the operand
1991      * @param operandOffset offset of the operand in its array
1992      * @param result array where result must be stored (for
1993      * shifted logarithm the result array <em>cannot</em> be the input array)
1994      * @param resultOffset offset of the result in its array
1995      */
1996     public void log1p(final double[] operand, final int operandOffset,
1997                       final double[] result, final int resultOffset) {
1998 
1999         // create the function value and derivatives
2000         double[] function = new double[1 + order];
2001         function[0] = FastMath.log1p(operand[operandOffset]);
2002         if (order > 0) {
2003             double inv = 1.0 / (1.0 + operand[operandOffset]);
2004             double xk  = inv;
2005             for (int i = 1; i <= order; ++i) {
2006                 function[i] = xk;
2007                 xk *= -i * inv;
2008             }
2009         }
2010 
2011         // apply function composition
2012         compose(operand, operandOffset, function, result, resultOffset);
2013 
2014     }
2015 
2016     /** Computes shifted logarithm of a derivative structure.
2017      * @param operand array holding the operand
2018      * @param operandOffset offset of the operand in its array
2019      * @param result array where result must be stored (for
2020      * shifted logarithm the result array <em>cannot</em> be the input array)
2021      * @param resultOffset offset of the result in its array
2022      * @param <T> the type of the function parameters and value
2023      */
2024     public <T extends CalculusFieldElement<T>> void log1p(final T[] operand, final int operandOffset,
2025                                                           final T[] result, final int resultOffset) {
2026 
2027         final Field<T> field = operand[operandOffset].getField();
2028 
2029         // create the function value and derivatives
2030         T[] function = MathArrays.buildArray(field, 1 + order);
2031         function[0] = operand[operandOffset].log1p();
2032         if (order > 0) {
2033             T inv = operand[operandOffset].add(1).reciprocal();
2034             T xk  = inv;
2035             for (int i = 1; i <= order; ++i) {
2036                 function[i] = xk;
2037                 xk = xk.multiply(inv.multiply(-i));
2038             }
2039         }
2040 
2041         // apply function composition
2042         compose(operand, operandOffset, function, result, resultOffset);
2043 
2044     }
2045 
2046     /** Computes base 10 logarithm of a derivative structure.
2047      * @param operand array holding the operand
2048      * @param operandOffset offset of the operand in its array
2049      * @param result array where result must be stored (for
2050      * base 10 logarithm the result array <em>cannot</em> be the input array)
2051      * @param resultOffset offset of the result in its array
2052      */
2053     public void log10(final double[] operand, final int operandOffset,
2054                       final double[] result, final int resultOffset) {
2055 
2056         // create the function value and derivatives
2057         double[] function = new double[1 + order];
2058         function[0] = FastMath.log10(operand[operandOffset]);
2059         if (order > 0) {
2060             double inv = 1.0 / operand[operandOffset];
2061             double xk  = inv / FastMath.log(10.0);
2062             for (int i = 1; i <= order; ++i) {
2063                 function[i] = xk;
2064                 xk *= -i * inv;
2065             }
2066         }
2067 
2068         // apply function composition
2069         compose(operand, operandOffset, function, result, resultOffset);
2070 
2071     }
2072 
2073     /** Computes base 10 logarithm of a derivative structure.
2074      * @param operand array holding the operand
2075      * @param operandOffset offset of the operand in its array
2076      * @param result array where result must be stored (for
2077      * base 10 logarithm the result array <em>cannot</em> be the input array)
2078      * @param resultOffset offset of the result in its array
2079      * @param <T> the type of the function parameters and value
2080      */
2081     public <T extends CalculusFieldElement<T>> void log10(final T[] operand, final int operandOffset,
2082                                                           final T[] result, final int resultOffset) {
2083 
2084         final Field<T> field = operand[operandOffset].getField();
2085 
2086         // create the function value and derivatives
2087         T[] function = MathArrays.buildArray(field, 1 + order);
2088         function[0] = operand[operandOffset].log10();
2089         if (order > 0) {
2090             T inv = operand[operandOffset].reciprocal();
2091             T xk  = inv.multiply(1.0 / FastMath.log(10.0));
2092             for (int i = 1; i <= order; ++i) {
2093                 function[i] = xk;
2094                 xk = xk.multiply(inv.multiply(-i));
2095             }
2096         }
2097 
2098         // apply function composition
2099         compose(operand, operandOffset, function, result, resultOffset);
2100 
2101     }
2102 
2103     /** Compute cosine of a derivative structure.
2104      * @param operand array holding the operand
2105      * @param operandOffset offset of the operand in its array
2106      * @param result array where result must be stored (for
2107      * cosine the result array <em>cannot</em> be the input
2108      * array)
2109      * @param resultOffset offset of the result in its array
2110      */
2111     public void cos(final double[] operand, final int operandOffset,
2112                     final double[] result, final int resultOffset) {
2113 
2114         // create the function value and derivatives
2115         double[] function = new double[1 + order];
2116         final SinCos sinCos = FastMath.sinCos(operand[operandOffset]);
2117         function[0] = sinCos.cos();
2118         if (order > 0) {
2119             function[1] = -sinCos.sin();
2120             for (int i = 2; i <= order; ++i) {
2121                 function[i] = -function[i - 2];
2122             }
2123         }
2124 
2125         // apply function composition
2126         compose(operand, operandOffset, function, result, resultOffset);
2127 
2128     }
2129 
2130     /** Compute cosine of a derivative structure.
2131      * @param operand array holding the operand
2132      * @param operandOffset offset of the operand in its array
2133      * @param result array where result must be stored (for
2134      * cosine the result array <em>cannot</em> be the input
2135      * array)
2136      * @param resultOffset offset of the result in its array
2137      * @param <T> the type of the function parameters and value
2138      */
2139     public <T extends CalculusFieldElement<T>> void cos(final T[] operand, final int operandOffset,
2140                                                         final T[] result, final int resultOffset) {
2141 
2142         final Field<T> field = operand[operandOffset].getField();
2143 
2144         // create the function value and derivatives
2145         T[] function = MathArrays.buildArray(field, 1 + order);
2146         final FieldSinCos<T> sinCos = FastMath.sinCos(operand[operandOffset]);
2147         function[0] = sinCos.cos();
2148         if (order > 0) {
2149             function[1] = sinCos.sin().negate();
2150             if (order > 1) {
2151                 function[2] = sinCos.cos().negate();
2152                 if (order > 2) {
2153                     function[3] = sinCos.sin();
2154                     for (int i = 4; i <= order; ++i) {
2155                         function[i] = function[i - 4];
2156                     }
2157                 }
2158             }
2159         }
2160 
2161         // apply function composition
2162         compose(operand, operandOffset, function, result, resultOffset);
2163 
2164     }
2165 
2166     /** Compute sine of a derivative structure.
2167      * @param operand array holding the operand
2168      * @param operandOffset offset of the operand in its array
2169      * @param result array where result must be stored (for
2170      * sine the result array <em>cannot</em> be the input
2171      * array)
2172      * @param resultOffset offset of the result in its array
2173      */
2174     public void sin(final double[] operand, final int operandOffset,
2175                     final double[] result, final int resultOffset) {
2176 
2177         // create the function value and derivatives
2178         double[] function = new double[1 + order];
2179         final SinCos sinCos = FastMath.sinCos(operand[operandOffset]);
2180         function[0] = sinCos.sin();
2181         if (order > 0) {
2182             function[1] = sinCos.cos();
2183             for (int i = 2; i <= order; ++i) {
2184                 function[i] = -function[i - 2];
2185             }
2186         }
2187 
2188         // apply function composition
2189         compose(operand, operandOffset, function, result, resultOffset);
2190 
2191     }
2192 
2193     /** Compute sine of a derivative structure.
2194      * @param operand array holding the operand
2195      * @param operandOffset offset of the operand in its array
2196      * @param result array where result must be stored (for
2197      * sine the result array <em>cannot</em> be the input
2198      * array)
2199      * @param resultOffset offset of the result in its array
2200      * @param <T> the type of the function parameters and value
2201      */
2202     public <T extends CalculusFieldElement<T>> void sin(final T[] operand, final int operandOffset,
2203                                                         final T[] result, final int resultOffset) {
2204 
2205         final Field<T> field = operand[operandOffset].getField();
2206 
2207         // create the function value and derivatives
2208         T[] function = MathArrays.buildArray(field, 1 + order);
2209         final FieldSinCos<T> sinCos = FastMath.sinCos(operand[operandOffset]);
2210         function[0] = sinCos.sin();
2211         if (order > 0) {
2212             function[1] = sinCos.cos();
2213             if (order > 1) {
2214                 function[2] = sinCos.sin().negate();
2215                 if (order > 2) {
2216                     function[3] = sinCos.cos().negate();
2217                     for (int i = 4; i <= order; ++i) {
2218                         function[i] = function[i - 4];
2219                     }
2220                 }
2221             }
2222         }
2223 
2224         // apply function composition
2225         compose(operand, operandOffset, function, result, resultOffset);
2226 
2227     }
2228 
2229     /** Compute combined sine and cosine of a derivative structure.
2230      * @param operand array holding the operand
2231      * @param operandOffset offset of the operand in its array
2232      * @param sin array where sine must be stored (for
2233      * sine the result array <em>cannot</em> be the input
2234      * array)
2235      * @param sinOffset offset of the result in its array
2236      * @param cos array where cosine must be stored (for
2237      * cosine the result array <em>cannot</em> be the input
2238      * array)
2239      * @param cosOffset offset of the result in its array
2240      * @since 1.4
2241      */
2242     public void sinCos(final double[] operand, final int operandOffset,
2243                        final double[] sin, final int sinOffset,
2244                        final double[] cos, final int cosOffset) {
2245 
2246         // create the function value and derivatives
2247         double[] functionSin = new double[1 + order];
2248         double[] functionCos = new double[1 + order];
2249         final SinCos sinCos = FastMath.sinCos(operand[operandOffset]);
2250         functionSin[0] = sinCos.sin();
2251         functionCos[0] = sinCos.cos();
2252         if (order > 0) {
2253             functionSin[1] =  sinCos.cos();
2254             functionCos[1] = -sinCos.sin();
2255             for (int i = 2; i <= order; ++i) {
2256                 functionSin[i] = -functionSin[i - 2];
2257                 functionCos[i] = -functionCos[i - 2];
2258             }
2259         }
2260 
2261         // apply function composition
2262         compose(operand, operandOffset, functionSin, sin, sinOffset);
2263         compose(operand, operandOffset, functionCos, cos, cosOffset);
2264 
2265     }
2266 
2267     /** Compute combined sine and cosine of a derivative structure.
2268      * @param operand array holding the operand
2269      * @param operandOffset offset of the operand in its array
2270      * @param sin array where sine must be stored (for
2271      * sine the result array <em>cannot</em> be the input
2272      * array)
2273      * @param sinOffset offset of the result in its array
2274      * @param cos array where cosine must be stored (for
2275      * cosine the result array <em>cannot</em> be the input
2276      * array)
2277      * @param cosOffset offset of the result in its array
2278      * @param <T> the type of the function parameters and value
2279      * @since 1.4
2280      */
2281     public <T extends CalculusFieldElement<T>> void sinCos(final T[] operand, final int operandOffset,
2282                                                            final T[] sin, final int sinOffset,
2283                                                            final T[] cos, final int cosOffset) {
2284 
2285         final Field<T> field = operand[operandOffset].getField();
2286 
2287         // create the function value and derivatives
2288         T[] functionSin = MathArrays.buildArray(field, 1 + order);
2289         T[] functionCos = MathArrays.buildArray(field, 1 + order);
2290         final FieldSinCos<T> sinCos = FastMath.sinCos(operand[operandOffset]);
2291         functionCos[0] = sinCos.cos();
2292         if (order > 0) {
2293             functionCos[1] = sinCos.sin().negate();
2294             if (order > 1) {
2295                 functionCos[2] = sinCos.cos().negate();
2296                 if (order > 2) {
2297                     functionCos[3] = sinCos.sin();
2298                     for (int i = 4; i <= order; ++i) {
2299                         functionCos[i] = functionCos[i - 4];
2300                     }
2301                 }
2302             }
2303         }
2304         functionSin[0] = sinCos.sin();
2305         System.arraycopy(functionCos, 0, functionSin, 1, order);
2306 
2307         // apply function composition
2308         compose(operand, operandOffset, functionSin, sin, sinOffset);
2309         compose(operand, operandOffset, functionCos, cos, cosOffset);
2310 
2311     }
2312 
2313     /** Compute tangent of a derivative structure.
2314      * @param operand array holding the operand
2315      * @param operandOffset offset of the operand in its array
2316      * @param result array where result must be stored (for
2317      * tangent the result array <em>cannot</em> be the input
2318      * array)
2319      * @param resultOffset offset of the result in its array
2320      */
2321     public void tan(final double[] operand, final int operandOffset,
2322                     final double[] result, final int resultOffset) {
2323 
2324         // create the function value and derivatives
2325         final double[] function = new double[1 + order];
2326         final double t = FastMath.tan(operand[operandOffset]);
2327         function[0] = t;
2328 
2329         if (order > 0) {
2330 
2331             // the nth order derivative of tan has the form:
2332             // dn(tan(x)/dxn = P_n(tan(x))
2333             // where P_n(t) is a degree n+1 polynomial with same parity as n+1
2334             // P_0(t) = t, P_1(t) = 1 + t^2, P_2(t) = 2 t (1 + t^2) ...
2335             // the general recurrence relation for P_n is:
2336             // P_n(x) = (1+t^2) P_(n-1)'(t)
2337             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
2338             final double[] p = new double[order + 2];
2339             p[1] = 1;
2340             final double t2 = t * t;
2341             for (int n = 1; n <= order; ++n) {
2342 
2343                 // update and evaluate polynomial P_n(t)
2344                 double v = 0;
2345                 p[n + 1] = n * p[n];
2346                 for (int k = n + 1; k >= 0; k -= 2) {
2347                     v = v * t2 + p[k];
2348                     if (k > 2) {
2349                         p[k - 2] = (k - 1) * p[k - 1] + (k - 3) * p[k - 3];
2350                     } else if (k == 2) {
2351                         p[0] = p[1];
2352                     }
2353                 }
2354                 if ((n & 0x1) == 0) {
2355                     v *= t;
2356                 }
2357 
2358                 function[n] = v;
2359 
2360             }
2361         }
2362 
2363         // apply function composition
2364         compose(operand, operandOffset, function, result, resultOffset);
2365 
2366     }
2367 
2368     /** Compute tangent of a derivative structure.
2369      * @param operand array holding the operand
2370      * @param operandOffset offset of the operand in its array
2371      * @param result array where result must be stored (for
2372      * tangent the result array <em>cannot</em> be the input
2373      * array)
2374      * @param resultOffset offset of the result in its array
2375      * @param <T> the type of the function parameters and value
2376      */
2377     public <T extends CalculusFieldElement<T>> void tan(final T[] operand, final int operandOffset,
2378                                                         final T[] result, final int resultOffset) {
2379 
2380         final Field<T> field = operand[operandOffset].getField();
2381 
2382         // create the function value and derivatives
2383         T[] function = MathArrays.buildArray(field, 1 + order);
2384         final T t = operand[operandOffset].tan();
2385         function[0] = t;
2386 
2387         if (order > 0) {
2388 
2389             // the nth order derivative of tan has the form:
2390             // dn(tan(x)/dxn = P_n(tan(x))
2391             // where P_n(t) is a degree n+1 polynomial with same parity as n+1
2392             // P_0(t) = t, P_1(t) = 1 + t^2, P_2(t) = 2 t (1 + t^2) ...
2393             // the general recurrence relation for P_n is:
2394             // P_n(x) = (1+t^2) P_(n-1)'(t)
2395             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
2396             final T[] p = MathArrays.buildArray(field, order + 2);
2397             p[1] = field.getOne();
2398             final T t2 = t.multiply(t);
2399             for (int n = 1; n <= order; ++n) {
2400 
2401                 // update and evaluate polynomial P_n(t)
2402                 T v = field.getZero();
2403                 p[n + 1] = p[n].multiply(n);
2404                 for (int k = n + 1; k >= 0; k -= 2) {
2405                     v = v.multiply(t2).add(p[k]);
2406                     if (k > 2) {
2407                         p[k - 2] = p[k - 1].multiply(k - 1).add(p[k - 3].multiply(k - 3));
2408                     } else if (k == 2) {
2409                         p[0] = p[1];
2410                     }
2411                 }
2412                 if ((n & 0x1) == 0) {
2413                     v = v.multiply(t);
2414                 }
2415 
2416                 function[n] = v;
2417 
2418             }
2419         }
2420 
2421         // apply function composition
2422         compose(operand, operandOffset, function, result, resultOffset);
2423 
2424     }
2425 
2426     /** Compute arc cosine of a derivative structure.
2427      * @param operand array holding the operand
2428      * @param operandOffset offset of the operand in its array
2429      * @param result array where result must be stored (for
2430      * arc cosine the result array <em>cannot</em> be the input
2431      * array)
2432      * @param resultOffset offset of the result in its array
2433      */
2434     public void acos(final double[] operand, final int operandOffset,
2435                      final double[] result, final int resultOffset) {
2436 
2437         // create the function value and derivatives
2438         double[] function = new double[1 + order];
2439         final double x = operand[operandOffset];
2440         function[0] = FastMath.acos(x);
2441         if (order > 0) {
2442             // the nth order derivative of acos has the form:
2443             // dn(acos(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
2444             // where P_n(x) is a degree n-1 polynomial with same parity as n-1
2445             // P_1(x) = -1, P_2(x) = -x, P_3(x) = -2x^2 - 1 ...
2446             // the general recurrence relation for P_n is:
2447             // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
2448             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
2449             final double[] p = new double[order];
2450             p[0] = -1;
2451             final double x2    = x * x;
2452             final double f     = 1.0 / (1 - x2);
2453             double coeff = FastMath.sqrt(f);
2454             function[1] = coeff * p[0];
2455             for (int n = 2; n <= order; ++n) {
2456 
2457                 // update and evaluate polynomial P_n(x)
2458                 double v = 0;
2459                 p[n - 1] = (n - 1) * p[n - 2];
2460                 for (int k = n - 1; k >= 0; k -= 2) {
2461                     v = v * x2 + p[k];
2462                     if (k > 2) {
2463                         p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
2464                     } else if (k == 2) {
2465                         p[0] = p[1];
2466                     }
2467                 }
2468                 if ((n & 0x1) == 0) {
2469                     v *= x;
2470                 }
2471 
2472                 coeff *= f;
2473                 function[n] = coeff * v;
2474 
2475             }
2476         }
2477 
2478         // apply function composition
2479         compose(operand, operandOffset, function, result, resultOffset);
2480 
2481     }
2482 
2483     /** Compute arc cosine of a derivative structure.
2484      * @param operand array holding the operand
2485      * @param operandOffset offset of the operand in its array
2486      * @param result array where result must be stored (for
2487      * arc cosine the result array <em>cannot</em> be the input
2488      * array)
2489      * @param resultOffset offset of the result in its array
2490      * @param <T> the type of the function parameters and value
2491      */
2492     public <T extends CalculusFieldElement<T>> void acos(final T[] operand, final int operandOffset,
2493                                                          final T[] result, final int resultOffset) {
2494 
2495         final Field<T> field = operand[operandOffset].getField();
2496 
2497         // create the function value and derivatives
2498         T[] function = MathArrays.buildArray(field, 1 + order);
2499         final T x = operand[operandOffset];
2500         function[0] = x.acos();
2501         if (order > 0) {
2502             // the nth order derivative of acos has the form:
2503             // dn(acos(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
2504             // where P_n(x) is a degree n-1 polynomial with same parity as n-1
2505             // P_1(x) = -1, P_2(x) = -x, P_3(x) = -2x^2 - 1 ...
2506             // the general recurrence relation for P_n is:
2507             // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
2508             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
2509             final T[] p = MathArrays.buildArray(field, order);
2510             p[0] = field.getOne().negate();
2511             final T x2    = x.square();
2512             final T f     = x2.subtract(1).negate().reciprocal();
2513             T coeff = f.sqrt();
2514             function[1] = coeff.multiply(p[0]);
2515             for (int n = 2; n <= order; ++n) {
2516 
2517                 // update and evaluate polynomial P_n(x)
2518                 T v = field.getZero();
2519                 p[n - 1] = p[n - 2].multiply(n - 1);
2520                 for (int k = n - 1; k >= 0; k -= 2) {
2521                     v = v.multiply(x2).add(p[k]);
2522                     if (k > 2) {
2523                         p[k - 2] = p[k - 1].multiply(k - 1).add(p[k - 3].multiply(2 * n - k));
2524                     } else if (k == 2) {
2525                         p[0] = p[1];
2526                     }
2527                 }
2528                 if ((n & 0x1) == 0) {
2529                     v = v.multiply(x);
2530                 }
2531 
2532                 coeff = coeff.multiply(f);
2533                 function[n] = coeff.multiply(v);
2534 
2535             }
2536         }
2537 
2538         // apply function composition
2539         compose(operand, operandOffset, function, result, resultOffset);
2540 
2541     }
2542 
2543     /** Compute arc sine of a derivative structure.
2544      * @param operand array holding the operand
2545      * @param operandOffset offset of the operand in its array
2546      * @param result array where result must be stored (for
2547      * arc sine the result array <em>cannot</em> be the input
2548      * array)
2549      * @param resultOffset offset of the result in its array
2550      */
2551     public void asin(final double[] operand, final int operandOffset,
2552                      final double[] result, final int resultOffset) {
2553 
2554         // create the function value and derivatives
2555         double[] function = new double[1 + order];
2556         final double x = operand[operandOffset];
2557         function[0] = FastMath.asin(x);
2558         if (order > 0) {
2559             // the nth order derivative of asin has the form:
2560             // dn(asin(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
2561             // where P_n(x) is a degree n-1 polynomial with same parity as n-1
2562             // P_1(x) = 1, P_2(x) = x, P_3(x) = 2x^2 + 1 ...
2563             // the general recurrence relation for P_n is:
2564             // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
2565             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
2566             final double[] p = new double[order];
2567             p[0] = 1;
2568             final double x2    = x * x;
2569             final double f     = 1.0 / (1 - x2);
2570             double coeff = FastMath.sqrt(f);
2571             function[1] = coeff * p[0];
2572             for (int n = 2; n <= order; ++n) {
2573 
2574                 // update and evaluate polynomial P_n(x)
2575                 double v = 0;
2576                 p[n - 1] = (n - 1) * p[n - 2];
2577                 for (int k = n - 1; k >= 0; k -= 2) {
2578                     v = v * x2 + p[k];
2579                     if (k > 2) {
2580                         p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
2581                     } else if (k == 2) {
2582                         p[0] = p[1];
2583                     }
2584                 }
2585                 if ((n & 0x1) == 0) {
2586                     v *= x;
2587                 }
2588 
2589                 coeff *= f;
2590                 function[n] = coeff * v;
2591 
2592             }
2593         }
2594 
2595         // apply function composition
2596         compose(operand, operandOffset, function, result, resultOffset);
2597 
2598     }
2599 
2600     /** Compute arc sine of a derivative structure.
2601      * @param operand array holding the operand
2602      * @param operandOffset offset of the operand in its array
2603      * @param result array where result must be stored (for
2604      * arc sine the result array <em>cannot</em> be the input
2605      * array)
2606      * @param resultOffset offset of the result in its array
2607      * @param <T> the type of the function parameters and value
2608      */
2609     public <T extends CalculusFieldElement<T>> void asin(final T[] operand, final int operandOffset,
2610                                                          final T[] result, final int resultOffset) {
2611 
2612         final Field<T> field = operand[operandOffset].getField();
2613 
2614         // create the function value and derivatives
2615         T[] function = MathArrays.buildArray(field, 1 + order);
2616         final T x = operand[operandOffset];
2617         function[0] = x.asin();
2618         if (order > 0) {
2619             // the nth order derivative of asin has the form:
2620             // dn(asin(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
2621             // where P_n(x) is a degree n-1 polynomial with same parity as n-1
2622             // P_1(x) = 1, P_2(x) = x, P_3(x) = 2x^2 + 1 ...
2623             // the general recurrence relation for P_n is:
2624             // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
2625             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
2626             final T[] p = MathArrays.buildArray(field, order);
2627             p[0] = field.getOne();
2628             final T x2    = x.square();
2629             final T f     = x2.subtract(1).negate().reciprocal();
2630             T coeff = f.sqrt();
2631             function[1] = coeff.multiply(p[0]);
2632             for (int n = 2; n <= order; ++n) {
2633 
2634                 // update and evaluate polynomial P_n(x)
2635                 T v = field.getZero();
2636                 p[n - 1] = p[n - 2].multiply(n - 1);
2637                 for (int k = n - 1; k >= 0; k -= 2) {
2638                     v = v.multiply(x2).add(p[k]);
2639                     if (k > 2) {
2640                         p[k - 2] = p[k - 1].multiply(k - 1).add(p[k - 3].multiply(2 * n - k));
2641                     } else if (k == 2) {
2642                         p[0] = p[1];
2643                     }
2644                 }
2645                 if ((n & 0x1) == 0) {
2646                     v = v.multiply(x);
2647                 }
2648 
2649                 coeff = coeff.multiply(f);
2650                 function[n] = coeff.multiply(v);
2651 
2652             }
2653         }
2654 
2655         // apply function composition
2656         compose(operand, operandOffset, function, result, resultOffset);
2657 
2658     }
2659 
2660     /** Compute arc tangent of a derivative structure.
2661      * @param operand array holding the operand
2662      * @param operandOffset offset of the operand in its array
2663      * @param result array where result must be stored (for
2664      * arc tangent the result array <em>cannot</em> be the input
2665      * array)
2666      * @param resultOffset offset of the result in its array
2667      */
2668     public void atan(final double[] operand, final int operandOffset,
2669                      final double[] result, final int resultOffset) {
2670 
2671         // create the function value and derivatives
2672         double[] function = new double[1 + order];
2673         final double x = operand[operandOffset];
2674         function[0] = FastMath.atan(x);
2675         if (order > 0) {
2676             // the nth order derivative of atan has the form:
2677             // dn(atan(x)/dxn = Q_n(x) / (1 + x^2)^n
2678             // where Q_n(x) is a degree n-1 polynomial with same parity as n-1
2679             // Q_1(x) = 1, Q_2(x) = -2x, Q_3(x) = 6x^2 - 2 ...
2680             // the general recurrence relation for Q_n is:
2681             // Q_n(x) = (1+x^2) Q_(n-1)'(x) - 2(n-1) x Q_(n-1)(x)
2682             // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array
2683             final double[] q = new double[order];
2684             q[0] = 1;
2685             final double x2    = x * x;
2686             final double f     = 1.0 / (1 + x2);
2687             double coeff = f;
2688             function[1] = coeff * q[0];
2689             for (int n = 2; n <= order; ++n) {
2690 
2691                 // update and evaluate polynomial Q_n(x)
2692                 double v = 0;
2693                 q[n - 1] = -n * q[n - 2];
2694                 for (int k = n - 1; k >= 0; k -= 2) {
2695                     v = v * x2 + q[k];
2696                     if (k > 2) {
2697                         q[k - 2] = (k - 1) * q[k - 1] + (k - 1 - 2 * n) * q[k - 3];
2698                     } else if (k == 2) {
2699                         q[0] = q[1];
2700                     }
2701                 }
2702                 if ((n & 0x1) == 0) {
2703                     v *= x;
2704                 }
2705 
2706                 coeff *= f;
2707                 function[n] = coeff * v;
2708 
2709             }
2710         }
2711 
2712         // apply function composition
2713         compose(operand, operandOffset, function, result, resultOffset);
2714 
2715     }
2716 
2717     /** Compute arc tangent of a derivative structure.
2718      * @param operand array holding the operand
2719      * @param operandOffset offset of the operand in its array
2720      * @param result array where result must be stored (for
2721      * arc tangent the result array <em>cannot</em> be the input
2722      * array)
2723      * @param resultOffset offset of the result in its array
2724      * @param <T> the type of the function parameters and value
2725      */
2726     public <T extends CalculusFieldElement<T>> void atan(final T[] operand, final int operandOffset,
2727                                                          final T[] result, final int resultOffset) {
2728 
2729         final Field<T> field = operand[operandOffset].getField();
2730 
2731         // create the function value and derivatives
2732         T[] function = MathArrays.buildArray(field, 1 + order);
2733         final T x = operand[operandOffset];
2734         function[0] = x.atan();
2735         if (order > 0) {
2736             // the nth order derivative of atan has the form:
2737             // dn(atan(x)/dxn = Q_n(x) / (1 + x^2)^n
2738             // where Q_n(x) is a degree n-1 polynomial with same parity as n-1
2739             // Q_1(x) = 1, Q_2(x) = -2x, Q_3(x) = 6x^2 - 2 ...
2740             // the general recurrence relation for Q_n is:
2741             // Q_n(x) = (1+x^2) Q_(n-1)'(x) - 2(n-1) x Q_(n-1)(x)
2742             // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array
2743             final T[] q = MathArrays.buildArray(field, order);
2744             q[0] = field.getOne();
2745             final T x2    = x.square();
2746             final T f     = x2.add(1).reciprocal();
2747             T coeff = f;
2748             function[1] = coeff.multiply(q[0]);
2749             for (int n = 2; n <= order; ++n) {
2750 
2751                 // update and evaluate polynomial Q_n(x)
2752                 T v = field.getZero();
2753                 q[n - 1] = q[n - 2].multiply(-n);
2754                 for (int k = n - 1; k >= 0; k -= 2) {
2755                     v = v.multiply(x2).add(q[k]);
2756                     if (k > 2) {
2757                         q[k - 2] = q[k - 1].multiply(k - 1).add(q[k - 3].multiply(k - 1 - 2 * n));
2758                     } else if (k == 2) {
2759                         q[0] = q[1];
2760                     }
2761                 }
2762                 if ((n & 0x1) == 0) {
2763                     v = v.multiply(x);
2764                 }
2765 
2766                 coeff = coeff.multiply(f);
2767                 function[n] = coeff.multiply(v);
2768 
2769             }
2770         }
2771 
2772         // apply function composition
2773         compose(operand, operandOffset, function, result, resultOffset);
2774 
2775     }
2776 
2777     /** Compute two arguments arc tangent of a derivative structure.
2778      * @param y array holding the first operand
2779      * @param yOffset offset of the first operand in its array
2780      * @param x array holding the second operand
2781      * @param xOffset offset of the second operand in its array
2782      * @param result array where result must be stored (for
2783      * two arguments arc tangent the result array <em>cannot</em>
2784      * be the input array)
2785      * @param resultOffset offset of the result in its array
2786      */
2787     public void atan2(final double[] y, final int yOffset,
2788                       final double[] x, final int xOffset,
2789                       final double[] result, final int resultOffset) {
2790 
2791         // compute r = sqrt(x^2+y^2)
2792         double[] tmp1 = new double[getSize()];
2793         multiply(x, xOffset, x, xOffset, tmp1, 0);      // x^2
2794         double[] tmp2 = new double[getSize()];
2795         multiply(y, yOffset, y, yOffset, tmp2, 0);      // y^2
2796         add(tmp1, 0, tmp2, 0, tmp2, 0);                 // x^2 + y^2
2797         rootN(tmp2, 0, 2, tmp1, 0);                     // r = sqrt(x^2 + y^2)
2798 
2799         if (x[xOffset] >= 0) {
2800 
2801             // compute atan2(y, x) = 2 atan(y / (r + x))
2802             add(tmp1, 0, x, xOffset, tmp2, 0);          // r + x
2803             divide(y, yOffset, tmp2, 0, tmp1, 0);       // y /(r + x)
2804             atan(tmp1, 0, tmp2, 0);                     // atan(y / (r + x))
2805             for (int i = 0; i < tmp2.length; ++i) {
2806                 result[resultOffset + i] = 2 * tmp2[i]; // 2 * atan(y / (r + x))
2807             }
2808 
2809         } else {
2810 
2811             // compute atan2(y, x) = +/- pi - 2 atan(y / (r - x))
2812             subtract(tmp1, 0, x, xOffset, tmp2, 0);     // r - x
2813             divide(y, yOffset, tmp2, 0, tmp1, 0);       // y /(r - x)
2814             atan(tmp1, 0, tmp2, 0);                     // atan(y / (r - x))
2815             result[resultOffset] =
2816                     ((tmp2[0] <= 0) ? -FastMath.PI : FastMath.PI) - 2 * tmp2[0]; // +/-pi - 2 * atan(y / (r - x))
2817             for (int i = 1; i < tmp2.length; ++i) {
2818                 result[resultOffset + i] = -2 * tmp2[i]; // +/-pi - 2 * atan(y / (r - x))
2819             }
2820 
2821         }
2822 
2823         // fix value to take special cases (+0/+0, +0/-0, -0/+0, -0/-0, +/-infinity) correctly
2824         result[resultOffset] = FastMath.atan2(y[yOffset], x[xOffset]);
2825 
2826     }
2827 
2828     /** Compute two arguments arc tangent of a derivative structure.
2829      * @param y array holding the first operand
2830      * @param yOffset offset of the first operand in its array
2831      * @param x array holding the second operand
2832      * @param xOffset offset of the second operand in its array
2833      * @param result array where result must be stored (for
2834      * two arguments arc tangent the result array <em>cannot</em>
2835      * be the input array)
2836      * @param resultOffset offset of the result in its array
2837      * @param <T> the type of the function parameters and value
2838      */
2839     public <T extends CalculusFieldElement<T>> void atan2(final T[] y, final int yOffset,
2840                                                           final T[] x, final int xOffset,
2841                                                           final T[] result, final int resultOffset) {
2842 
2843         final Field<T> field = y[yOffset].getField();
2844 
2845         // compute r = sqrt(x^2+y^2)
2846         T[] tmp1 = MathArrays.buildArray(field, getSize());
2847         multiply(x, xOffset, x, xOffset, tmp1, 0);      // x^2
2848         T[] tmp2 = MathArrays.buildArray(field, getSize());
2849         multiply(y, yOffset, y, yOffset, tmp2, 0);      // y^2
2850         add(tmp1, 0, tmp2, 0, tmp2, 0);                 // x^2 + y^2
2851         rootN(tmp2, 0, 2, tmp1, 0);                     // r = sqrt(x^2 + y^2)
2852 
2853         if (x[xOffset].getReal() >= 0) {
2854 
2855             // compute atan2(y, x) = 2 atan(y / (r + x))
2856             add(tmp1, 0, x, xOffset, tmp2, 0);          // r + x
2857             divide(y, yOffset, tmp2, 0, tmp1, 0);       // y /(r + x)
2858             atan(tmp1, 0, tmp2, 0);                     // atan(y / (r + x))
2859             for (int i = 0; i < tmp2.length; ++i) {
2860                 result[resultOffset + i] = tmp2[i].add(tmp2[i]); // 2 * atan(y / (r + x))
2861             }
2862 
2863         } else {
2864 
2865             // compute atan2(y, x) = +/- pi - 2 atan(y / (r - x))
2866             subtract(tmp1, 0, x, xOffset, tmp2, 0);     // r - x
2867             divide(y, yOffset, tmp2, 0, tmp1, 0);       // y /(r - x)
2868             atan(tmp1, 0, tmp2, 0);                     // atan(y / (r - x))
2869             result[resultOffset] = tmp2[0].add(tmp2[0]).negate().
2870                                    add((tmp2[0].getReal() <= 0) ? -FastMath.PI : FastMath.PI); // +/-pi - 2 * atan(y / (r - x))
2871             for (int i = 1; i < tmp2.length; ++i) {
2872                 result[resultOffset + i] = tmp2[i].add(tmp2[i]).negate(); // +/-pi - 2 * atan(y / (r - x))
2873             }
2874 
2875         }
2876 
2877         // fix value to take special cases (+0/+0, +0/-0, -0/+0, -0/-0, +/-infinity) correctly
2878         result[resultOffset] = y[yOffset].atan2(x[xOffset]);
2879 
2880     }
2881 
2882     /** Compute hyperbolic cosine of a derivative structure.
2883      * @param operand array holding the operand
2884      * @param operandOffset offset of the operand in its array
2885      * @param result array where result must be stored (for
2886      * hyperbolic cosine the result array <em>cannot</em> be the input
2887      * array)
2888      * @param resultOffset offset of the result in its array
2889      */
2890     public void cosh(final double[] operand, final int operandOffset,
2891                      final double[] result, final int resultOffset) {
2892 
2893         // create the function value and derivatives
2894         double[] function = new double[1 + order];
2895         function[0] = FastMath.cosh(operand[operandOffset]);
2896         if (order > 0) {
2897             function[1] = FastMath.sinh(operand[operandOffset]);
2898             for (int i = 2; i <= order; ++i) {
2899                 function[i] = function[i - 2];
2900             }
2901         }
2902 
2903         // apply function composition
2904         compose(operand, operandOffset, function, result, resultOffset);
2905 
2906     }
2907 
2908     /** Compute hyperbolic cosine of a derivative structure.
2909      * @param operand array holding the operand
2910      * @param operandOffset offset of the operand in its array
2911      * @param result array where result must be stored (for
2912      * hyperbolic cosine the result array <em>cannot</em> be the input
2913      * array)
2914      * @param resultOffset offset of the result in its array
2915      * @param <T> the type of the function parameters and value
2916      */
2917     public <T extends CalculusFieldElement<T>> void cosh(final T[] operand, final int operandOffset,
2918                                                          final T[] result, final int resultOffset) {
2919 
2920         final Field<T> field = operand[operandOffset].getField();
2921 
2922         // create the function value and derivatives
2923         T[] function = MathArrays.buildArray(field, 1 + order);
2924         function[0] = operand[operandOffset].cosh();
2925         if (order > 0) {
2926             function[1] = operand[operandOffset].sinh();
2927             for (int i = 2; i <= order; ++i) {
2928                 function[i] = function[i - 2];
2929             }
2930         }
2931 
2932         // apply function composition
2933         compose(operand, operandOffset, function, result, resultOffset);
2934 
2935     }
2936 
2937     /** Compute hyperbolic sine of a derivative structure.
2938      * @param operand array holding the operand
2939      * @param operandOffset offset of the operand in its array
2940      * @param result array where result must be stored (for
2941      * hyperbolic sine the result array <em>cannot</em> be the input
2942      * array)
2943      * @param resultOffset offset of the result in its array
2944      */
2945     public void sinh(final double[] operand, final int operandOffset,
2946                      final double[] result, final int resultOffset) {
2947 
2948         // create the function value and derivatives
2949         double[] function = new double[1 + order];
2950         function[0] = FastMath.sinh(operand[operandOffset]);
2951         if (order > 0) {
2952             function[1] = FastMath.cosh(operand[operandOffset]);
2953             for (int i = 2; i <= order; ++i) {
2954                 function[i] = function[i - 2];
2955             }
2956         }
2957 
2958         // apply function composition
2959         compose(operand, operandOffset, function, result, resultOffset);
2960 
2961     }
2962 
2963     /** Compute hyperbolic sine of a derivative structure.
2964      * @param operand array holding the operand
2965      * @param operandOffset offset of the operand in its array
2966      * @param result array where result must be stored (for
2967      * hyperbolic sine the result array <em>cannot</em> be the input
2968      * array)
2969      * @param resultOffset offset of the result in its array
2970      * @param <T> the type of the function parameters and value
2971      */
2972     public <T extends CalculusFieldElement<T>> void sinh(final T[] operand, final int operandOffset,
2973                                                          final T[] result, final int resultOffset) {
2974 
2975         final Field<T> field = operand[operandOffset].getField();
2976 
2977         // create the function value and derivatives
2978         T[] function = MathArrays.buildArray(field, 1 + order);
2979         function[0] = operand[operandOffset].sinh();
2980         if (order > 0) {
2981             function[1] = operand[operandOffset].cosh();
2982             for (int i = 2; i <= order; ++i) {
2983                 function[i] = function[i - 2];
2984             }
2985         }
2986 
2987         // apply function composition
2988         compose(operand, operandOffset, function, result, resultOffset);
2989 
2990     }
2991 
2992     /** Compute combined hyperbolic sine and cosine of a derivative structure.
2993      * @param operand array holding the operand
2994      * @param operandOffset offset of the operand in its array
2995      * @param sinh array where hyperbolic sine must be stored (for
2996      * sine the result array <em>cannot</em> be the input
2997      * array)
2998      * @param sinhOffset offset of the result in its array
2999      * @param cosh array where hyperbolic <em>cannot</em> be the input
3000      * array)
3001      * @param coshOffset offset of the result in its array
3002      * @since 2.0
3003      */
3004     public void sinhCosh(final double[] operand, final int operandOffset,
3005                          final double[] sinh, final int sinhOffset,
3006                          final double[] cosh, final int coshOffset) {
3007 
3008         // create the function value and derivatives
3009         double[] functionSinh = new double[1 + order];
3010         double[] functionCosh = new double[1 + order];
3011         final SinhCosh sinhCosh = FastMath.sinhCosh(operand[operandOffset]);
3012         functionSinh[0] = sinhCosh.sinh();
3013         functionCosh[0] = sinhCosh.cosh();
3014         if (order > 0) {
3015             functionSinh[1] = sinhCosh.cosh();
3016             functionCosh[1] = sinhCosh.sinh();
3017             for (int i = 2; i <= order; ++i) {
3018                 functionSinh[i] = functionSinh[i - 2];
3019                 functionCosh[i] = functionCosh[i - 2];
3020             }
3021         }
3022 
3023         // apply function composition
3024         compose(operand, operandOffset, functionSinh, sinh, sinhOffset);
3025         compose(operand, operandOffset, functionCosh, cosh, coshOffset);
3026 
3027     }
3028 
3029     /** Compute combined hyperbolic sine and cosine of a derivative structure.
3030      * @param operand array holding the operand
3031      * @param operandOffset offset of the operand in its array
3032      * @param sinh array where hyperbolic sine must be stored (for
3033      * sine the result array <em>cannot</em> be the input
3034      * array)
3035      * @param sinhOffset offset of the result in its array
3036      * @param cosh array where hyperbolic cosine must be stored (for
3037      * cosine the result array <em>cannot</em> be the input
3038      * array)
3039      * @param coshOffset offset of the result in its array
3040      * @param <T> the type of the function parameters and value
3041      * @since 1.4
3042      */
3043     public <T extends CalculusFieldElement<T>> void sinhCosh(final T[] operand, final int operandOffset,
3044                                                              final T[] sinh, final int sinhOffset,
3045                                                              final T[] cosh, final int coshOffset) {
3046 
3047         final Field<T> field = operand[operandOffset].getField();
3048 
3049         // create the function value and derivatives
3050         T[] functionSinh = MathArrays.buildArray(field, 1 + order);
3051         T[] functionCosh = MathArrays.buildArray(field, 1 + order);
3052         final FieldSinhCosh<T> sinhCosh = FastMath.sinhCosh(operand[operandOffset]);
3053         functionSinh[0] = sinhCosh.sinh();
3054         functionCosh[0] = sinhCosh.cosh();
3055         for (int i = 1; i <= order; ++i) {
3056             functionSinh[i] = functionCosh[i - 1];
3057             functionCosh[i] = functionSinh[i - 1];
3058         }
3059 
3060         // apply function composition
3061         compose(operand, operandOffset, functionSinh, sinh, sinhOffset);
3062         compose(operand, operandOffset, functionCosh, cosh, coshOffset);
3063 
3064     }
3065 
3066     /** Compute hyperbolic tangent of a derivative structure.
3067      * @param operand array holding the operand
3068      * @param operandOffset offset of the operand in its array
3069      * @param result array where result must be stored (for
3070      * hyperbolic tangent the result array <em>cannot</em> be the input
3071      * array)
3072      * @param resultOffset offset of the result in its array
3073      */
3074     public void tanh(final double[] operand, final int operandOffset,
3075                      final double[] result, final int resultOffset) {
3076 
3077         // create the function value and derivatives
3078         final double[] function = new double[1 + order];
3079         final double t = FastMath.tanh(operand[operandOffset]);
3080         function[0] = t;
3081 
3082         if (order > 0) {
3083 
3084             // the nth order derivative of tanh has the form:
3085             // dn(tanh(x)/dxn = P_n(tanh(x))
3086             // where P_n(t) is a degree n+1 polynomial with same parity as n+1
3087             // P_0(t) = t, P_1(t) = 1 - t^2, P_2(t) = -2 t (1 - t^2) ...
3088             // the general recurrence relation for P_n is:
3089             // P_n(x) = (1-t^2) P_(n-1)'(t)
3090             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
3091             final double[] p = new double[order + 2];
3092             p[1] = 1;
3093             final double t2 = t * t;
3094             for (int n = 1; n <= order; ++n) {
3095 
3096                 // update and evaluate polynomial P_n(t)
3097                 double v = 0;
3098                 p[n + 1] = -n * p[n];
3099                 for (int k = n + 1; k >= 0; k -= 2) {
3100                     v = v * t2 + p[k];
3101                     if (k > 2) {
3102                         p[k - 2] = (k - 1) * p[k - 1] - (k - 3) * p[k - 3];
3103                     } else if (k == 2) {
3104                         p[0] = p[1];
3105                     }
3106                 }
3107                 if ((n & 0x1) == 0) {
3108                     v *= t;
3109                 }
3110 
3111                 function[n] = v;
3112 
3113             }
3114         }
3115 
3116         // apply function composition
3117         compose(operand, operandOffset, function, result, resultOffset);
3118 
3119     }
3120 
3121     /** Compute hyperbolic tangent of a derivative structure.
3122      * @param operand array holding the operand
3123      * @param operandOffset offset of the operand in its array
3124      * @param result array where result must be stored (for
3125      * hyperbolic tangent the result array <em>cannot</em> be the input
3126      * array)
3127      * @param resultOffset offset of the result in its array
3128      * @param <T> the type of the function parameters and value
3129      */
3130     public <T extends CalculusFieldElement<T>> void tanh(final T[] operand, final int operandOffset,
3131                                                          final T[] result, final int resultOffset) {
3132 
3133         final Field<T> field = operand[operandOffset].getField();
3134 
3135         // create the function value and derivatives
3136         T[] function = MathArrays.buildArray(field, 1 + order);
3137         final T t = operand[operandOffset].tanh();
3138         function[0] = t;
3139 
3140         if (order > 0) {
3141 
3142             // the nth order derivative of tanh has the form:
3143             // dn(tanh(x)/dxn = P_n(tanh(x))
3144             // where P_n(t) is a degree n+1 polynomial with same parity as n+1
3145             // P_0(t) = t, P_1(t) = 1 - t^2, P_2(t) = -2 t (1 - t^2) ...
3146             // the general recurrence relation for P_n is:
3147             // P_n(x) = (1-t^2) P_(n-1)'(t)
3148             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
3149             final T[] p = MathArrays.buildArray(field, order + 2);
3150             p[1] = field.getOne();
3151             final T t2 = t.multiply(t);
3152             for (int n = 1; n <= order; ++n) {
3153 
3154                 // update and evaluate polynomial P_n(t)
3155                 T v = field.getZero();
3156                 p[n + 1] = p[n].multiply(-n);
3157                 for (int k = n + 1; k >= 0; k -= 2) {
3158                     v = v.multiply(t2).add(p[k]);
3159                     if (k > 2) {
3160                         p[k - 2] = p[k - 1].multiply(k - 1).subtract(p[k - 3].multiply(k - 3));
3161                     } else if (k == 2) {
3162                         p[0] = p[1];
3163                     }
3164                 }
3165                 if ((n & 0x1) == 0) {
3166                     v = v.multiply(t);
3167                 }
3168 
3169                 function[n] = v;
3170 
3171             }
3172         }
3173 
3174         // apply function composition
3175         compose(operand, operandOffset, function, result, resultOffset);
3176 
3177     }
3178 
3179     /** Compute inverse hyperbolic cosine of a derivative structure.
3180      * @param operand array holding the operand
3181      * @param operandOffset offset of the operand in its array
3182      * @param result array where result must be stored (for
3183      * inverse hyperbolic cosine the result array <em>cannot</em> be the input
3184      * array)
3185      * @param resultOffset offset of the result in its array
3186      */
3187     public void acosh(final double[] operand, final int operandOffset,
3188                       final double[] result, final int resultOffset) {
3189 
3190         // create the function value and derivatives
3191         double[] function = new double[1 + order];
3192         final double x = operand[operandOffset];
3193         function[0] = FastMath.acosh(x);
3194         if (order > 0) {
3195             // the nth order derivative of acosh has the form:
3196             // dn(acosh(x)/dxn = P_n(x) / [x^2 - 1]^((2n-1)/2)
3197             // where P_n(x) is a degree n-1 polynomial with same parity as n-1
3198             // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 + 1 ...
3199             // the general recurrence relation for P_n is:
3200             // P_n(x) = (x^2-1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x)
3201             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
3202             final double[] p = new double[order];
3203             p[0] = 1;
3204             final double x2  = x * x;
3205             final double f   = 1.0 / (x2 - 1);
3206             double coeff = FastMath.sqrt(f);
3207             function[1] = coeff * p[0];
3208             for (int n = 2; n <= order; ++n) {
3209 
3210                 // update and evaluate polynomial P_n(x)
3211                 double v = 0;
3212                 p[n - 1] = (1 - n) * p[n - 2];
3213                 for (int k = n - 1; k >= 0; k -= 2) {
3214                     v = v * x2 + p[k];
3215                     if (k > 2) {
3216                         p[k - 2] = (1 - k) * p[k - 1] + (k - 2 * n) * p[k - 3];
3217                     } else if (k == 2) {
3218                         p[0] = -p[1];
3219                     }
3220                 }
3221                 if ((n & 0x1) == 0) {
3222                     v *= x;
3223                 }
3224 
3225                 coeff *= f;
3226                 function[n] = coeff * v;
3227 
3228             }
3229         }
3230 
3231         // apply function composition
3232         compose(operand, operandOffset, function, result, resultOffset);
3233 
3234     }
3235 
3236     /** Compute inverse hyperbolic cosine of a derivative structure.
3237      * @param operand array holding the operand
3238      * @param operandOffset offset of the operand in its array
3239      * @param result array where result must be stored (for
3240      * inverse hyperbolic cosine the result array <em>cannot</em> be the input
3241      * array)
3242      * @param resultOffset offset of the result in its array
3243      * @param <T> the type of the function parameters and value
3244      */
3245     public <T extends CalculusFieldElement<T>> void acosh(final T[] operand, final int operandOffset,
3246                                                           final T[] result, final int resultOffset) {
3247 
3248         final Field<T> field = operand[operandOffset].getField();
3249 
3250         // create the function value and derivatives
3251         T[] function = MathArrays.buildArray(field, 1 + order);
3252         final T x = operand[operandOffset];
3253         function[0] = x.acosh();
3254         if (order > 0) {
3255             // the nth order derivative of acosh has the form:
3256             // dn(acosh(x)/dxn = P_n(x) / [x^2 - 1]^((2n-1)/2)
3257             // where P_n(x) is a degree n-1 polynomial with same parity as n-1
3258             // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 + 1 ...
3259             // the general recurrence relation for P_n is:
3260             // P_n(x) = (x^2-1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x)
3261             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
3262             final T[] p = MathArrays.buildArray(field, order);
3263             p[0] = field.getOne();
3264             final T x2  = x.square();
3265             final T f   = x2.subtract(1).reciprocal();
3266             T coeff = f.sqrt();
3267             function[1] = coeff.multiply(p[0]);
3268             for (int n = 2; n <= order; ++n) {
3269 
3270                 // update and evaluate polynomial P_n(x)
3271                 T v = field.getZero();
3272                 p[n - 1] = p[n - 2].multiply(1 - n);
3273                 for (int k = n - 1; k >= 0; k -= 2) {
3274                     v = v.multiply(x2).add(p[k]);
3275                     if (k > 2) {
3276                         p[k - 2] = p[k - 1].multiply(1 - k).add(p[k - 3].multiply(k - 2 * n));
3277                     } else if (k == 2) {
3278                         p[0] = p[1].negate();
3279                     }
3280                 }
3281                 if ((n & 0x1) == 0) {
3282                     v = v.multiply(x);
3283                 }
3284 
3285                 coeff = coeff.multiply(f);
3286                 function[n] = coeff.multiply(v);
3287 
3288             }
3289         }
3290 
3291         // apply function composition
3292         compose(operand, operandOffset, function, result, resultOffset);
3293 
3294     }
3295 
3296     /** Compute inverse hyperbolic sine of a derivative structure.
3297      * @param operand array holding the operand
3298      * @param operandOffset offset of the operand in its array
3299      * @param result array where result must be stored (for
3300      * inverse hyperbolic sine the result array <em>cannot</em> be the input
3301      * array)
3302      * @param resultOffset offset of the result in its array
3303      */
3304     public void asinh(final double[] operand, final int operandOffset,
3305                       final double[] result, final int resultOffset) {
3306 
3307         // create the function value and derivatives
3308         double[] function = new double[1 + order];
3309         final double x = operand[operandOffset];
3310         function[0] = FastMath.asinh(x);
3311         if (order > 0) {
3312             // the nth order derivative of asinh has the form:
3313             // dn(asinh(x)/dxn = P_n(x) / [x^2 + 1]^((2n-1)/2)
3314             // where P_n(x) is a degree n-1 polynomial with same parity as n-1
3315             // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 - 1 ...
3316             // the general recurrence relation for P_n is:
3317             // P_n(x) = (x^2+1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x)
3318             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
3319             final double[] p = new double[order];
3320             p[0] = 1;
3321             final double x2    = x * x;
3322             final double f     = 1.0 / (1 + x2);
3323             double coeff = FastMath.sqrt(f);
3324             function[1] = coeff * p[0];
3325             for (int n = 2; n <= order; ++n) {
3326 
3327                 // update and evaluate polynomial P_n(x)
3328                 double v = 0;
3329                 p[n - 1] = (1 - n) * p[n - 2];
3330                 for (int k = n - 1; k >= 0; k -= 2) {
3331                     v = v * x2 + p[k];
3332                     if (k > 2) {
3333                         p[k - 2] = (k - 1) * p[k - 1] + (k - 2 * n) * p[k - 3];
3334                     } else if (k == 2) {
3335                         p[0] = p[1];
3336                     }
3337                 }
3338                 if ((n & 0x1) == 0) {
3339                     v *= x;
3340                 }
3341 
3342                 coeff *= f;
3343                 function[n] = coeff * v;
3344 
3345             }
3346         }
3347 
3348         // apply function composition
3349         compose(operand, operandOffset, function, result, resultOffset);
3350 
3351     }
3352 
3353     /** Compute inverse hyperbolic sine of a derivative structure.
3354      * @param operand array holding the operand
3355      * @param operandOffset offset of the operand in its array
3356      * @param result array where result must be stored (for
3357      * inverse hyperbolic sine the result array <em>cannot</em> be the input
3358      * array)
3359      * @param resultOffset offset of the result in its array
3360      * @param <T> the type of the function parameters and value
3361      */
3362     public <T extends CalculusFieldElement<T>> void asinh(final T[] operand, final int operandOffset,
3363                                                           final T[] result, final int resultOffset) {
3364 
3365         final Field<T> field = operand[operandOffset].getField();
3366 
3367         // create the function value and derivatives
3368         T[] function = MathArrays.buildArray(field, 1 + order);
3369         final T x = operand[operandOffset];
3370         function[0] = x.asinh();
3371         if (order > 0) {
3372             // the nth order derivative of asinh has the form:
3373             // dn(asinh(x)/dxn = P_n(x) / [x^2 + 1]^((2n-1)/2)
3374             // where P_n(x) is a degree n-1 polynomial with same parity as n-1
3375             // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 - 1 ...
3376             // the general recurrence relation for P_n is:
3377             // P_n(x) = (x^2+1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x)
3378             // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
3379             final T[] p = MathArrays.buildArray(field, order);
3380             p[0] = field.getOne();
3381             final T x2    = x.square();
3382             final T f     = x2.add(1).reciprocal();
3383             T coeff = f.sqrt();
3384             function[1] = coeff.multiply(p[0]);
3385             for (int n = 2; n <= order; ++n) {
3386 
3387                 // update and evaluate polynomial P_n(x)
3388                 T v = field.getZero();
3389                 p[n - 1] = p[n - 2].multiply(1 - n);
3390                 for (int k = n - 1; k >= 0; k -= 2) {
3391                     v = v.multiply(x2).add(p[k]);
3392                     if (k > 2) {
3393                         p[k - 2] = p[k - 1].multiply(k - 1).add(p[k - 3].multiply(k - 2 * n));
3394                     } else if (k == 2) {
3395                         p[0] = p[1];
3396                     }
3397                 }
3398                 if ((n & 0x1) == 0) {
3399                     v = v.multiply(x);
3400                 }
3401 
3402                 coeff = coeff.multiply(f);
3403                 function[n] = coeff.multiply(v);
3404 
3405             }
3406         }
3407 
3408         // apply function composition
3409         compose(operand, operandOffset, function, result, resultOffset);
3410 
3411     }
3412 
3413     /** Compute inverse hyperbolic tangent of a derivative structure.
3414      * @param operand array holding the operand
3415      * @param operandOffset offset of the operand in its array
3416      * @param result array where result must be stored (for
3417      * inverse hyperbolic tangent the result array <em>cannot</em> be the input
3418      * array)
3419      * @param resultOffset offset of the result in its array
3420      */
3421     public void atanh(final double[] operand, final int operandOffset,
3422                       final double[] result, final int resultOffset) {
3423 
3424         // create the function value and derivatives
3425         double[] function = new double[1 + order];
3426         final double x = operand[operandOffset];
3427         function[0] = FastMath.atanh(x);
3428         if (order > 0) {
3429             // the nth order derivative of atanh has the form:
3430             // dn(atanh(x)/dxn = Q_n(x) / (1 - x^2)^n
3431             // where Q_n(x) is a degree n-1 polynomial with same parity as n-1
3432             // Q_1(x) = 1, Q_2(x) = 2x, Q_3(x) = 6x^2 + 2 ...
3433             // the general recurrence relation for Q_n is:
3434             // Q_n(x) = (1-x^2) Q_(n-1)'(x) + 2(n-1) x Q_(n-1)(x)
3435             // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array
3436             final double[] q = new double[order];
3437             q[0] = 1;
3438             final double x2 = x * x;
3439             final double f  = 1.0 / (1 - x2);
3440             double coeff = f;
3441             function[1] = coeff * q[0];
3442             for (int n = 2; n <= order; ++n) {
3443 
3444                 // update and evaluate polynomial Q_n(x)
3445                 double v = 0;
3446                 q[n - 1] = n * q[n - 2];
3447                 for (int k = n - 1; k >= 0; k -= 2) {
3448                     v = v * x2 + q[k];
3449                     if (k > 2) {
3450                         q[k - 2] = (k - 1) * q[k - 1] + (2 * n - k + 1) * q[k - 3];
3451                     } else if (k == 2) {
3452                         q[0] = q[1];
3453                     }
3454                 }
3455                 if ((n & 0x1) == 0) {
3456                     v *= x;
3457                 }
3458 
3459                 coeff *= f;
3460                 function[n] = coeff * v;
3461 
3462             }
3463         }
3464 
3465         // apply function composition
3466         compose(operand, operandOffset, function, result, resultOffset);
3467 
3468     }
3469 
3470     /** Compute inverse hyperbolic tangent of a derivative structure.
3471      * @param operand array holding the operand
3472      * @param operandOffset offset of the operand in its array
3473      * @param result array where result must be stored (for
3474      * inverse hyperbolic tangent the result array <em>cannot</em> be the input
3475      * array)
3476      * @param resultOffset offset of the result in its array
3477      * @param <T> the type of the function parameters and value
3478      */
3479     public <T extends CalculusFieldElement<T>> void atanh(final T[] operand, final int operandOffset,
3480                                                           final T[] result, final int resultOffset) {
3481 
3482         final Field<T> field = operand[operandOffset].getField();
3483 
3484         // create the function value and derivatives
3485         T[] function = MathArrays.buildArray(field, 1 + order);
3486         final T x = operand[operandOffset];
3487         function[0] = x.atanh();
3488         if (order > 0) {
3489             // the nth order derivative of atanh has the form:
3490             // dn(atanh(x)/dxn = Q_n(x) / (1 - x^2)^n
3491             // where Q_n(x) is a degree n-1 polynomial with same parity as n-1
3492             // Q_1(x) = 1, Q_2(x) = 2x, Q_3(x) = 6x^2 + 2 ...
3493             // the general recurrence relation for Q_n is:
3494             // Q_n(x) = (1-x^2) Q_(n-1)'(x) + 2(n-1) x Q_(n-1)(x)
3495             // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array
3496             final T[] q = MathArrays.buildArray(field, order);
3497             q[0] = field.getOne();
3498             final T x2 = x.square();
3499             final T f  =x2.subtract(1).negate().reciprocal();
3500             T coeff = f;
3501             function[1] = coeff.multiply(q[0]);
3502             for (int n = 2; n <= order; ++n) {
3503 
3504                 // update and evaluate polynomial Q_n(x)
3505                 T v = field.getZero();
3506                 q[n - 1] = q[n - 2].multiply(n);
3507                 for (int k = n - 1; k >= 0; k -= 2) {
3508                     v = v.multiply(x2).add(q[k]);
3509                     if (k > 2) {
3510                         q[k - 2] = q[k - 1].multiply(k - 1).add(q[k - 3].multiply(2 * n - k + 1));
3511                     } else if (k == 2) {
3512                         q[0] = q[1];
3513                     }
3514                 }
3515                 if ((n & 0x1) == 0) {
3516                     v = v.multiply(x);
3517                 }
3518 
3519                 coeff = coeff.multiply(f);
3520                 function[n] = coeff.multiply(v);
3521 
3522             }
3523         }
3524 
3525         // apply function composition
3526         compose(operand, operandOffset, function, result, resultOffset);
3527 
3528     }
3529 
3530     /** Compute composition of a derivative structure by a function.
3531      * @param operand array holding the operand
3532      * @param operandOffset offset of the operand in its array
3533      * @param f array of value and derivatives of the function at
3534      * the current point (i.e. at {@code operand[operandOffset]}).
3535      * @param result array where result must be stored (for
3536      * composition the result array <em>cannot</em> be the input
3537      * array)
3538      * @param resultOffset offset of the result in its array
3539      */
3540     public void compose(final double[] operand, final int operandOffset, final double[] f,
3541                         final double[] result, final int resultOffset) {
3542         for (int i = 0; i < compIndirection.length; ++i) {
3543             final UnivariateCompositionMapper[] mappingI = compIndirection[i];
3544             double r = 0;
3545             for (UnivariateCompositionMapper mapping : mappingI) {
3546                 double product = mapping.getCoeff() * f[mapping.fIndex];
3547                 for (int k = 0; k < mapping.dsIndices.length; ++k) {
3548                     product *= operand[operandOffset + mapping.dsIndices[k]];
3549                 }
3550                 r += product;
3551             }
3552             result[resultOffset + i] = r;
3553         }
3554     }
3555 
3556     /** Compute composition of a derivative structure by a function.
3557      * @param operand array holding the operand
3558      * @param operandOffset offset of the operand in its array
3559      * @param f array of value and derivatives of the function at
3560      * the current point (i.e. at {@code operand[operandOffset]}).
3561      * @param result array where result must be stored (for
3562      * composition the result array <em>cannot</em> be the input
3563      * array)
3564      * @param resultOffset offset of the result in its array
3565      * @param <T> the type of the function parameters and value
3566      */
3567     public <T extends CalculusFieldElement<T>> void compose(final T[] operand, final int operandOffset, final T[] f,
3568                                                             final T[] result, final int resultOffset) {
3569         final T zero = f[0].getField().getZero();
3570         for (int i = 0; i < compIndirection.length; ++i) {
3571             final UnivariateCompositionMapper[] mappingI = compIndirection[i];
3572             T r = zero;
3573             for (UnivariateCompositionMapper mapping : mappingI) {
3574                 T product = f[mapping.fIndex].multiply(mapping.getCoeff());
3575                 for (int k = 0; k < mapping.dsIndices.length; ++k) {
3576                     product = product.multiply(operand[operandOffset + mapping.dsIndices[k]]);
3577                 }
3578                 r = r.add(product);
3579             }
3580             result[resultOffset + i] = r;
3581         }
3582     }
3583 
3584     /** Compute composition of a derivative structure by a function.
3585      * @param operand array holding the operand
3586      * @param operandOffset offset of the operand in its array
3587      * @param f array of value and derivatives of the function at
3588      * the current point (i.e. at {@code operand[operandOffset]}).
3589      * @param result array where result must be stored (for
3590      * composition the result array <em>cannot</em> be the input
3591      * array)
3592      * @param resultOffset offset of the result in its array
3593      * @param <T> the type of the function parameters and value
3594      */
3595     public <T extends CalculusFieldElement<T>> void compose(final T[] operand, final int operandOffset, final double[] f,
3596                                                             final T[] result, final int resultOffset) {
3597         final T zero = operand[operandOffset].getField().getZero();
3598         for (int i = 0; i < compIndirection.length; ++i) {
3599             final UnivariateCompositionMapper[] mappingI = compIndirection[i];
3600             T r = zero;
3601             for (UnivariateCompositionMapper mapping : mappingI) {
3602                 T product = zero.add(f[mapping.fIndex] * mapping.getCoeff());
3603                 for (int k = 0; k < mapping.dsIndices.length; ++k) {
3604                     product = product.multiply(operand[operandOffset + mapping.dsIndices[k]]);
3605                 }
3606                 r = r.add(product);
3607             }
3608             result[resultOffset + i] = r;
3609         }
3610     }
3611 
3612     /** Evaluate Taylor expansion of a derivative structure.
3613      * @param ds array holding the derivative structure
3614      * @param dsOffset offset of the derivative structure in its array
3615      * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
3616      * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
3617      * @throws MathRuntimeException if factorials becomes too large
3618      */
3619     public double taylor(final double[] ds, final int dsOffset, final double ... delta)
3620        throws MathRuntimeException {
3621         double value = 0;
3622         for (int i = getSize() - 1; i >= 0; --i) {
3623             final int[] orders = derivativesOrders[i];
3624             double term = ds[dsOffset + i];
3625             for (int k = 0; k < orders.length; ++k) {
3626                 if (orders[k] > 0) {
3627                     term *= FastMath.pow(delta[k], orders[k]) /
3628                             CombinatoricsUtils.factorial(orders[k]);
3629                 }
3630             }
3631             value += term;
3632         }
3633         return value;
3634     }
3635 
3636     /** Evaluate Taylor expansion of a derivative structure.
3637      * @param ds array holding the derivative structure
3638      * @param dsOffset offset of the derivative structure in its array
3639      * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
3640      * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
3641      * @throws MathRuntimeException if factorials becomes too large
3642      * @param <T> the type of the function parameters and value
3643      */
3644     @SafeVarargs
3645     public final <T extends CalculusFieldElement<T>> T taylor(final T[] ds, final int dsOffset,
3646                                                               final T ... delta)
3647        throws MathRuntimeException {
3648         final Field<T> field = ds[dsOffset].getField();
3649         T value = field.getZero();
3650         for (int i = getSize() - 1; i >= 0; --i) {
3651             final int[] orders = derivativesOrders[i];
3652             T term = ds[dsOffset + i];
3653             for (int k = 0; k < orders.length; ++k) {
3654                 if (orders[k] > 0) {
3655                     term = term.multiply(delta[k].pow(orders[k]).
3656                                          divide(CombinatoricsUtils.factorial(orders[k])));
3657                 }
3658             }
3659             value = value.add(term);
3660         }
3661         return value;
3662     }
3663 
3664     /** Evaluate Taylor expansion of a derivative structure.
3665      * @param ds array holding the derivative structure
3666      * @param dsOffset offset of the derivative structure in its array
3667      * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
3668      * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
3669      * @throws MathRuntimeException if factorials becomes too large
3670      * @param <T> the type of the function parameters and value
3671      */
3672     public <T extends CalculusFieldElement<T>> T taylor(final T[] ds, final int dsOffset,
3673                                                         final double ... delta)
3674        throws MathRuntimeException {
3675         final Field<T> field = ds[dsOffset].getField();
3676         T value = field.getZero();
3677         for (int i = getSize() - 1; i >= 0; --i) {
3678             final int[] orders = derivativesOrders[i];
3679             T term = ds[dsOffset + i];
3680             for (int k = 0; k < orders.length; ++k) {
3681                 if (orders[k] > 0) {
3682                     term = term.multiply(field.getZero().newInstance(delta[k]).pow(orders[k]).
3683                                          divide(CombinatoricsUtils.factorial(orders[k])));
3684                 }
3685             }
3686             value = value.add(term);
3687         }
3688         return value;
3689     }
3690 
3691     /** Rebase derivative structure with respect to low level parameter functions.
3692      * @param ds array holding the derivative structure
3693      * @param dsOffset offset of the derivative structure in its array
3694      * @param baseCompiler compiler associated with the low level parameter functions
3695      * @param p array holding the low level parameter functions (one flat array)
3696      * @param result array where result must be stored (for
3697      * composition the result array <em>cannot</em> be the input
3698      * @param resultOffset offset of the result in its array
3699      * @since 2.2
3700      */
3701     public void rebase(final double[] ds, final int dsOffset,
3702                        final DSCompiler baseCompiler, double[] p,
3703                        final double[] result, final int resultOffset) {
3704         final MultivariateCompositionMapper[][] rebaser = getRebaser(baseCompiler);
3705         for (int i = 0; i < rebaser.length; ++i) {
3706             final MultivariateCompositionMapper[] mappingI = rebaser[i];
3707             double r = 0;
3708             for (MultivariateCompositionMapper mapping : mappingI) {
3709                 double product = mapping.getCoeff() * ds[dsOffset + mapping.dsIndex];
3710                 for (int k = 0; k < mapping.productIndices.length; ++k) {
3711                     product *= p[mapping.productIndices[k]];
3712                 }
3713                 r += product;
3714             }
3715             result[resultOffset + i] = r;
3716         }
3717     }
3718 
3719     /** Rebase derivative structure with respect to low level parameter functions.
3720      * @param <T> type of the field elements
3721      * @param ds array holding the derivative structure
3722      * @param dsOffset offset of the derivative structure in its array
3723      * @param baseCompiler compiler associated with the low level parameter functions
3724      * @param p array holding the low level parameter functions (one flat array)
3725      * @param result array where result must be stored (for
3726      * composition the result array <em>cannot</em> be the input
3727      * @param resultOffset offset of the result in its array
3728      * @since 2.2
3729      */
3730     public <T extends CalculusFieldElement<T>> void rebase(final T[] ds, final int dsOffset,
3731                                                            final DSCompiler baseCompiler, T[] p,
3732                                                            final T[] result, final int resultOffset) {
3733         final MultivariateCompositionMapper[][] rebaser = getRebaser(baseCompiler);
3734         for (int i = 0; i < rebaser.length; ++i) {
3735             final MultivariateCompositionMapper[] mappingI = rebaser[i];
3736             T r = ds[0].getField().getZero();
3737             for (MultivariateCompositionMapper mapping : mappingI) {
3738                 T product =  ds[dsOffset + mapping.dsIndex].multiply(mapping.getCoeff());
3739                 for (int k = 0; k < mapping.productIndices.length; ++k) {
3740                     product = product.multiply(p[mapping.productIndices[k]]);
3741                 }
3742                 r = r.add(product);
3743             }
3744             result[resultOffset + i] = r;
3745         }
3746     }
3747 
3748     /** Check rules set compatibility.
3749      * @param compiler other compiler to check against instance
3750      * @exception MathIllegalArgumentException if number of free parameters or orders are inconsistent
3751      */
3752     public void checkCompatibility(final DSCompiler compiler)
3753         throws MathIllegalArgumentException {
3754         MathUtils.checkDimension(parameters, compiler.parameters);
3755         MathUtils.checkDimension(order, compiler.order);
3756     }
3757 
3758     /** Combine terms with similar derivation orders.
3759      * @param <T> type of the field elements
3760      * @param terms list of terms
3761      * @return combined array
3762      */
3763     @SuppressWarnings("unchecked")
3764     private static <T extends AbstractMapper<T>> T[] combineSimilarTerms(final List<T> terms) {
3765 
3766         final List<T> combined = new ArrayList<>(terms.size());
3767 
3768         for (int j = 0; j < terms.size(); ++j) {
3769             final T termJ = terms.get(j);
3770             if (termJ.getCoeff() > 0) {
3771                 for (int k = j + 1; k < terms.size(); ++k) {
3772                     final T termK = terms.get(k);
3773                     if (termJ.isSimilar(termK)) {
3774                         // combine terms
3775                         termJ.setCoeff(termJ.getCoeff() + termK.getCoeff());
3776                         // make sure we will skip other term later on in the outer loop
3777                         termK.setCoeff(0);
3778                     }
3779                 }
3780                 combined.add(termJ);
3781             }
3782         }
3783 
3784         return combined.toArray((T[]) Array.newInstance(terms.get(0).getClass(), combined.size()));
3785 
3786     }
3787 
3788     /** Base mapper.
3789      * @param <T> type of the field elements
3790      * @since 2.2
3791      */
3792     private abstract static class AbstractMapper<T extends AbstractMapper<T>> {
3793 
3794         /** Multiplication coefficient. */
3795         private int coeff;
3796 
3797         /** Simple constructor.
3798          * @param coeff multiplication coefficient
3799          */
3800         AbstractMapper(final int coeff) {
3801             this.coeff    = coeff;
3802         }
3803 
3804         /** Set the multiplication coefficient.
3805          * @param coeff new coefficient
3806          */
3807         public void setCoeff(final int coeff) {
3808             this.coeff = coeff;
3809         }
3810 
3811         /** Get the multiplication coefficient.
3812          * @return multiplication coefficient
3813          */
3814         public int getCoeff() {
3815             return coeff;
3816         }
3817 
3818         /** Check if another instance if correspond to term with similar derivation orders.
3819          * @param other other instance to check
3820          * @return true if instances are similar
3821          */
3822         protected abstract boolean isSimilar(T other);
3823 
3824     }
3825 
3826     /** Multiplication mapper.
3827      * @since 2.2
3828      */
3829     private static class MultiplicationMapper extends AbstractMapper<MultiplicationMapper> {
3830 
3831         /** Left hand side index. */
3832         private final int lhsIndex;
3833 
3834         /** Right hand side index. */
3835         private final int rhsIndex;
3836 
3837         /** Simple constructor.
3838          * @param coeff multiplication coefficient
3839          * @param lhsIndex left hand side index
3840          * @param rhsIndex right hand side index
3841          */
3842         MultiplicationMapper(final int coeff, final int lhsIndex, final int rhsIndex) {
3843             super(coeff);
3844             this.lhsIndex = lhsIndex;
3845             this.rhsIndex = rhsIndex;
3846         }
3847 
3848         /** {@inheritDoc} */
3849         @Override
3850         public boolean isSimilar(final MultiplicationMapper other) {
3851             return lhsIndex == other.lhsIndex && rhsIndex == other.rhsIndex;
3852         }
3853 
3854     }
3855 
3856     /** Univariate composition mapper.
3857      * @since 2.2
3858      */
3859     private static class UnivariateCompositionMapper extends AbstractMapper<UnivariateCompositionMapper> {
3860 
3861         /** Univariate derivative index. */
3862         private final int fIndex;
3863 
3864         /** Derivative structure indices. */
3865         private final int[] dsIndices;
3866 
3867         /** Simple constructor.
3868          * @param coeff multiplication coefficient
3869          * @param fIndex univariate derivative index
3870          * @param dsIndices derivative structure indices
3871          */
3872         UnivariateCompositionMapper(final int coeff, final int fIndex, final int[] dsIndices) {
3873             super(coeff);
3874             this.fIndex    = fIndex;
3875             this.dsIndices = dsIndices.clone();
3876         }
3877 
3878         /** Sort the derivatives structures indices.
3879          */
3880         public void sort() {
3881             Arrays.sort(dsIndices);
3882         }
3883 
3884         /** {@inheritDoc} */
3885         @Override
3886         public boolean isSimilar(final UnivariateCompositionMapper other) {
3887 
3888             if (fIndex == other.fIndex && dsIndices.length == other.dsIndices.length) {
3889 
3890                 for (int j = 0; j < dsIndices.length; ++j) {
3891                     if (dsIndices[j] != other.dsIndices[j]) {
3892                         return false;
3893                     }
3894                 }
3895 
3896                 return true;
3897 
3898             }
3899 
3900             return false;
3901 
3902         }
3903 
3904     }
3905 
3906     /** Multivariate composition mapper.
3907      * @since 2.2
3908      */
3909     private static class MultivariateCompositionMapper extends AbstractMapper<MultivariateCompositionMapper> {
3910 
3911         /** Multivariate derivative index. */
3912         private final int dsIndex;
3913 
3914         /** Indices of the intermediate variables derivatives products. */
3915         private final int[] productIndices;
3916 
3917         /** Simple constructor.
3918          * @param coeff multiplication coefficient
3919          * @param dsIndex multivariate derivative index of ∂ₘf/∂pᵢ⋯∂pⱼ
3920          * @param productIndices indices of intermediate partial derivatives ∂pᵢ/∂qₘ⋯∂qₙ
3921          */
3922         MultivariateCompositionMapper(final int coeff, final int dsIndex, final int[] productIndices) {
3923             super(coeff);
3924             this.dsIndex        = dsIndex;
3925             this.productIndices = productIndices.clone();
3926         }
3927 
3928         /** Sort the indices of the intermediate variables derivatives products.
3929          */
3930         public void sort() {
3931             Arrays.sort(productIndices);
3932         }
3933 
3934         /** {@inheritDoc} */
3935         @Override
3936         public boolean isSimilar(final MultivariateCompositionMapper other) {
3937 
3938             if (dsIndex == other.dsIndex && productIndices.length == other.productIndices.length) {
3939 
3940                 for (int j = 0; j < productIndices.length; ++j) {
3941                     if (productIndices[j] != other.productIndices[j]) {
3942                         return false;
3943                     }
3944                 }
3945 
3946                 return true;
3947 
3948             }
3949 
3950             return false;
3951 
3952         }
3953 
3954     }
3955 
3956 }