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