DSCompiler.java

  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.  */

  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;

  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;

  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;

  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 {

  138.     /** Array of all compilers created so far. */
  139.     private static AtomicReference<DSCompiler[][]> compilers = new AtomicReference<>(null);

  141.     /** Number of free parameters. */
  142.     private final int parameters;

  144.     /** Derivation order. */
  145.     private final int order;

  147.     /** Number of partial derivatives (including the single 0 order derivative element). */
  148.     private final int[][] sizes;

  150.     /** Orders array for partial derivatives. */
  151.     private final int[][] derivativesOrders;

  153.     /** Sum of orders array for partial derivatives. */
  154.     private final int[] derivativesOrdersSum;

  156.     /** Indirection array of the lower derivative elements. */
  157.     private final int[] lowerIndirection;

  159.     /** Indirection arrays for multiplication. */
  160.     private final MultiplicationMapper[][] multIndirection;

  162.     /** Indirection arrays for univariate function composition. */
  163.     private final UnivariateCompositionMapper[][] compIndirection;

  165.     /** Indirection arrays for multivariate function rebasing. */
  166.     private final List<MultivariateCompositionMapper[][]> rebaseIndirection;

  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 {

  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);

  194.         this.rebaseIndirection = new ArrayList<>();
  195.     }

  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 {

  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.         }

  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];

  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.         }

  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.         }

  238.         // atomically reset the cached compilers array
  239.         compilers.compareAndSet(cache, newCache);

  241.         return newCache[parameters][order];

  243.     }

  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) {

  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.         }

  265.         return sizes;

  267.     }

  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) {

  280.         if (parameters == 0 || order == 0) {
  281.             return new int[1][parameters];
  282.         }

  284.         final int vSize = valueCompiler.derivativesOrders.length;
  285.         final int dSize = derivativeCompiler.derivativesOrders.length;
  286.         final int[][] derivativesOrders = new int[vSize + dSize][parameters];

  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.         }

  296.         // set up the indices for the derivative part
  297.         for (int i = 0; i < dSize; ++i) {

  299.             // copy the indices
  300.             System.arraycopy(derivativeCompiler.derivativesOrders[i], 0,
  301.                              derivativesOrders[vSize + i], 0,
  302.                              parameters);

  304.             // increment the derivation order for the last parameter
  305.             derivativesOrders[vSize + i][parameters - 1]++;

  307.         }

  309.         return derivativesOrders;

  311.     }

  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) {

  319.         final int[] derivativesOrdersSum = new int[derivativesOrders.length];

  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.         }

  328.         return derivativesOrdersSum;

  330.     }

  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) {

  347.         if (parameters == 0 || order <= 1) {
  348.             return new int[] { 0 };
  349.         }

  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.         }

  360.         return lowerIndirection;

  362.     }

  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) {

  382.         if (parameters == 0 || order == 0) {
  383.             return new MultiplicationMapper[][] { { new MultiplicationMapper(1, 0, 0) } };
  384.         }

  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][];

  391.         System.arraycopy(valueCompiler.multIndirection, 0, multIndirection, 0, vSize);

  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.         }

  403.         return multIndirection;

  405.     }

  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 {

  429.         if (parameters == 0 || order == 0) {
  430.             return new UnivariateCompositionMapper[][] { { new UnivariateCompositionMapper(1, 0, new int[0]) } };
  431.         }

  433.         final int vSize = valueCompiler.compIndirection.length;
  434.         final int dSize = derivativeCompiler.compIndirection.length;
  435.         final UnivariateCompositionMapper[][] compIndirection = new UnivariateCompositionMapper[vSize + dSize][];

  437.         // the composition rules from the value part can be reused as is
  438.         System.arraycopy(valueCompiler.compIndirection, 0, compIndirection, 0, vSize);

  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]) {

  448.                 // handle term p * f_k(g(x)) * g_l1(x) * g_l2(x) * ... * g_lp(x)

  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);

  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.                 }

  489.             }

  491.             // combine terms with similar derivation orders
  492.             compIndirection[vSize + i] = combineSimilarTerms(row);

  494.         }

  496.         return compIndirection;

  498.     }

  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) {

  508.             final int m = baseCompiler.getFreeParameters();
  509.             while (rebaseIndirection.size() <= m) {
  510.                 rebaseIndirection.add(null);
  511.             }

  513.             if (rebaseIndirection.get(m) == null) {
  514.                 // we need to create the rebaser from instance to m base variables

  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.                 }

  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);

  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<>();

  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.                         }

  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)];

  568.                         // apply recursion formula
  569.                         for (final MultivariateCompositionMapper lowerTerm : lowerRow) {

  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.                             }

  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.                             }

  581.                         }

  583.                         // simplifies and store the completed entry
  584.                         rebaser[k] = combineSimilarTerms(row);

  586.                     }

  588.                 }

  590.                 rebaseIndirection.set(m, rebaser);

  592.             }

  594.             return rebaseIndirection.get(m);

  596.         }
  597.     }

  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) {

  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);

  612.         // allocate array for rebaser at current order
  613.         final int baseSize = baseCompiler.getSize();
  614.         final MultivariateCompositionMapper[][] rebaser = new MultivariateCompositionMapper[baseSize][];

  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.         }

  639.         return rebaser;

  641.     }

  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) {

  654.         // differentiate the term ∂fⁿ⁻¹/∂pᵤ⋯∂pᵥ with respect to pi
  655.         final int[] termOrders = derivativesOrders[lowerTerm.dsIndex].clone();
  656.         termOrders[i]++;

  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);

  667.         // generate the differentiated term
  668.         final MultivariateCompositionMapper termF =
  669.                         new MultivariateCompositionMapper(lowerTerm.getCoeff(), fDSIndex, productIndicesF);
  670.         termF.sort();
  671.         return termF;

  673.     }

  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) {

  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);

  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;

  698.         // create new term
  699.         final MultivariateCompositionMapper termP =
  700.                         new MultivariateCompositionMapper(lowerTerm.getCoeff(), lowerTerm.dsIndex, productIndicesP);
  701.         termP.sort();
  702.         return termP;

  704.     }

  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 {

  741.         // safety check
  742.         MathUtils.checkDimension(orders.length, getFreeParameters());
  743.         return getPartialDerivativeIndex(parameters, order, sizes, orders);

  745.     }

  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 {

  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) {

  768.             // derivative order for current free parameter
  769.             int derivativeOrder = orders[i];

  771.             // safety check
  772.             ordersSum += derivativeOrder;
  773.             if (ordersSum > order) {
  774.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
  775.                                                        ordersSum, order);
  776.             }

  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.             }

  786.         }

  788.         return index;

  790.     }

  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.     }

  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.     }

  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.     }

  840.     /** Get the number of free parameters.
  841.      * @return number of free parameters
  842.      */
  843.     public int getFreeParameters() {
  844.         return parameters;
  845.     }

  847.     /** Get the derivation order.
  848.      * @return derivation order
  849.      */
  850.     public int getOrder() {
  851.         return order;
  852.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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) {

  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]);

  1346.         // set up value
  1347.         result[resultOffset] = rem;

  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.         }

  1354.     }

  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) {

  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());

  1374.         // set up value
  1375.         result[resultOffset] = rem;

  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.         }

  1382.     }

  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) {

  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.         }


  1420.         // apply function composition
  1421.         compose(operand, operandOffset, function, result, resultOffset);

  1423.     }

  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) {

  1439.         final T zero = operand[operandOffset].getField().getZero();

  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.         }


  1464.         // apply function composition
  1465.         compose(operand, operandOffset, function, result, resultOffset);

  1467.     }

  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) {

  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.         }

  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.         }

  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.         }

  1509.         // apply function composition
  1510.         compose(operand, operandOffset, function, result, resultOffset);

  1512.     }

  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) {

  1527.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  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.         }

  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.         }

  1557.         // apply function composition
  1558.         compose(operand, operandOffset, function, result, resultOffset);

  1560.     }

  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) {

  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.         }

  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];

  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.         }

  1604.         double coefficient = n;
  1605.         for (int i = 1; i <= order; ++i) {
  1606.             function[i] *= coefficient;
  1607.             coefficient *= n - i;
  1608.         }

  1610.         // apply function composition
  1611.         compose(operand, operandOffset, function, result, resultOffset);

  1613.     }

  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) {

  1628.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  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);

  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.         }

  1660.         double coefficient = n;
  1661.         for (int i = 1; i <= order; ++i) {
  1662.             function[i]  = function[i].multiply(coefficient);
  1663.             coefficient *= n - i;
  1664.         }

  1666.         // apply function composition
  1667.         compose(operand, operandOffset, function, result, resultOffset);

  1669.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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) {

  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.         }

  1796.         // apply function composition
  1797.         compose(operand, operandOffset, function, result, resultOffset);

  1799.     }

  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) {

  1814.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  1837.         // apply function composition
  1838.         compose(operand, operandOffset, function, result, resultOffset);

  1840.     }

  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) {

  1853.         // create the function value and derivatives
  1854.         double[] function = new double[1 + order];
  1855.         Arrays.fill(function, FastMath.exp(operand[operandOffset]));

  1857.         // apply function composition
  1858.         compose(operand, operandOffset, function, result, resultOffset);

  1860.     }

  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) {

  1874.         final Field<T> field = operand[operandOffset].getField();

  1876.         // create the function value and derivatives
  1877.         T[] function = MathArrays.buildArray(field, 1 + order);
  1878.         Arrays.fill(function, operand[operandOffset].exp());

  1880.         // apply function composition
  1881.         compose(operand, operandOffset, function, result, resultOffset);

  1883.     }

  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) {

  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]));

  1901.         // apply function composition
  1902.         compose(operand, operandOffset, function, result, resultOffset);

  1904.     }

  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) {

  1918.         final Field<T> field = operand[operandOffset].getField();

  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());

  1925.         // apply function composition
  1926.         compose(operand, operandOffset, function, result, resultOffset);

  1928.     }

  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) {

  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.         }

  1953.         // apply function composition
  1954.         compose(operand, operandOffset, function, result, resultOffset);

  1956.     }

  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) {

  1970.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  1984.         // apply function composition
  1985.         compose(operand, operandOffset, function, result, resultOffset);

  1987.     }

  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) {

  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.         }

  2011.         // apply function composition
  2012.         compose(operand, operandOffset, function, result, resultOffset);

  2014.     }

  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) {

  2027.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  2041.         // apply function composition
  2042.         compose(operand, operandOffset, function, result, resultOffset);

  2044.     }

  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) {

  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.         }

  2068.         // apply function composition
  2069.         compose(operand, operandOffset, function, result, resultOffset);

  2071.     }

  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) {

  2084.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  2098.         // apply function composition
  2099.         compose(operand, operandOffset, function, result, resultOffset);

  2101.     }

  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) {

  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.         }

  2125.         // apply function composition
  2126.         compose(operand, operandOffset, function, result, resultOffset);

  2128.     }

  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) {

  2142.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  2161.         // apply function composition
  2162.         compose(operand, operandOffset, function, result, resultOffset);

  2164.     }

  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) {

  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.         }

  2188.         // apply function composition
  2189.         compose(operand, operandOffset, function, result, resultOffset);

  2191.     }

  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) {

  2205.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  2224.         // apply function composition
  2225.         compose(operand, operandOffset, function, result, resultOffset);

  2227.     }

  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) {

  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.         }

  2261.         // apply function composition
  2262.         compose(operand, operandOffset, functionSin, sin, sinOffset);
  2263.         compose(operand, operandOffset, functionCos, cos, cosOffset);

  2265.     }

  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) {

  2285.         final Field<T> field = operand[operandOffset].getField();

  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);

  2307.         // apply function composition
  2308.         compose(operand, operandOffset, functionSin, sin, sinOffset);
  2309.         compose(operand, operandOffset, functionCos, cos, cosOffset);

  2311.     }

  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) {

  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;

  2329.         if (order > 0) {

  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) {

  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.                 }

  2358.                 function[n] = v;

  2360.             }
  2361.         }

  2363.         // apply function composition
  2364.         compose(operand, operandOffset, function, result, resultOffset);

  2366.     }

  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) {

  2380.         final Field<T> field = operand[operandOffset].getField();

  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;

  2387.         if (order > 0) {

  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) {

  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.                 }

  2416.                 function[n] = v;

  2418.             }
  2419.         }

  2421.         // apply function composition
  2422.         compose(operand, operandOffset, function, result, resultOffset);

  2424.     }

  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) {

  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) {

  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.                 }

  2472.                 coeff *= f;
  2473.                 function[n] = coeff * v;

  2475.             }
  2476.         }

  2478.         // apply function composition
  2479.         compose(operand, operandOffset, function, result, resultOffset);

  2481.     }

  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) {

  2495.         final Field<T> field = operand[operandOffset].getField();

  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) {

  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.                 }

  2532.                 coeff = coeff.multiply(f);
  2533.                 function[n] = coeff.multiply(v);

  2535.             }
  2536.         }

  2538.         // apply function composition
  2539.         compose(operand, operandOffset, function, result, resultOffset);

  2541.     }

  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) {

  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) {

  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.                 }

  2589.                 coeff *= f;
  2590.                 function[n] = coeff * v;

  2592.             }
  2593.         }

  2595.         // apply function composition
  2596.         compose(operand, operandOffset, function, result, resultOffset);

  2598.     }

  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) {

  2612.         final Field<T> field = operand[operandOffset].getField();

  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) {

  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.                 }

  2649.                 coeff = coeff.multiply(f);
  2650.                 function[n] = coeff.multiply(v);

  2652.             }
  2653.         }

  2655.         // apply function composition
  2656.         compose(operand, operandOffset, function, result, resultOffset);

  2658.     }

  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) {

  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) {

  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.                 }

  2706.                 coeff *= f;
  2707.                 function[n] = coeff * v;

  2709.             }
  2710.         }

  2712.         // apply function composition
  2713.         compose(operand, operandOffset, function, result, resultOffset);

  2715.     }

  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) {

  2729.         final Field<T> field = operand[operandOffset].getField();

  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) {

  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.                 }

  2766.                 coeff = coeff.multiply(f);
  2767.                 function[n] = coeff.multiply(v);

  2769.             }
  2770.         }

  2772.         // apply function composition
  2773.         compose(operand, operandOffset, function, result, resultOffset);

  2775.     }

  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) {

  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)

  2799.         if (x[xOffset] >= 0) {

  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.             }

  2809.         } else {

  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.             }

  2821.         }

  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]);

  2826.     }

  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) {

  2843.         final Field<T> field = y[yOffset].getField();

  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)

  2853.         if (x[xOffset].getReal() >= 0) {

  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.             }

  2863.         } else {

  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.             }

  2875.         }

  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]);

  2880.     }

  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) {

  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.         }

  2903.         // apply function composition
  2904.         compose(operand, operandOffset, function, result, resultOffset);

  2906.     }

  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) {

  2920.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  2932.         // apply function composition
  2933.         compose(operand, operandOffset, function, result, resultOffset);

  2935.     }

  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) {

  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.         }

  2958.         // apply function composition
  2959.         compose(operand, operandOffset, function, result, resultOffset);

  2961.     }

  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) {

  2975.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  2987.         // apply function composition
  2988.         compose(operand, operandOffset, function, result, resultOffset);

  2990.     }

  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) {

  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.         }

  3023.         // apply function composition
  3024.         compose(operand, operandOffset, functionSinh, sinh, sinhOffset);
  3025.         compose(operand, operandOffset, functionCosh, cosh, coshOffset);

  3027.     }

  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) {

  3047.         final Field<T> field = operand[operandOffset].getField();

  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.         }

  3060.         // apply function composition
  3061.         compose(operand, operandOffset, functionSinh, sinh, sinhOffset);
  3062.         compose(operand, operandOffset, functionCosh, cosh, coshOffset);

  3064.     }

  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) {

  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;

  3082.         if (order > 0) {

  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) {

  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.                 }

  3111.                 function[n] = v;

  3113.             }
  3114.         }

  3116.         // apply function composition
  3117.         compose(operand, operandOffset, function, result, resultOffset);

  3119.     }

  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) {

  3133.         final Field<T> field = operand[operandOffset].getField();

  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;

  3140.         if (order > 0) {

  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) {

  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.                 }

  3169.                 function[n] = v;

  3171.             }
  3172.         }

  3174.         // apply function composition
  3175.         compose(operand, operandOffset, function, result, resultOffset);

  3177.     }

  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) {

  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) {

  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.                 }

  3225.                 coeff *= f;
  3226.                 function[n] = coeff * v;

  3228.             }
  3229.         }

  3231.         // apply function composition
  3232.         compose(operand, operandOffset, function, result, resultOffset);

  3234.     }

  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) {

  3248.         final Field<T> field = operand[operandOffset].getField();

  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) {

  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.                 }

  3285.                 coeff = coeff.multiply(f);
  3286.                 function[n] = coeff.multiply(v);

  3288.             }
  3289.         }

  3291.         // apply function composition
  3292.         compose(operand, operandOffset, function, result, resultOffset);

  3294.     }

  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) {

  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) {

  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.                 }

  3342.                 coeff *= f;
  3343.                 function[n] = coeff * v;

  3345.             }
  3346.         }

  3348.         // apply function composition
  3349.         compose(operand, operandOffset, function, result, resultOffset);

  3351.     }

  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) {

  3365.         final Field<T> field = operand[operandOffset].getField();

  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) {

  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.                 }

  3402.                 coeff = coeff.multiply(f);
  3403.                 function[n] = coeff.multiply(v);

  3405.             }
  3406.         }

  3408.         // apply function composition
  3409.         compose(operand, operandOffset, function, result, resultOffset);

  3411.     }

  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) {

  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) {

  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.                 }

  3459.                 coeff *= f;
  3460.                 function[n] = coeff * v;

  3462.             }
  3463.         }

  3465.         // apply function composition
  3466.         compose(operand, operandOffset, function, result, resultOffset);

  3468.     }

  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) {

  3482.         final Field<T> field = operand[operandOffset].getField();

  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) {

  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.                 }

  3519.                 coeff = coeff.multiply(f);
  3520.                 function[n] = coeff.multiply(v);

  3522.             }
  3523.         }

  3525.         // apply function composition
  3526.         compose(operand, operandOffset, function, result, resultOffset);

  3528.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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) {

  3766.         final List<T> combined = new ArrayList<>(terms.size());

  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.         }

  3784.         return combined.toArray((T[]) Array.newInstance(terms.get(0).getClass(), combined.size()));

  3786.     }

  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>> {

  3794.         /** Multiplication coefficient. */
  3795.         private int coeff;

  3797.         /** Simple constructor.
  3798.          * @param coeff multiplication coefficient
  3799.          */
  3800.         AbstractMapper(final int coeff) {
  3801.             this.coeff    = coeff;
  3802.         }

  3804.         /** Set the multiplication coefficient.
  3805.          * @param coeff new coefficient
  3806.          */
  3807.         public void setCoeff(final int coeff) {
  3808.             this.coeff = coeff;
  3809.         }

  3811.         /** Get the multiplication coefficient.
  3812.          * @return multiplication coefficient
  3813.          */
  3814.         public int getCoeff() {
  3815.             return coeff;
  3816.         }

  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);

  3824.     }

  3826.     /** Multiplication mapper.
  3827.      * @since 2.2
  3828.      */
  3829.     private static class MultiplicationMapper extends AbstractMapper<MultiplicationMapper> {

  3831.         /** Left hand side index. */
  3832.         private final int lhsIndex;

  3834.         /** Right hand side index. */
  3835.         private final int rhsIndex;

  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.         }

  3848.         /** {@inheritDoc} */
  3849.         @Override
  3850.         public boolean isSimilar(final MultiplicationMapper other) {
  3851.             return lhsIndex == other.lhsIndex && rhsIndex == other.rhsIndex;
  3852.         }

  3854.     }

  3856.     /** Univariate composition mapper.
  3857.      * @since 2.2
  3858.      */
  3859.     private static class UnivariateCompositionMapper extends AbstractMapper<UnivariateCompositionMapper> {

  3861.         /** Univariate derivative index. */
  3862.         private final int fIndex;

  3864.         /** Derivative structure indices. */
  3865.         private final int[] dsIndices;

  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.         }

  3878.         /** Sort the derivatives structures indices.
  3879.          */
  3880.         public void sort() {
  3881.             Arrays.sort(dsIndices);
  3882.         }

  3884.         /** {@inheritDoc} */
  3885.         @Override
  3886.         public boolean isSimilar(final UnivariateCompositionMapper other) {

  3888.             if (fIndex == other.fIndex && dsIndices.length == other.dsIndices.length) {

  3890.                 for (int j = 0; j < dsIndices.length; ++j) {
  3891.                     if (dsIndices[j] != other.dsIndices[j]) {
  3892.                         return false;
  3893.                     }
  3894.                 }

  3896.                 return true;

  3898.             }

  3900.             return false;

  3902.         }

  3904.     }

  3906.     /** Multivariate composition mapper.
  3907.      * @since 2.2
  3908.      */
  3909.     private static class MultivariateCompositionMapper extends AbstractMapper<MultivariateCompositionMapper> {

  3911.         /** Multivariate derivative index. */
  3912.         private final int dsIndex;

  3914.         /** Indices of the intermediate variables derivatives products. */
  3915.         private final int[] productIndices;

  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.         }

  3928.         /** Sort the indices of the intermediate variables derivatives products.
  3929.          */
  3930.         public void sort() {
  3931.             Arrays.sort(productIndices);
  3932.         }

  3934.         /** {@inheritDoc} */
  3935.         @Override
  3936.         public boolean isSimilar(final MultivariateCompositionMapper other) {

  3938.             if (dsIndex == other.dsIndex && productIndices.length == other.productIndices.length) {

  3940.                 for (int j = 0; j < productIndices.length; ++j) {
  3941.                     if (productIndices[j] != other.productIndices[j]) {
  3942.                         return false;
  3943.                     }
  3944.                 }

  3946.                 return true;

  3948.             }

  3950.             return false;

  3952.         }

  3954.     }

  3956. }
Created with JaCoCo 0.8.12.202403310830