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    }