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