View Javadoc
1   /*
2    * Licensed to the Hipparchus project 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 Hipparchus project 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  package org.hipparchus.analysis.differentiation;
19  
20  import org.hipparchus.CalculusFieldElement;
21  import org.hipparchus.CalculusFieldElementAbstractTest;
22  import org.hipparchus.Field;
23  import org.hipparchus.analysis.CalculusFieldMultivariateFunction;
24  import org.hipparchus.analysis.CalculusFieldMultivariateVectorFunction;
25  import org.hipparchus.analysis.polynomials.FieldPolynomialFunction;
26  import org.hipparchus.analysis.polynomials.PolynomialFunction;
27  import org.hipparchus.dfp.Dfp;
28  import org.hipparchus.dfp.DfpField;
29  import org.hipparchus.exception.LocalizedCoreFormats;
30  import org.hipparchus.exception.MathIllegalArgumentException;
31  import org.hipparchus.random.RandomGenerator;
32  import org.hipparchus.random.Well1024a;
33  import org.hipparchus.random.Well19937a;
34  import org.hipparchus.util.ArithmeticUtils;
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.Precision;
41  import org.junit.jupiter.api.Assertions;
42  import org.junit.jupiter.api.Test;
43  
44  import java.lang.reflect.Array;
45  import java.util.ArrayList;
46  import java.util.Arrays;
47  import java.util.HashMap;
48  import java.util.List;
49  import java.util.Map;
50  
51  import static org.junit.jupiter.api.Assertions.assertEquals;
52  import static org.junit.jupiter.api.Assertions.assertNotEquals;
53  import static org.junit.jupiter.api.Assertions.assertSame;
54  import static org.junit.jupiter.api.Assertions.assertThrows;
55  import static org.junit.jupiter.api.Assertions.assertTrue;
56  import static org.junit.jupiter.api.Assertions.fail;
57  
58  /**
59   * Abstract test for class {@link FieldDerivativeStructure}.
60   */
61  public abstract class FieldDerivativeStructureAbstractTest<T extends CalculusFieldElement<T>>
62      extends CalculusFieldElementAbstractTest<FieldDerivativeStructure<T>> {
63  
64      protected abstract Field<T> getField();
65  
66      protected T buildScalar(double value) {
67          return getField().getZero().newInstance(value);
68      }
69  
70      protected FDSFactory<T> buildFactory(int parameters, int order) {
71          return new FDSFactory<>(getField(), parameters, order);
72      }
73  
74      @Override
75      protected FieldDerivativeStructure<T> build(final double x) {
76          return buildFactory(2, 1).variable(0, x);
77      }
78  
79      @Test
80      public void testWrongFieldVariableIndex() {
81          assertThrows(MathIllegalArgumentException.class, () -> {
82              buildFactory(3, 1).variable(3, buildScalar(1.0));
83          });
84      }
85  
86      @Test
87      public void testWrongPrimitiveVariableIndex() {
88          assertThrows(MathIllegalArgumentException.class, () -> {
89              final FDSFactory<T> factory = buildFactory(3, 1);
90              factory.variable(3, 1.0);
91          });
92      }
93  
94      @Test
95      public void testMissingOrders() {
96          assertThrows(MathIllegalArgumentException.class, () -> {
97              final FDSFactory<T> factory = buildFactory(3, 1);
98              factory.variable(0, 1.0).getPartialDerivative(0, 1);
99          });
100     }
101 
102     @Test
103     public void testWrongDimensionField() {
104         assertThrows(MathIllegalArgumentException.class, () -> {
105             final FDSFactory<T> factory = buildFactory(3, 1);
106             factory.build(buildScalar(1.0), buildScalar(1.0), buildScalar(1.0), buildScalar(1.0), buildScalar(1.0));
107         });
108     }
109 
110     @Test
111     public void testWrongDimensionPrimitive() {
112         assertThrows(MathIllegalArgumentException.class, () -> {
113             final FDSFactory<T> factory = buildFactory(3, 1);
114             factory.build(1.0, 1.0, 1.0, 1.0, 1.0);
115         });
116     }
117 
118     @Test
119     public void testTooLargeOrder() {
120         assertThrows(MathIllegalArgumentException.class, () -> {
121             final FDSFactory<T> factory = buildFactory(3, 1);
122             factory.variable(0, 1.0).getPartialDerivative(1, 1, 2);
123         });
124     }
125 
126     @Test
127     public void testVariableWithoutDerivativeField() {
128         final FDSFactory<T> factory = buildFactory(1, 0);
129         FieldDerivativeStructure<T> v = factory.variable(0, buildScalar(1.0));
130         assertEquals(1.0, v.getReal(), 1.0e-15);
131     }
132 
133     @Test
134     public void testVariableWithoutDerivativePrimitive() {
135         final FDSFactory<T> factory = buildFactory(1, 0);
136         FieldDerivativeStructure<T> v = factory.variable(0, 1.0);
137         assertEquals(1.0, v.getReal(), 1.0e-15);
138     }
139 
140     @Test
141     public void testVariableWithoutDerivative1() {
142         assertThrows(MathIllegalArgumentException.class, () -> {
143             final FDSFactory<T> factory = buildFactory(1, 0);
144             FieldDerivativeStructure<T> v = factory.variable(0, 1.0);
145             assertEquals(1.0, v.getPartialDerivative(1).getReal(), 1.0e-15);
146         });
147     }
148 
149     @Test
150     public void testVariable() {
151         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
152             final FDSFactory<T> factory = buildFactory(3, maxOrder);
153             checkF0F1(factory.variable(0, 1.0),
154                       1.0, 1.0, 0.0, 0.0);
155             checkF0F1(factory.variable(1, 2.0),
156                       2.0, 0.0, 1.0, 0.0);
157             checkF0F1(factory.variable(2, 3.0),
158                       3.0, 0.0, 0.0, 1.0);
159         }
160     }
161 
162     @Test
163     public void testConstant() {
164         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
165             final FDSFactory<T> factory = buildFactory(3, maxOrder);
166             checkF0F1(factory.constant(FastMath.PI),
167                       FastMath.PI, 0.0, 0.0, 0.0);
168         }
169     }
170 
171     @Test
172     public void testFieldAdd() {
173         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
174             final FDSFactory<T> factory = buildFactory(3, maxOrder);
175             checkF0F1(factory.variable(0, 1.0).add(buildScalar(5)), 6.0, 1.0, 0.0, 0.0);
176             checkF0F1(factory.variable(1, 2.0).add(buildScalar(5)), 7.0, 0.0, 1.0, 0.0);
177             checkF0F1(factory.variable(2, 3.0).add(buildScalar(5)), 8.0, 0.0, 0.0, 1.0);
178         }
179     }
180 
181     @Test
182     public void testPrimitiveAdd() {
183         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
184             final FDSFactory<T> factory = buildFactory(3, maxOrder);
185             checkF0F1(factory.variable(0, 1.0).add(5), 6.0, 1.0, 0.0, 0.0);
186             checkF0F1(factory.variable(1, 2.0).add(5), 7.0, 0.0, 1.0, 0.0);
187             checkF0F1(factory.variable(2, 3.0).add(5), 8.0, 0.0, 0.0, 1.0);
188         }
189     }
190 
191     @Test
192     public void testAdd() {
193         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
194             final FDSFactory<T> factory = buildFactory(3, maxOrder);
195             FieldDerivativeStructure<T> x = factory.variable(0, 1.0);
196             FieldDerivativeStructure<T> y = factory.variable(1, 2.0);
197             FieldDerivativeStructure<T> z = factory.variable(2, 3.0);
198             FieldDerivativeStructure<T> xyz = x.add(y.add(z));
199             checkF0F1(xyz, x.getValue().getReal() + y.getValue().getReal() + z.getValue().getReal(), 1.0, 1.0, 1.0);
200         }
201     }
202 
203     @Test
204     public void testFieldSubtract() {
205         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
206             final FDSFactory<T> factory = buildFactory(3, maxOrder);
207             checkF0F1(factory.variable(0, 1.0).subtract(buildScalar(5)), -4.0, 1.0, 0.0, 0.0);
208             checkF0F1(factory.variable(1, 2.0).subtract(buildScalar(5)), -3.0, 0.0, 1.0, 0.0);
209             checkF0F1(factory.variable(2, 3.0).subtract(buildScalar(5)), -2.0, 0.0, 0.0, 1.0);
210         }
211     }
212 
213     @Test
214     public void testPrimitiveSubtract() {
215         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
216             final FDSFactory<T> factory = buildFactory(3, maxOrder);
217             checkF0F1(factory.variable(0, 1.0).subtract(5), -4.0, 1.0, 0.0, 0.0);
218             checkF0F1(factory.variable(1, 2.0).subtract(5), -3.0, 0.0, 1.0, 0.0);
219             checkF0F1(factory.variable(2, 3.0).subtract(5), -2.0, 0.0, 0.0, 1.0);
220         }
221     }
222 
223     @Test
224     public void testSubtract() {
225         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
226             final FDSFactory<T> factory = buildFactory(3, maxOrder);
227             FieldDerivativeStructure<T> x = factory.variable(0, 1.0);
228             FieldDerivativeStructure<T> y = factory.variable(1, 2.0);
229             FieldDerivativeStructure<T> z = factory.variable(2, 3.0);
230             FieldDerivativeStructure<T> xyz = x.subtract(y.subtract(z));
231             checkF0F1(xyz, x.getReal() - (y.getReal() - z.getReal()), 1.0, -1.0, 1.0);
232         }
233     }
234 
235     @Test
236     public void testFieldMultiply() {
237         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
238             final FDSFactory<T> factory = buildFactory(3, maxOrder);
239             checkF0F1(factory.variable(0, 1.0).multiply(buildScalar(5)),  5.0, 5.0, 0.0, 0.0);
240             checkF0F1(factory.variable(1, 2.0).multiply(buildScalar(5)), 10.0, 0.0, 5.0, 0.0);
241             checkF0F1(factory.variable(2, 3.0).multiply(buildScalar(5)), 15.0, 0.0, 0.0, 5.0);
242         }
243     }
244 
245     @Test
246     public void testPrimitiveMultiply() {
247         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
248             final FDSFactory<T> factory = buildFactory(3, maxOrder);
249             checkF0F1(factory.variable(0, 1.0).multiply(5),  5.0, 5.0, 0.0, 0.0);
250             checkF0F1(factory.variable(1, 2.0).multiply(5), 10.0, 0.0, 5.0, 0.0);
251             checkF0F1(factory.variable(2, 3.0).multiply(5), 15.0, 0.0, 0.0, 5.0);
252         }
253     }
254 
255     @Test
256     public void testMultiply() {
257         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
258             final FDSFactory<T> factory = buildFactory(3, maxOrder);
259             FieldDerivativeStructure<T> x = factory.variable(0, 1.0);
260             FieldDerivativeStructure<T> y = factory.variable(1, 2.0);
261             FieldDerivativeStructure<T> z = factory.variable(2, 3.0);
262             FieldDerivativeStructure<T> xyz = x.multiply(y.multiply(z));
263             for (int i = 0; i <= maxOrder; ++i) {
264                 for (int j = 0; j <= maxOrder; ++j) {
265                     for (int k = 0; k <= maxOrder; ++k) {
266                         if (i + j + k <= maxOrder) {
267                             assertEquals((i == 0 ? x.getReal() : (i == 1 ? 1.0 : 0.0)) *
268                                                 (j == 0 ? y.getReal() : (j == 1 ? 1.0 : 0.0)) *
269                                                 (k == 0 ? z.getReal() : (k == 1 ? 1.0 : 0.0)),
270                                                 xyz.getPartialDerivative(i, j, k).getReal(),
271                                                 1.0e-15);
272                         }
273                     }
274                 }
275             }
276         }
277     }
278 
279     @Test
280     public void testFieldDivide() {
281         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
282             final FDSFactory<T> factory = buildFactory(3, maxOrder);
283             checkF0F1(factory.variable(0, 1.0).divide(buildScalar(2)),  0.5, 0.5, 0.0, 0.0);
284             checkF0F1(factory.variable(1, 2.0).divide(buildScalar(2)),  1.0, 0.0, 0.5, 0.0);
285             checkF0F1(factory.variable(2, 3.0).divide(buildScalar(2)),  1.5, 0.0, 0.0, 0.5);
286         }
287     }
288 
289     @Test
290     public void testPrimitiveDivide() {
291         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
292             final FDSFactory<T> factory = buildFactory(3, maxOrder);
293             checkF0F1(factory.variable(0, 1.0).divide(2),  0.5, 0.5, 0.0, 0.0);
294             checkF0F1(factory.variable(1, 2.0).divide(2),  1.0, 0.0, 0.5, 0.0);
295             checkF0F1(factory.variable(2, 3.0).divide(2),  1.5, 0.0, 0.0, 0.5);
296         }
297     }
298 
299     @Test
300     public void testNegate() {
301         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
302             final FDSFactory<T> factory = buildFactory(3, maxOrder);
303             checkF0F1(factory.variable(0, 1.0).negate(), -1.0, -1.0, 0.0, 0.0);
304             checkF0F1(factory.variable(1, 2.0).negate(), -2.0, 0.0, -1.0, 0.0);
305             checkF0F1(factory.variable(2, 3.0).negate(), -3.0, 0.0, 0.0, -1.0);
306         }
307     }
308 
309     @Test
310     public void testReciprocal() {
311         final FDSFactory<T> factory = buildFactory(1, 6);
312         for (double x = 0.1; x < 1.2; x += 0.1) {
313             FieldDerivativeStructure<T> r = factory.variable(0, x).reciprocal();
314             assertEquals(1 / x, r.getReal(), 1.0e-15);
315             for (int i = 1; i < r.getOrder(); ++i) {
316                 double expected = ArithmeticUtils.pow(-1, i) * CombinatoricsUtils.factorial(i) /
317                                   FastMath.pow(x, i + 1);
318                 assertEquals(expected, r.getPartialDerivative(i).getReal(), 1.0e-15 * FastMath.abs(expected));
319             }
320         }
321     }
322 
323     @Test
324     public void testPow() {
325         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
326             final FDSFactory<T> factory = buildFactory(3, maxOrder);
327             for (int n = 0; n < 10; ++n) {
328 
329                 FieldDerivativeStructure<T> x = factory.variable(0, 1.0);
330                 FieldDerivativeStructure<T> y = factory.variable(1, 2.0);
331                 FieldDerivativeStructure<T> z = factory.variable(2, 3.0);
332                 List<FieldDerivativeStructure<T>> list = Arrays.asList(x, y, z,
333                                                                                x.add(y).add(z),
334                                                                                x.multiply(y).multiply(z));
335 
336                 if (n == 0) {
337                     for (FieldDerivativeStructure<T> ds : list) {
338                         checkEquals(ds.getField().getOne(), ds.pow(n), 1.0e-15);
339                     }
340                 } else if (n == 1) {
341                     for (FieldDerivativeStructure<T> ds : list) {
342                         checkEquals(ds, ds.pow(n), 1.0e-15);
343                     }
344                 } else {
345                     for (FieldDerivativeStructure<T> ds : list) {
346                         FieldDerivativeStructure<T> p = ds.getField().getOne();
347                         for (int i = 0; i < n; ++i) {
348                             p = p.multiply(ds);
349                         }
350                         checkEquals(p, ds.pow(n), 1.0e-15);
351                     }
352                 }
353             }
354         }
355     }
356 
357     @Test
358     public void testPowDoubleDS() {
359         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
360 
361             final FDSFactory<T> factory = buildFactory(3, maxOrder);
362             FieldDerivativeStructure<T> x = factory.variable(0, 0.1);
363             FieldDerivativeStructure<T> y = factory.variable(1, 0.2);
364             FieldDerivativeStructure<T> z = factory.variable(2, 0.3);
365             List<FieldDerivativeStructure<T>> list = Arrays.asList(x, y, z,
366                                                                            x.add(y).add(z),
367                                                                            x.multiply(y).multiply(z));
368 
369             for (FieldDerivativeStructure<T> ds : list) {
370                 // the special case a = 0 is included here
371                 for (double a : new double[] { 0.0, 0.1, 1.0, 2.0, 5.0 }) {
372                     FieldDerivativeStructure<T> reference = (a == 0) ?
373                                                     x.getField().getZero() :
374                                                     factory.constant(a).pow(ds);
375                     FieldDerivativeStructure<T> result = FieldDerivativeStructure.pow(a, ds);
376                     checkEquals(reference, result, 2.0e-14 * FastMath.abs(reference.getReal()));
377                 }
378 
379             }
380 
381             // negative base: -1^x can be evaluated for integers only, so value is sometimes OK, derivatives are always NaN
382             FieldDerivativeStructure<T> negEvenInteger = FieldDerivativeStructure.pow(-2.0, factory.variable(0, 2.0));
383             assertEquals(4.0, negEvenInteger.getReal(), 1.0e-15);
384             assertTrue(Double.isNaN(negEvenInteger.getPartialDerivative(1, 0, 0).getReal()));
385             FieldDerivativeStructure<T> negOddInteger = FieldDerivativeStructure.pow(-2.0, factory.variable(0, 3.0));
386             assertEquals(-8.0, negOddInteger.getReal(), 1.0e-15);
387             assertTrue(Double.isNaN(negOddInteger.getPartialDerivative(1, 0, 0).getReal()));
388             FieldDerivativeStructure<T> negNonInteger = FieldDerivativeStructure.pow(-2.0, factory.variable(0, 2.001));
389             assertTrue(Double.isNaN(negNonInteger.getReal()));
390             assertTrue(Double.isNaN(negNonInteger.getPartialDerivative(1, 0, 0).getReal()));
391 
392             FieldDerivativeStructure<T> zeroNeg = FieldDerivativeStructure.pow(0.0, factory.variable(0, -1.0));
393             assertTrue(Double.isNaN(zeroNeg.getReal()));
394             assertTrue(Double.isNaN(zeroNeg.getPartialDerivative(1, 0, 0).getReal()));
395             FieldDerivativeStructure<T> posNeg = FieldDerivativeStructure.pow(2.0, factory.variable(0, -2.0));
396             assertEquals(1.0 / 4.0, posNeg.getReal(), 1.0e-15);
397             assertEquals(FastMath.log(2.0) / 4.0, posNeg.getPartialDerivative(1, 0, 0).getReal(), 1.0e-15);
398 
399             // very special case: a = 0 and power = 0
400             FieldDerivativeStructure<T> zeroZero = FieldDerivativeStructure.pow(0.0, factory.variable(0, 0.0));
401 
402             // this should be OK for simple first derivative with one variable only ...
403             assertEquals(1.0, zeroZero.getReal(), 1.0e-15);
404             assertEquals(Double.NEGATIVE_INFINITY, zeroZero.getPartialDerivative(1, 0, 0).getReal(), 1.0e-15);
405 
406             // the following checks show a LIMITATION of the current implementation
407             // we have no way to tell x is a pure linear variable x = 0
408             // we only say: "x is a structure with value = 0.0,
409             // first derivative with respect to x = 1.0, and all other derivatives
410             // (first order with respect to y and z and higher derivatives) all 0.0.
411             // We have function f(x) = a^x and x = 0 so we compute:
412             // f(0) = 1, f'(0) = ln(a), f''(0) = ln(a)^2. The limit of these values
413             // when a converges to 0 implies all derivatives keep switching between
414             // +infinity and -infinity.
415             //
416             // Function composition rule for first derivatives is:
417             // d[f(g(x,y,z))]/dy = f'(g(x,y,z)) * dg(x,y,z)/dy
418             // so given that in our case x represents g and does not depend
419             // on y or z, we have dg(x,y,z)/dy = 0
420             // applying the composition rules gives:
421             // d[f(g(x,y,z))]/dy = f'(g(x,y,z)) * dg(x,y,z)/dy
422             //                 = -infinity * 0
423             //                 = NaN
424             // if we knew x is really the x variable and not the identity
425             // function applied to x, we would not have computed f'(g(x,y,z)) * dg(x,y,z)/dy
426             // and we would have found that the result was 0 and not NaN
427             assertTrue(Double.isNaN(zeroZero.getPartialDerivative(0, 1, 0).getReal()));
428             assertTrue(Double.isNaN(zeroZero.getPartialDerivative(0, 0, 1).getReal()));
429 
430             // Function composition rule for second derivatives is:
431             // d2[f(g(x))]/dx2 = f''(g(x)) * [g'(x)]^2 + f'(g(x)) * g''(x)
432             // when function f is the a^x root and x = 0 we have:
433             // f(0) = 1, f'(0) = ln(a), f''(0) = ln(a)^2 which for a = 0 implies
434             // all derivatives keep switching between +infinity and -infinity
435             // so given that in our case x represents g, we have g(x) = 0,
436             // g'(x) = 1 and g''(x) = 0
437             // applying the composition rules gives:
438             // d2[f(g(x))]/dx2 = f''(g(x)) * [g'(x)]^2 + f'(g(x)) * g''(x)
439             //                 = +infinity * 1^2 + -infinity * 0
440             //                 = +infinity + NaN
441             //                 = NaN
442             // if we knew x is really the x variable and not the identity
443             // function applied to x, we would not have computed f'(g(x)) * g''(x)
444             // and we would have found that the result was +infinity and not NaN
445             if (maxOrder > 1) {
446                 assertTrue(Double.isNaN(zeroZero.getPartialDerivative(2, 0, 0).getReal()));
447                 assertTrue(Double.isNaN(zeroZero.getPartialDerivative(0, 2, 0).getReal()));
448                 assertTrue(Double.isNaN(zeroZero.getPartialDerivative(0, 0, 2).getReal()));
449                 assertTrue(Double.isNaN(zeroZero.getPartialDerivative(1, 1, 0).getReal()));
450                 assertTrue(Double.isNaN(zeroZero.getPartialDerivative(0, 1, 1).getReal()));
451                 assertTrue(Double.isNaN(zeroZero.getPartialDerivative(1, 1, 0).getReal()));
452             }
453 
454             // very special case: 0^0 where the power is a primitive
455             FieldDerivativeStructure<T> zeroDsZeroDouble = factory.variable(0, 0.0).pow(0.0);
456             boolean first = true;
457             for (final T d : zeroDsZeroDouble.getAllDerivatives()) {
458                 if (first) {
459                     assertEquals(1.0, d.getReal(), Precision.EPSILON);
460                     first = false;
461                 } else {
462                     assertEquals(0.0, d.getReal(), Precision.SAFE_MIN);
463                 }
464             }
465             FieldDerivativeStructure<T> zeroDsZeroInt = factory.variable(0, 0.0).pow(0);
466             first = true;
467             for (final T d : zeroDsZeroInt.getAllDerivatives()) {
468                 if (first) {
469                     assertEquals(1.0, d.getReal(), Precision.EPSILON);
470                     first = false;
471                 } else {
472                     assertEquals(0.0, d.getReal(), Precision.SAFE_MIN);
473                 }
474             }
475 
476             // 0^p with p smaller than 1.0
477             FieldDerivativeStructure<T> u = factory.variable(1, -0.0).pow(0.25);
478             for (int i0 = 0; i0 <= maxOrder; ++i0) {
479                 for (int i1 = 0; i1 <= maxOrder; ++i1) {
480                     for (int i2 = 0; i2 <= maxOrder; ++i2) {
481                         if (i0 + i1 + i2 <= maxOrder) {
482                             assertEquals(0.0, u.getPartialDerivative(i0, i1, i2).getReal(), 1.0e-10);
483                         }
484                     }
485                 }
486             }
487         }
488 
489     }
490 
491     @Test
492     public void testExpression() {
493         final FDSFactory<T> factory = buildFactory(3, 5);
494         double epsilon = 2.5e-13;
495         for (double x = 0; x < 2; x += 0.2) {
496             FieldDerivativeStructure<T> dsX = factory.variable(0, x);
497             for (double y = 0; y < 2; y += 0.2) {
498                 FieldDerivativeStructure<T> dsY = factory.variable(1, y);
499                 for (double z = 0; z >- 2; z -= 0.2) {
500                     FieldDerivativeStructure<T> dsZ = factory.variable(2, z);
501 
502                     // f(x, y, z) = x + 5 x y - 2 z + (8 z x - y)^3
503                     FieldDerivativeStructure<T> ds =
504                             dsX.linearCombination(1, dsX,
505                                                     5, dsX.multiply(dsY),
506                                                     -2, dsZ,
507                                                     1, dsX.linearCombination(8, dsZ.multiply(dsX),
508                                                                                -1, dsY).pow(3));
509                     FieldDerivativeStructure<T> dsOther =
510                                     dsX.linearCombination(1, dsX,
511                                                     5, dsX.multiply(dsY),
512                                                     -2, dsZ).add(dsX.linearCombination   (8, dsZ.multiply(dsX),
513                                                                                          -1, dsY).pow(3));
514                     double f = x + 5 * x * y - 2 * z + FastMath.pow(8 * z * x - y, 3);
515                     assertEquals(f, ds.getReal(),
516                                         FastMath.abs(epsilon * f));
517                     assertEquals(f, dsOther.getReal(),
518                                         FastMath.abs(epsilon * f));
519 
520                     // df/dx = 1 + 5 y + 24 (8 z x - y)^2 z
521                     double dfdx = 1 + 5 * y + 24 * z * FastMath.pow(8 * z * x - y, 2);
522                     assertEquals(dfdx, ds.getPartialDerivative(1, 0, 0).getReal(),
523                                         FastMath.abs(epsilon * dfdx));
524                     assertEquals(dfdx, dsOther.getPartialDerivative(1, 0, 0).getReal(),
525                                         FastMath.abs(epsilon * dfdx));
526 
527                     // df/dxdy = 5 + 48 z*(y - 8 z x)
528                     double dfdxdy = 5 + 48 * z * (y - 8 * z * x);
529                     assertEquals(dfdxdy, ds.getPartialDerivative(1, 1, 0).getReal(),
530                                         FastMath.abs(epsilon * dfdxdy));
531                     assertEquals(dfdxdy, dsOther.getPartialDerivative(1, 1, 0).getReal(),
532                                         FastMath.abs(epsilon * dfdxdy));
533 
534                     // df/dxdydz = 48 (y - 16 z x)
535                     double dfdxdydz = 48 * (y - 16 * z * x);
536                     assertEquals(dfdxdydz, ds.getPartialDerivative(1, 1, 1).getReal(),
537                                         FastMath.abs(epsilon * dfdxdydz));
538                     assertEquals(dfdxdydz, dsOther.getPartialDerivative(1, 1, 1).getReal(),
539                                         FastMath.abs(epsilon * dfdxdydz));
540 
541                 }
542 
543             }
544         }
545     }
546 
547     @Test
548     public void testCompositionOneVariableX() {
549         double epsilon = 1.0e-13;
550         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
551             final FDSFactory<T> factory = buildFactory(1, maxOrder);
552             for (double x = 0.1; x < 1.2; x += 0.1) {
553                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
554                 for (double y = 0.1; y < 1.2; y += 0.1) {
555                     FieldDerivativeStructure<T> dsY = factory.constant(y);
556                     FieldDerivativeStructure<T> f = dsX.divide(dsY).sqrt();
557                     double f0 = FastMath.sqrt(x / y);
558                     assertEquals(f0, f.getReal(), FastMath.abs(epsilon * f0));
559                     if (f.getOrder() > 0) {
560                         double f1 = 1 / (2 * FastMath.sqrt(x * y));
561                         assertEquals(f1, f.getPartialDerivative(1).getReal(), FastMath.abs(epsilon * f1));
562                         if (f.getOrder() > 1) {
563                             double f2 = -f1 / (2 * x);
564                             assertEquals(f2, f.getPartialDerivative(2).getReal(), FastMath.abs(epsilon * f2));
565                             if (f.getOrder() > 2) {
566                                 double f3 = (f0 + x / (2 * y * f0)) / (4 * x * x * x);
567                                 assertEquals(f3, f.getPartialDerivative(3).getReal(), FastMath.abs(epsilon * f3));
568                             }
569                         }
570                     }
571                 }
572             }
573         }
574     }
575 
576     @Test
577     public void testTrigo() {
578         double epsilon = 2.0e-12;
579         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
580             final FDSFactory<T> factory = buildFactory(3, maxOrder);
581             for (double x = 0.1; x < 1.2; x += 0.1) {
582                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
583                 for (double y = 0.1; y < 1.2; y += 0.1) {
584                     FieldDerivativeStructure<T> dsY = factory.variable(1, y);
585                     for (double z = 0.1; z < 1.2; z += 0.1) {
586                         FieldDerivativeStructure<T> dsZ = factory.variable(2, z);
587                         FieldDerivativeStructure<T> f = dsX.divide(dsY.cos().add(dsZ.tan())).sin();
588                         double a = FastMath.cos(y) + FastMath.tan(z);
589                         double f0 = FastMath.sin(x / a);
590                         assertEquals(f0, f.getReal(), FastMath.abs(epsilon * f0));
591                         if (f.getOrder() > 0) {
592                             double dfdx = FastMath.cos(x / a) / a;
593                             assertEquals(dfdx, f.getPartialDerivative(1, 0, 0).getReal(), FastMath.abs(epsilon * dfdx));
594                             double dfdy =  x * FastMath.sin(y) * dfdx / a;
595                             assertEquals(dfdy, f.getPartialDerivative(0, 1, 0).getReal(), FastMath.abs(epsilon * dfdy));
596                             double cz = FastMath.cos(z);
597                             double cz2 = cz * cz;
598                             double dfdz = -x * dfdx / (a * cz2);
599                             assertEquals(dfdz, f.getPartialDerivative(0, 0, 1).getReal(), FastMath.abs(epsilon * dfdz));
600                             if (f.getOrder() > 1) {
601                                 double df2dx2 = -(f0 / (a * a));
602                                 assertEquals(df2dx2, f.getPartialDerivative(2, 0, 0).getReal(), FastMath.abs(epsilon * df2dx2));
603                                 double df2dy2 = x * FastMath.cos(y) * dfdx / a -
604                                                 x * x * FastMath.sin(y) * FastMath.sin(y) * f0 / (a * a * a * a) +
605                                                 2 * FastMath.sin(y) * dfdy / a;
606                                 assertEquals(df2dy2, f.getPartialDerivative(0, 2, 0).getReal(), FastMath.abs(epsilon * df2dy2));
607                                 double c4 = cz2 * cz2;
608                                 double df2dz2 = x * (2 * a * (1 - a * cz * FastMath.sin(z)) * dfdx - x * f0 / a ) / (a * a * a * c4);
609                                 assertEquals(df2dz2, f.getPartialDerivative(0, 0, 2).getReal(), FastMath.abs(epsilon * df2dz2));
610                                 double df2dxdy = dfdy / x  - x * FastMath.sin(y) * f0 / (a * a * a);
611                                 assertEquals(df2dxdy, f.getPartialDerivative(1, 1, 0).getReal(), FastMath.abs(epsilon * df2dxdy));
612                             }
613                         }
614                     }
615                 }
616             }
617         }
618     }
619 
620     @Test
621     public void testSqrtDefinition() {
622         double[] epsilon = new double[] { 5.0e-16, 5.0e-16, 2.7e-15, 5.7e-14, 2.0e-12 };
623         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
624             final FDSFactory<T> factory = buildFactory(1, maxOrder);
625             for (double x = 0.1; x < 1.2; x += 0.001) {
626                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
627                 FieldDerivativeStructure<T> sqrt1 = dsX.pow(0.5);
628                 FieldDerivativeStructure<T> sqrt2 = dsX.sqrt();
629                 FieldDerivativeStructure<T> zero = sqrt1.subtract(sqrt2);
630                 for (int n = 0; n <= maxOrder; ++n) {
631                     assertEquals(0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
632                 }
633             }
634         }
635     }
636 
637     @Test
638     public void testRootNSingularity() {
639         doTestRootNSingularity(true);
640     }
641 
642     protected void doTestRootNSingularity(final boolean signedInfinities) {
643         for (int n = 2; n < 10; ++n) {
644             for (int maxOrder = 0; maxOrder < 12; ++maxOrder) {
645                 final FDSFactory<T> factory = buildFactory(1, maxOrder);
646                 FieldDerivativeStructure<T> dsZero = factory.variable(0, 0.0);
647                 FieldDerivativeStructure<T> rootN  = dsZero.rootN(n);
648                 assertEquals(0.0, rootN.getReal(), 1.0e-20);
649                 if (maxOrder > 0) {
650                     assertTrue(Double.isInfinite(rootN.getPartialDerivative(1).getReal()));
651                     assertTrue(rootN.getPartialDerivative(1).getReal() > 0);
652                     for (int order = 2; order <= maxOrder; ++order) {
653                         // the following checks shows a LIMITATION of the current implementation
654                         // we have no way to tell dsZero is a pure linear variable x = 0
655                         // we only say: "dsZero is a structure with value = 0.0,
656                         // first derivative = 1.0, second and higher derivatives = 0.0".
657                         // Function composition rule for second derivatives is:
658                         // d2[f(g(x))]/dx2 = f''(g(x)) * [g'(x)]^2 + f'(g(x)) * g''(x)
659                         // when function f is the nth root and x = 0 we have:
660                         // f(0) = 0, f'(0) = +infinity, f''(0) = -infinity (and higher
661                         // derivatives keep switching between +infinity and -infinity)
662                         // so given that in our case dsZero represents g, we have g(x) = 0,
663                         // g'(x) = 1 and g''(x) = 0
664                         // applying the composition rules gives:
665                         // d2[f(g(x))]/dx2 = f''(g(x)) * [g'(x)]^2 + f'(g(x)) * g''(x)
666                         //                 = -infinity * 1^2 + +infinity * 0
667                         //                 = -infinity + NaN
668                         //                 = NaN
669                         // if we knew dsZero is really the x variable and not the identity
670                         // function applied to x, we would not have computed f'(g(x)) * g''(x)
671                         // and we would have found that the result was -infinity and not NaN
672                         final double d = rootN.getPartialDerivative(order).getReal();
673                         assertTrue(Double.isNaN(d) || Double.isInfinite(d));
674                     }
675                 }
676 
677                 // the following shows that the limitation explained above is NOT a bug...
678                 // if we set up the higher order derivatives for g appropriately, we do
679                 // compute the higher order derivatives of the composition correctly
680                 double[] gDerivatives = new double[ 1 + maxOrder];
681                 gDerivatives[0] = 0.0;
682                 for (int k = 1; k <= maxOrder; ++k) {
683                     gDerivatives[k] = FastMath.pow(-1.0, k + 1);
684                 }
685                 FieldDerivativeStructure<T> correctRoot = factory.build(gDerivatives).rootN(n);
686                 assertEquals(0.0, correctRoot.getReal(), 1.0e-20);
687                 if (maxOrder > 0) {
688                     assertTrue(Double.isInfinite(correctRoot.getPartialDerivative(1).getReal()));
689                     assertTrue(correctRoot.getPartialDerivative(1).getReal() > 0);
690                     for (int order = 2; order <= maxOrder; ++order) {
691                         assertTrue(Double.isInfinite(correctRoot.getPartialDerivative(order).getReal()));
692                         if (signedInfinities) {
693                             if ((order % 2) == 0) {
694                                 assertTrue(correctRoot.getPartialDerivative(order).getReal() < 0);
695                             } else {
696                                 assertTrue(correctRoot.getPartialDerivative(order).getReal() > 0);
697                             }
698                         }
699                     }
700                 }
701 
702             }
703 
704         }
705 
706     }
707 
708     @Test
709     public void testSqrtPow2() {
710         double[] epsilon = new double[] { 1.0e-16, 3.0e-16, 2.0e-15, 6.0e-14, 6.0e-12 };
711         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
712             final FDSFactory<T> factory = buildFactory(1, maxOrder);
713             for (double x = 0.1; x < 1.2; x += 0.001) {
714                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
715                 FieldDerivativeStructure<T> rebuiltX = dsX.multiply(dsX).sqrt();
716                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
717                 for (int n = 0; n <= maxOrder; ++n) {
718                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
719                 }
720             }
721         }
722     }
723 
724     @Test
725     public void testCbrtDefinition() {
726         double[] epsilon = new double[] { 4.0e-16, 9.0e-16, 6.0e-15, 2.0e-13, 4.0e-12 };
727         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
728             final FDSFactory<T> factory = buildFactory(1, maxOrder);
729             for (double x = 0.1; x < 1.2; x += 0.001) {
730                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
731                 FieldDerivativeStructure<T> cbrt1 = dsX.pow(1.0 / 3.0);
732                 FieldDerivativeStructure<T> cbrt2 = dsX.cbrt();
733                 FieldDerivativeStructure<T> zero = cbrt1.subtract(cbrt2);
734                 for (int n = 0; n <= maxOrder; ++n) {
735                     assertEquals(0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
736                 }
737             }
738         }
739     }
740 
741     @Test
742     public void testCbrtPow3() {
743         double[] epsilon = new double[] { 1.0e-16, 5.0e-16, 8.0e-15, 4.0e-13, 3.0e-11 };
744         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
745             final FDSFactory<T> factory = buildFactory(1, maxOrder);
746             for (double x = 0.1; x < 1.2; x += 0.001) {
747                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
748                 FieldDerivativeStructure<T> rebuiltX = dsX.multiply(dsX.multiply(dsX)).cbrt();
749                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
750                 for (int n = 0; n <= maxOrder; ++n) {
751                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
752                 }
753             }
754         }
755     }
756 
757     @Test
758     public void testPowReciprocalPow() {
759         double[] epsilon = new double[] { 2.0e-15, 2.0e-14, 3.0e-13, 8.0e-12, 3.0e-10 };
760         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
761             final FDSFactory<T> factory = buildFactory(2, maxOrder);
762             for (double x = 0.1; x < 1.2; x += 0.01) {
763                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
764                 for (double y = 0.1; y < 1.2; y += 0.01) {
765                     FieldDerivativeStructure<T> dsY = factory.variable(1, y);
766                     FieldDerivativeStructure<T> rebuiltX = dsX.pow(dsY).pow(dsY.reciprocal());
767                     FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
768                     for (int n = 0; n <= maxOrder; ++n) {
769                         for (int m = 0; m <= maxOrder; ++m) {
770                             if (n + m <= maxOrder) {
771                                 assertEquals(0.0, zero.getPartialDerivative(n, m).getReal(), epsilon[n + m]);
772                             }
773                         }
774                     }
775                 }
776             }
777         }
778     }
779 
780     @Test
781     public void testHypotDefinition() {
782         double epsilon = 1.0e-20;
783         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
784             final FDSFactory<T> factory = buildFactory(2, maxOrder);
785             for (double x = -1.7; x < 2; x += 0.2) {
786                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
787                 for (double y = -1.7; y < 2; y += 0.2) {
788                     FieldDerivativeStructure<T> dsY = factory.variable(1, y);
789                     FieldDerivativeStructure<T> hypot = FieldDerivativeStructure.hypot(dsY, dsX);
790                     FieldDerivativeStructure<T> ref = dsX.multiply(dsX).add(dsY.multiply(dsY)).sqrt();
791                     FieldDerivativeStructure<T> zero = hypot.subtract(ref);
792                     for (int n = 0; n <= maxOrder; ++n) {
793                         for (int m = 0; m <= maxOrder; ++m) {
794                             if (n + m <= maxOrder) {
795                                 assertEquals(0, zero.getPartialDerivative(n, m).getReal(), epsilon);
796                             }
797                         }
798                     }
799                 }
800             }
801         }
802     }
803 
804     @Test
805     public abstract void testHypotNoOverflow();
806 
807     protected void doTestHypotNoOverflow(int tenPower) {
808 
809         final FDSFactory<T> factory = buildFactory(2, 5);
810         FieldDerivativeStructure<T> dsX = factory.variable(0, +3.0);
811         FieldDerivativeStructure<T> dsY = factory.variable(1, -4.0);
812         T scaling = factory.getValueField().getOne();
813         for (int i = 0; i < tenPower; ++i) {
814             scaling = scaling.multiply(10);
815         }
816         dsX = dsX.multiply(scaling);
817         dsY = dsY.multiply(scaling);
818         FieldDerivativeStructure<T> hypot = FieldDerivativeStructure.hypot(dsX, dsY);
819         FieldDerivativeStructure<T> scaledDownHypot = hypot;
820         scaledDownHypot = scaledDownHypot.divide(scaling);
821         assertEquals(5.0, scaledDownHypot.getReal(), 5.0e-15);
822         assertEquals(dsX.divide(hypot).getReal(), scaledDownHypot.getPartialDerivative(1, 0).getReal(), 1.0e-10);
823         assertEquals(dsY.divide(hypot).getReal(), scaledDownHypot.getPartialDerivative(0, 1).getReal(), 1.0e-10);
824 
825         FieldDerivativeStructure<T> sqrt  = dsX.multiply(dsX).add(dsY.multiply(dsY)).sqrt();
826         assertTrue(sqrt.getValue().isInfinite() || sqrt.getValue().isNaN());
827 
828     }
829 
830     @Test
831     public void testHypotNeglectible() {
832 
833         final FDSFactory<T> factory = buildFactory(2, 5);
834         FieldDerivativeStructure<T> dsSmall = factory.variable(0, +3.0e-10);
835         FieldDerivativeStructure<T> dsLarge = factory.variable(1, -4.0e25);
836 
837         assertEquals(dsLarge.norm(),
838                             FieldDerivativeStructure.hypot(dsSmall, dsLarge).getReal(),
839                             1.0e-10);
840         assertEquals(0,
841                             FieldDerivativeStructure.hypot(dsSmall, dsLarge).getPartialDerivative(1, 0).getReal(),
842                             1.0e-10);
843         assertEquals(-1,
844                             FieldDerivativeStructure.hypot(dsSmall, dsLarge).getPartialDerivative(0, 1).getReal(),
845                             1.0e-10);
846 
847         assertEquals(dsLarge.norm(),
848                             FieldDerivativeStructure.hypot(dsLarge, dsSmall).getReal(),
849                             1.0e-10);
850         assertEquals(0,
851                             FieldDerivativeStructure.hypot(dsLarge, dsSmall).getPartialDerivative(1, 0).getReal(),
852                             1.0e-10);
853         assertEquals(-1,
854                             FieldDerivativeStructure.hypot(dsLarge, dsSmall).getPartialDerivative(0, 1).getReal(),
855                             1.0e-10);
856 
857     }
858 
859     @Test
860     public void testHypotSpecial() {
861         final FDSFactory<T> factory = buildFactory(2, 5);
862         assertTrue(Double.isNaN(FieldDerivativeStructure.hypot(factory.variable(0, Double.NaN),
863                                                                  factory.variable(0, +3.0e250)).getReal()));
864         assertTrue(Double.isNaN(FieldDerivativeStructure.hypot(factory.variable(0, +3.0e250),
865                                                                  factory.variable(0, Double.NaN)).getReal()));
866         assertTrue(Double.isInfinite(FieldDerivativeStructure.hypot(factory.variable(0, Double.POSITIVE_INFINITY),
867                                                                       factory.variable(0, +3.0e250)).getReal()));
868         assertTrue(Double.isInfinite(FieldDerivativeStructure.hypot(factory.variable(0, +3.0e250),
869                                                                       factory.variable(0, Double.POSITIVE_INFINITY)).getReal()));
870     }
871 
872     @Test
873     public void testFieldRemainder() {
874         double epsilon = 1.0e-15;
875         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
876             final FDSFactory<T> factory = buildFactory(2, maxOrder);
877             for (double x = -1.7; x < 2; x += 0.2) {
878                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
879                 for (double y = -1.7; y < 2; y += 0.2) {
880                     FieldDerivativeStructure<T> remainder = dsX.remainder(buildScalar(y));
881                     FieldDerivativeStructure<T> ref = dsX.subtract(x - FastMath.IEEEremainder(x, y));
882                     FieldDerivativeStructure<T> zero = remainder.subtract(ref);
883                     for (int n = 0; n <= maxOrder; ++n) {
884                         for (int m = 0; m <= maxOrder; ++m) {
885                             if (n + m <= maxOrder) {
886                                 assertEquals(0, zero.getPartialDerivative(n, m).getReal(), epsilon);
887                             }
888                         }
889                     }
890                 }
891             }
892         }
893     }
894 
895     @Test
896     public void testPrimitiveRemainder() {
897         double epsilon = 1.0e-15;
898         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
899             final FDSFactory<T> factory = buildFactory(2, maxOrder);
900             for (double x = -1.7; x < 2; x += 0.2) {
901                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
902                 for (double y = -1.7; y < 2; y += 0.2) {
903                     FieldDerivativeStructure<T> remainder = dsX.remainder(y);
904                     FieldDerivativeStructure<T> ref = dsX.subtract(x - FastMath.IEEEremainder(x, y));
905                     FieldDerivativeStructure<T> zero = remainder.subtract(ref);
906                     for (int n = 0; n <= maxOrder; ++n) {
907                         for (int m = 0; m <= maxOrder; ++m) {
908                             if (n + m <= maxOrder) {
909                                 assertEquals(0, zero.getPartialDerivative(n, m).getReal(), epsilon);
910                             }
911                         }
912                     }
913                 }
914             }
915         }
916     }
917 
918     @Test
919     public void testRemainder() {
920         double epsilon = 2.0e-15;
921         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
922             final FDSFactory<T> factory = buildFactory(2, maxOrder);
923             for (double x = -1.7; x < 2; x += 0.2) {
924                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
925                 for (double y = -1.7; y < 2; y += 0.2) {
926                     FieldDerivativeStructure<T> dsY = factory.variable(1, y);
927                     FieldDerivativeStructure<T> remainder = dsX.remainder(dsY);
928                     FieldDerivativeStructure<T> ref = dsX.subtract(dsY.multiply((x - FastMath.IEEEremainder(x, y)) / y));
929                     FieldDerivativeStructure<T> zero = remainder.subtract(ref);
930                     for (int n = 0; n <= maxOrder; ++n) {
931                         for (int m = 0; m <= maxOrder; ++m) {
932                             if (n + m <= maxOrder) {
933                                 assertEquals(0, zero.getPartialDerivative(n, m).getReal(), epsilon);
934                             }
935                         }
936                     }
937                 }
938             }
939         }
940     }
941 
942     @Override
943     @Test
944     public void testExp() {
945         double[] epsilon = new double[] { 1.0e-16, 1.0e-16, 1.0e-16, 1.0e-16, 1.0e-16 };
946         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
947             final FDSFactory<T> factory = buildFactory(1, maxOrder);
948             for (double x = 0.1; x < 1.2; x += 0.001) {
949                 double refExp = FastMath.exp(x);
950                 FieldDerivativeStructure<T> exp = factory.variable(0, x).exp();
951                 for (int n = 0; n <= maxOrder; ++n) {
952                     assertEquals(refExp, exp.getPartialDerivative(n).getReal(), epsilon[n]);
953                 }
954             }
955         }
956     }
957 
958     @Test
959     public void testExpm1Definition() {
960         double epsilon = 3.0e-16;
961         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
962             final FDSFactory<T> factory = buildFactory(1, maxOrder);
963             for (double x = 0.1; x < 1.2; x += 0.001) {
964                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
965                 FieldDerivativeStructure<T> expm11 = dsX.expm1();
966                 FieldDerivativeStructure<T> expm12 = dsX.exp().subtract(dsX.getField().getOne());
967                 FieldDerivativeStructure<T> zero = expm11.subtract(expm12);
968                 for (int n = 0; n <= maxOrder; ++n) {
969                     assertEquals(0, zero.getPartialDerivative(n).getReal(), epsilon);
970                 }
971             }
972         }
973     }
974 
975     @Override
976     @Test
977     public void testLog() {
978         double[] epsilon = new double[] { 1.0e-16, 1.0e-16, 3.0e-14, 7.0e-13, 3.0e-11 };
979         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
980             final FDSFactory<T> factory = buildFactory(1, maxOrder);
981             for (double x = 0.1; x < 1.2; x += 0.001) {
982                 FieldDerivativeStructure<T> log = factory.variable(0, x).log();
983                 assertEquals(FastMath.log(x), log.getReal(), epsilon[0]);
984                 for (int n = 1; n <= maxOrder; ++n) {
985                     double refDer = -CombinatoricsUtils.factorial(n - 1) / FastMath.pow(-x, n);
986                     assertEquals(refDer, log.getPartialDerivative(n).getReal(), epsilon[n]);
987                 }
988             }
989         }
990     }
991 
992     @Test
993     public void testLog1pDefinition() {
994         double epsilon = 3.0e-16;
995         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
996             final FDSFactory<T> factory = buildFactory(1, maxOrder);
997             for (double x = 0.1; x < 1.2; x += 0.001) {
998                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
999                 FieldDerivativeStructure<T> log1p1 = dsX.log1p();
1000                 FieldDerivativeStructure<T> log1p2 = dsX.add(dsX.getField().getOne()).log();
1001                 FieldDerivativeStructure<T> zero = log1p1.subtract(log1p2);
1002                 for (int n = 0; n <= maxOrder; ++n) {
1003                     assertEquals(0, zero.getPartialDerivative(n).getReal(), epsilon);
1004                 }
1005             }
1006         }
1007     }
1008 
1009     @Test
1010     public void testLog10Definition() {
1011         double[] epsilon = new double[] { 3.0e-16, 9.0e-16, 8.0e-15, 3.0e-13, 8.0e-12 };
1012         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1013             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1014             for (double x = 0.1; x < 1.2; x += 0.001) {
1015                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1016                 FieldDerivativeStructure<T> log101 = dsX.log10();
1017                 FieldDerivativeStructure<T> log102 = dsX.log().divide(FastMath.log(10.0));
1018                 FieldDerivativeStructure<T> zero = log101.subtract(log102);
1019                 for (int n = 0; n <= maxOrder; ++n) {
1020                     assertEquals(0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1021                 }
1022             }
1023         }
1024     }
1025 
1026     @Test
1027     public void testLogExp() {
1028         double[] epsilon = new double[] { 2.0e-16, 2.0e-16, 3.0e-16, 2.0e-15, 6.0e-15 };
1029         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1030             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1031             for (double x = 0.1; x < 1.2; x += 0.001) {
1032                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1033                 FieldDerivativeStructure<T> rebuiltX = dsX.exp().log();
1034                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
1035                 for (int n = 0; n <= maxOrder; ++n) {
1036                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1037                 }
1038             }
1039         }
1040     }
1041 
1042     @Test
1043     public void testLog1pExpm1() {
1044         double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 5.0e-16, 9.0e-16, 6.0e-15 };
1045         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1046             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1047             for (double x = 0.1; x < 1.2; x += 0.001) {
1048                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1049                 FieldDerivativeStructure<T> rebuiltX = dsX.expm1().log1p();
1050                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
1051                 for (int n = 0; n <= maxOrder; ++n) {
1052                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1053                 }
1054             }
1055         }
1056     }
1057 
1058     @Test
1059     public void testLog10Power() {
1060         double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 9.0e-16, 6.0e-15, 7.0e-14 };
1061         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1062             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1063             for (double x = 0.1; x < 1.2; x += 0.001) {
1064                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1065                 FieldDerivativeStructure<T> rebuiltX = factory.constant(10.0).pow(dsX).log10();
1066                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
1067                 for (int n = 0; n <= maxOrder; ++n) {
1068                     assertEquals(0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1069                 }
1070             }
1071         }
1072     }
1073 
1074     @Test
1075     public void testSinCosSeparated() {
1076         double epsilon = 5.0e-16;
1077         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1078             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1079             for (double x = 0.1; x < 1.2; x += 0.001) {
1080                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1081                 FieldDerivativeStructure<T> sin = dsX.sin();
1082                 FieldDerivativeStructure<T> cos = dsX.cos();
1083                 double s = FastMath.sin(x);
1084                 double c = FastMath.cos(x);
1085                 for (int n = 0; n <= maxOrder; ++n) {
1086                     switch (n % 4) {
1087                     case 0 :
1088                         assertEquals( s, sin.getPartialDerivative(n).getReal(), epsilon);
1089                         assertEquals( c, cos.getPartialDerivative(n).getReal(), epsilon);
1090                         break;
1091                     case 1 :
1092                         assertEquals( c, sin.getPartialDerivative(n).getReal(), epsilon);
1093                         assertEquals(-s, cos.getPartialDerivative(n).getReal(), epsilon);
1094                         break;
1095                     case 2 :
1096                         assertEquals(-s, sin.getPartialDerivative(n).getReal(), epsilon);
1097                         assertEquals(-c, cos.getPartialDerivative(n).getReal(), epsilon);
1098                         break;
1099                     default :
1100                         assertEquals(-c, sin.getPartialDerivative(n).getReal(), epsilon);
1101                         assertEquals( s, cos.getPartialDerivative(n).getReal(), epsilon);
1102                         break;
1103                     }
1104                 }
1105             }
1106         }
1107     }
1108 
1109     @Test
1110     public void testSinCosCombined() {
1111         double epsilon = 5.0e-16;
1112         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1113             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1114             for (double x = 0.1; x < 1.2; x += 0.001) {
1115                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1116                 FieldSinCos<FieldDerivativeStructure<T>> sinCos = dsX.sinCos();
1117                 double s = FastMath.sin(x);
1118                 double c = FastMath.cos(x);
1119                 for (int n = 0; n <= maxOrder; ++n) {
1120                     switch (n % 4) {
1121                     case 0 :
1122                         assertEquals( s, sinCos.sin().getPartialDerivative(n).getReal(), epsilon);
1123                         assertEquals( c, sinCos.cos().getPartialDerivative(n).getReal(), epsilon);
1124                         break;
1125                     case 1 :
1126                         assertEquals( c, sinCos.sin().getPartialDerivative(n).getReal(), epsilon);
1127                         assertEquals(-s, sinCos.cos().getPartialDerivative(n).getReal(), epsilon);
1128                         break;
1129                     case 2 :
1130                         assertEquals(-s, sinCos.sin().getPartialDerivative(n).getReal(), epsilon);
1131                         assertEquals(-c, sinCos.cos().getPartialDerivative(n).getReal(), epsilon);
1132                         break;
1133                     default :
1134                         assertEquals(-c, sinCos.sin().getPartialDerivative(n).getReal(), epsilon);
1135                         assertEquals( s, sinCos.cos().getPartialDerivative(n).getReal(), epsilon);
1136                         break;
1137                     }
1138                 }
1139             }
1140         }
1141     }
1142 
1143     @Test
1144     public void testSinAsin() {
1145         double[] epsilon = new double[] { 3.0e-16, 5.0e-16, 3.0e-15, 2.0e-14, 4.0e-13 };
1146         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1147             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1148             for (double x = 0.1; x < 1.2; x += 0.001) {
1149                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1150                 FieldDerivativeStructure<T> rebuiltX = dsX.sin().asin();
1151                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
1152                 for (int n = 0; n <= maxOrder; ++n) {
1153                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1154                 }
1155             }
1156         }
1157     }
1158 
1159     @Test
1160     public void testCosAcos() {
1161         double[] epsilon = new double[] { 7.0e-16, 6.0e-15, 2.0e-13, 4.0e-12, 2.0e-10 };
1162         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1163             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1164             for (double x = 0.1; x < 1.2; x += 0.001) {
1165                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1166                 FieldDerivativeStructure<T> rebuiltX = dsX.cos().acos();
1167                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
1168                 for (int n = 0; n <= maxOrder; ++n) {
1169                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1170                 }
1171             }
1172         }
1173     }
1174 
1175     @Test
1176     public void testTanAtan() {
1177         double[] epsilon = new double[] { 3.0e-16, 2.0e-16, 2.0e-15, 4.0e-14, 2.0e-12 };
1178         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1179             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1180             for (double x = 0.1; x < 1.2; x += 0.001) {
1181                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1182                 FieldDerivativeStructure<T> rebuiltX = dsX.tan().atan();
1183                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
1184                 for (int n = 0; n <= maxOrder; ++n) {
1185                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1186                 }
1187             }
1188         }
1189     }
1190 
1191     @Test
1192     public void testTangentDefinition() {
1193         double[] epsilon = new double[] { 9.0e-16, 4.0e-15, 4.0e-14, 5.0e-13, 2.0e-11 };
1194         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1195             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1196             for (double x = 0.1; x < 1.2; x += 0.001) {
1197                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1198                 FieldDerivativeStructure<T> tan1 = dsX.sin().divide(dsX.cos());
1199                 FieldDerivativeStructure<T> tan2 = dsX.tan();
1200                 FieldDerivativeStructure<T> zero = tan1.subtract(tan2);
1201                 for (int n = 0; n <= maxOrder; ++n) {
1202                     assertEquals(0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1203                 }
1204             }
1205         }
1206     }
1207 
1208     @Override
1209     @Test
1210     public void testAtan2() {
1211         double[] epsilon = new double[] { 5.0e-16, 3.0e-15, 2.9e-14, 1.0e-12, 8.0e-11 };
1212         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1213             final FDSFactory<T> factory = buildFactory(2, maxOrder);
1214             for (double x = -1.7; x < 2; x += 0.2) {
1215                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1216                 for (double y = -1.7; y < 2; y += 0.2) {
1217                     FieldDerivativeStructure<T> dsY = factory.variable(1, y);
1218                     FieldDerivativeStructure<T> atan2 = FieldDerivativeStructure.atan2(dsY, dsX);
1219                     FieldDerivativeStructure<T> ref = dsY.divide(dsX).atan();
1220                     if (x < 0) {
1221                         ref = (y < 0) ? ref.subtract(FastMath.PI) : ref.add(FastMath.PI);
1222                     }
1223                     FieldDerivativeStructure<T> zero = atan2.subtract(ref);
1224                     for (int n = 0; n <= maxOrder; ++n) {
1225                         for (int m = 0; m <= maxOrder; ++m) {
1226                             if (n + m <= maxOrder) {
1227                                 assertEquals(0, zero.getPartialDerivative(n, m).getReal(), epsilon[n + m]);
1228                             }
1229                         }
1230                     }
1231                 }
1232             }
1233         }
1234     }
1235 
1236     @Test
1237     public void testAtan2SpecialCasesDerivatives() {
1238 
1239         final FDSFactory<T> factory = buildFactory(2, 2);
1240         FieldDerivativeStructure<T> pp =
1241                 FieldDerivativeStructure.atan2(factory.variable(1, buildScalar(+0.0)), factory.variable(1, buildScalar(+0.0)));
1242         assertEquals(0, pp.getReal(), 1.0e-15);
1243         assertEquals(+1, FastMath.copySign(1, pp.getReal()), 1.0e-15);
1244 
1245         FieldDerivativeStructure<T> pn =
1246                 FieldDerivativeStructure.atan2(factory.variable(1, buildScalar(+0.0)), factory.variable(1, buildScalar(-0.0)));
1247         assertEquals(FastMath.PI, pn.getReal(), 1.0e-15);
1248 
1249         FieldDerivativeStructure<T> np =
1250                 FieldDerivativeStructure.atan2(factory.variable(1, buildScalar(-0.0)), factory.variable(1, buildScalar(+0.0)));
1251         assertEquals(0, np.getReal(), 1.0e-15);
1252         assertEquals(-1, FastMath.copySign(1, np.getReal()), 1.0e-15);
1253 
1254         FieldDerivativeStructure<T> nn =
1255                 FieldDerivativeStructure.atan2(factory.variable(1, buildScalar(-0.0)), factory.variable(1, buildScalar(-0.0)));
1256         assertEquals(-FastMath.PI, nn.getReal(), 1.0e-15);
1257 
1258     }
1259 
1260     @Test
1261     public void testSinhCoshCombined() {
1262         double epsilon = 5.0e-16;
1263         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1264             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1265             for (double x = 0.1; x < 1.2; x += 0.001) {
1266                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1267                 FieldSinhCosh<FieldDerivativeStructure<T>> sinhCosh = dsX.sinhCosh();
1268                 double sh = FastMath.sinh(x);
1269                 double ch = FastMath.cosh(x);
1270                 for (int n = 0; n <= maxOrder; ++n) {
1271                     if (n % 2 == 0) {
1272                         assertEquals(sh, sinhCosh.sinh().getPartialDerivative(n).getReal(), epsilon);
1273                         assertEquals(ch, sinhCosh.cosh().getPartialDerivative(n).getReal(), epsilon);
1274                     } else {
1275                         assertEquals(ch, sinhCosh.sinh().getPartialDerivative(n).getReal(), epsilon);
1276                         assertEquals(sh, sinhCosh.cosh().getPartialDerivative(n).getReal(), epsilon);
1277                     }
1278                 }
1279             }
1280         }
1281     }
1282 
1283     @Test
1284     public void testSinhDefinition() {
1285         double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 5.0e-16, 2.0e-15, 6.0e-15 };
1286         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1287             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1288             for (double x = 0.1; x < 1.2; x += 0.001) {
1289                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1290                 FieldDerivativeStructure<T> sinh1 = dsX.exp().subtract(dsX.exp().reciprocal()).multiply(0.5);
1291                 FieldDerivativeStructure<T> sinh2 = dsX.sinh();
1292                 FieldDerivativeStructure<T> zero = sinh1.subtract(sinh2);
1293                 for (int n = 0; n <= maxOrder; ++n) {
1294                     assertEquals(0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1295                 }
1296             }
1297         }
1298     }
1299 
1300     @Test
1301     public void testCoshDefinition() {
1302         double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 5.0e-16, 2.0e-15, 6.0e-15 };
1303         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1304             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1305             for (double x = 0.1; x < 1.2; x += 0.001) {
1306                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1307                 FieldDerivativeStructure<T> cosh1 = dsX.exp().add(dsX.exp().reciprocal()).multiply(0.5);
1308                 FieldDerivativeStructure<T> cosh2 = dsX.cosh();
1309                 FieldDerivativeStructure<T> zero = cosh1.subtract(cosh2);
1310                 for (int n = 0; n <= maxOrder; ++n) {
1311                     assertEquals(0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1312                 }
1313             }
1314         }
1315     }
1316 
1317     @Test
1318     public void testTanhDefinition() {
1319         double[] epsilon = new double[] { 3.0e-16, 5.0e-16, 7.0e-16, 3.0e-15, 2.0e-14 };
1320         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1321             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1322             for (double x = 0.1; x < 1.2; x += 0.001) {
1323                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1324                 FieldDerivativeStructure<T> tanh1 = dsX.exp().subtract(dsX.exp().reciprocal()).divide(dsX.exp().add(dsX.exp().reciprocal()));
1325                 FieldDerivativeStructure<T> tanh2 = dsX.tanh();
1326                 FieldDerivativeStructure<T> zero = tanh1.subtract(tanh2);
1327                 for (int n = 0; n <= maxOrder; ++n) {
1328                     assertEquals(0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1329                 }
1330             }
1331         }
1332     }
1333 
1334     @Test
1335     public void testSinhAsinh() {
1336         double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 4.0e-16, 7.0e-16, 3.0e-15, 8.0e-15 };
1337         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1338             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1339             for (double x = 0.1; x < 1.2; x += 0.001) {
1340                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1341                 FieldDerivativeStructure<T> rebuiltX = dsX.sinh().asinh();
1342                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
1343                 for (int n = 0; n <= maxOrder; ++n) {
1344                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1345                 }
1346             }
1347         }
1348     }
1349 
1350     @Test
1351     public void testCoshAcosh() {
1352         double[] epsilon = new double[] { 2.0e-15, 1.0e-14, 2.0e-13, 6.0e-12, 3.0e-10, 2.0e-8 };
1353         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1354             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1355             for (double x = 0.1; x < 1.2; x += 0.001) {
1356                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1357                 FieldDerivativeStructure<T> rebuiltX = dsX.cosh().acosh();
1358                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
1359                 for (int n = 0; n <= maxOrder; ++n) {
1360                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1361                 }
1362             }
1363         }
1364     }
1365 
1366     @Test
1367     public void testTanhAtanh() {
1368         double[] epsilon = new double[] { 5.0e-16, 2.0e-16, 7.0e-16, 4.0e-15, 3.0e-14, 4.0e-13 };
1369         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1370             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1371             for (double x = 0.1; x < 1.2; x += 0.001) {
1372                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1373                 FieldDerivativeStructure<T> rebuiltX = dsX.tanh().atanh();
1374                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
1375                 for (int n = 0; n <= maxOrder; ++n) {
1376                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon[n]);
1377                 }
1378             }
1379         }
1380     }
1381 
1382     @Test
1383     public void testCompositionOneVariableY() {
1384         double epsilon = 1.0e-13;
1385         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1386             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1387             for (double x = 0.1; x < 1.2; x += 0.1) {
1388                 FieldDerivativeStructure<T> dsX = factory.constant(x);
1389                 for (double y = 0.1; y < 1.2; y += 0.1) {
1390                     FieldDerivativeStructure<T> dsY = factory.variable(0, y);
1391                     FieldDerivativeStructure<T> f = dsX.divide(dsY).sqrt();
1392                     double f0 = FastMath.sqrt(x / y);
1393                     assertEquals(f0, f.getReal(), FastMath.abs(epsilon * f0));
1394                     if (f.getOrder() > 0) {
1395                         double f1 = -x / (2 * y * y * f0);
1396                         assertEquals(f1, f.getPartialDerivative(1).getReal(), FastMath.abs(epsilon * f1));
1397                         if (f.getOrder() > 1) {
1398                             double f2 = (f0 - x / (4 * y * f0)) / (y * y);
1399                             assertEquals(f2, f.getPartialDerivative(2).getReal(), FastMath.abs(epsilon * f2));
1400                             if (f.getOrder() > 2) {
1401                                 double f3 = (x / (8 * y * f0) - 2 * f0) / (y * y * y);
1402                                 assertEquals(f3, f.getPartialDerivative(3).getReal(), FastMath.abs(epsilon * f3));
1403                             }
1404                         }
1405                     }
1406                 }
1407             }
1408         }
1409     }
1410 
1411     @Test
1412     public void testTaylorPrimitivePolynomial() {
1413         final FDSFactory<T> factory = buildFactory(3, 4);
1414         for (double x = 0; x < 1.2; x += 0.1) {
1415             FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1416             for (double y = 0; y < 1.2; y += 0.2) {
1417                 FieldDerivativeStructure<T> dsY = factory.variable(1, y);
1418                 for (double z = 0; z < 1.2; z += 0.2) {
1419                     FieldDerivativeStructure<T> dsZ = factory.variable(2, z);
1420                     FieldDerivativeStructure<T> f = dsX.multiply(dsY).add(dsZ).multiply(dsX).multiply(dsY);
1421                     for (double dx = -0.2; dx < 0.2; dx += 0.2) {
1422                         for (double dy = -0.2; dy < 0.2; dy += 0.1) {
1423                             for (double dz = -0.2; dz < 0.2; dz += 0.1) {
1424                                 double ref = (x + dx) * (y + dy) * ((x + dx) * (y + dy) + (z + dz));
1425                                 assertEquals(ref, f.taylor(dx, dy, dz).getReal(), 2.0e-15);
1426                             }
1427                         }
1428                     }
1429                 }
1430             }
1431         }
1432     }
1433 
1434     @Test
1435     public void testTaylorFieldPolynomial() {
1436         final FDSFactory<T> factory = buildFactory(3, 4);
1437         for (double x = 0; x < 1.2; x += 0.1) {
1438             FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1439             for (double y = 0; y < 1.2; y += 0.2) {
1440                 FieldDerivativeStructure<T> dsY = factory.variable(1, y);
1441                 for (double z = 0; z < 1.2; z += 0.2) {
1442                     FieldDerivativeStructure<T> dsZ = factory.variable(2, z);
1443                     FieldDerivativeStructure<T> f = dsX.multiply(dsY).add(dsZ).multiply(dsX).multiply(dsY);
1444                     for (double dx = -0.2; dx < 0.2; dx += 0.2) {
1445                         T dxF = buildScalar(dx);
1446                         for (double dy = -0.2; dy < 0.2; dy += 0.1) {
1447                             T dyF = buildScalar(dy);
1448                             for (double dz = -0.2; dz < 0.2; dz += 0.1) {
1449                                 T dzF = buildScalar(dz);
1450                                 double ref = (x + dx) * (y + dy) * ((x + dx) * (y + dy) + (z + dz));
1451                                 assertEquals(ref, f.taylor(dxF, dyF, dzF).getReal(), 2.0e-15);
1452                             }
1453                         }
1454                     }
1455                 }
1456             }
1457         }
1458     }
1459 
1460     @Test
1461     public void testTaylorAtan2() {
1462         double[] expected = new double[] { 0.214, 0.0241, 0.00422, 6.48e-4, 8.04e-5 };
1463         double x0 =  0.1;
1464         double y0 = -0.3;
1465         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1466             final FDSFactory<T> factory = buildFactory(2, maxOrder);
1467             FieldDerivativeStructure<T> dsX   = factory.variable(0, x0);
1468             FieldDerivativeStructure<T> dsY   = factory.variable(1, y0);
1469             FieldDerivativeStructure<T> atan2 = FieldDerivativeStructure.atan2(dsY, dsX);
1470             double maxError = 0;
1471             for (double dx = -0.05; dx < 0.05; dx += 0.001) {
1472                 for (double dy = -0.05; dy < 0.05; dy += 0.001) {
1473                     double ref = FastMath.atan2(y0 + dy, x0 + dx);
1474                     maxError = FastMath.max(maxError, FastMath.abs(ref - atan2.taylor(dx, dy).getReal()));
1475                 }
1476             }
1477             assertEquals(0.0, expected[maxOrder] - maxError, 0.01 * expected[maxOrder]);
1478         }
1479     }
1480 
1481     @Test
1482     public void testNorm() {
1483 
1484         final FDSFactory<T> factory = buildFactory(1, 1);
1485         FieldDerivativeStructure<T> minusOne = factory.variable(0, -1.0);
1486         assertEquals(+1.0, minusOne.abs().getPartialDerivative(0).getReal(), 1.0e-15);
1487         assertEquals(-1.0, minusOne.abs().getPartialDerivative(1).getReal(), 1.0e-15);
1488 
1489         FieldDerivativeStructure<T> plusOne = factory.variable(0, +1.0);
1490         assertEquals(+1.0, plusOne.abs().getPartialDerivative(0).getReal(), 1.0e-15);
1491         assertEquals(+1.0, plusOne.abs().getPartialDerivative(1).getReal(), 1.0e-15);
1492 
1493         FieldDerivativeStructure<T> minusZero = factory.variable(0, buildScalar(-0.0));
1494         assertEquals(+0.0, minusZero.abs().getPartialDerivative(0).getReal(), 1.0e-15);
1495         assertEquals(-1.0, minusZero.abs().getPartialDerivative(1).getReal(), 1.0e-15);
1496 
1497         FieldDerivativeStructure<T> plusZero = factory.variable(0, buildScalar(+0.0));
1498         assertEquals(+0.0, plusZero.abs().getPartialDerivative(0).getReal(), 1.0e-15);
1499         assertEquals(+1.0, plusZero.abs().getPartialDerivative(1).getReal(), 1.0e-15);
1500 
1501     }
1502 
1503     @Override
1504     @Test
1505     public void testSign() {
1506 
1507         final FDSFactory<T> factory = buildFactory(1, 1);
1508         FieldDerivativeStructure<T> minusOne = factory.variable(0, -1.0);
1509         assertEquals(-1.0, minusOne.sign().getPartialDerivative(0).getReal(), 1.0e-15);
1510         assertEquals( 0.0, minusOne.sign().getPartialDerivative(1).getReal(), 1.0e-15);
1511 
1512         FieldDerivativeStructure<T> plusOne = factory.variable(0, +1.0);
1513         assertEquals(+1.0, plusOne.sign().getPartialDerivative(0).getReal(), 1.0e-15);
1514         assertEquals( 0.0, plusOne.sign().getPartialDerivative(1).getReal(), 1.0e-15);
1515 
1516         FieldDerivativeStructure<T> minusZero = factory.variable(0, buildScalar(-0.0));
1517         assertEquals(-0.0, minusZero.sign().getPartialDerivative(0).getReal(), 1.0e-15);
1518         assertTrue(Double.doubleToLongBits(minusZero.sign().getReal()) < 0);
1519         assertEquals( 0.0, minusZero.sign().getPartialDerivative(1).getReal(), 1.0e-15);
1520 
1521         FieldDerivativeStructure<T> plusZero = factory.variable(0, buildScalar(+0.0));
1522         assertEquals(+0.0, plusZero.sign().getPartialDerivative(0).getReal(), 1.0e-15);
1523         assertEquals(0, Double.doubleToLongBits(plusZero.sign().getReal()));
1524         assertEquals( 0.0, plusZero.sign().getPartialDerivative(1).getReal(), 1.0e-15);
1525 
1526     }
1527 
1528     @Test
1529     public void testCeilFloorRintLong() {
1530 
1531         final FDSFactory<T> factory = buildFactory(1, 1);
1532         FieldDerivativeStructure<T> x = factory.variable(0, -1.5);
1533         assertEquals(-1.5, x.getPartialDerivative(0).getReal(), 1.0e-15);
1534         assertEquals(+1.0, x.getPartialDerivative(1).getReal(), 1.0e-15);
1535         assertEquals(-1.0, x.ceil().getPartialDerivative(0).getReal(), 1.0e-15);
1536         assertEquals(+0.0, x.ceil().getPartialDerivative(1).getReal(), 1.0e-15);
1537         assertEquals(-2.0, x.floor().getPartialDerivative(0).getReal(), 1.0e-15);
1538         assertEquals(+0.0, x.floor().getPartialDerivative(1).getReal(), 1.0e-15);
1539         assertEquals(-2.0, x.rint().getPartialDerivative(0).getReal(), 1.0e-15);
1540         assertEquals(+0.0, x.rint().getPartialDerivative(1).getReal(), 1.0e-15);
1541         assertEquals(-2.0, x.subtract(x.getField().getOne()).rint().getPartialDerivative(0).getReal(), 1.0e-15);
1542 
1543     }
1544 
1545     @Test
1546     public void testCopySign() {
1547 
1548         final FDSFactory<T> factory = buildFactory(1, 1);
1549         FieldDerivativeStructure<T> minusOne = factory.variable(0, -1.0);
1550         assertEquals(+1.0, minusOne.copySign(+1.0).getPartialDerivative(0).getReal(), 1.0e-15);
1551         assertEquals(-1.0, minusOne.copySign(+1.0).getPartialDerivative(1).getReal(), 1.0e-15);
1552         assertEquals(-1.0, minusOne.copySign(-1.0).getPartialDerivative(0).getReal(), 1.0e-15);
1553         assertEquals(+1.0, minusOne.copySign(-1.0).getPartialDerivative(1).getReal(), 1.0e-15);
1554         assertEquals(+1.0, minusOne.copySign(+0.0).getPartialDerivative(0).getReal(), 1.0e-15);
1555         assertEquals(-1.0, minusOne.copySign(+0.0).getPartialDerivative(1).getReal(), 1.0e-15);
1556         assertEquals(-1.0, minusOne.copySign(-0.0).getPartialDerivative(0).getReal(), 1.0e-15);
1557         assertEquals(+1.0, minusOne.copySign(-0.0).getPartialDerivative(1).getReal(), 1.0e-15);
1558         assertEquals(+1.0, minusOne.copySign(Double.NaN).getPartialDerivative(0).getReal(), 1.0e-15);
1559         assertEquals(-1.0, minusOne.copySign(Double.NaN).getPartialDerivative(1).getReal(), 1.0e-15);
1560         assertEquals(+1.0, minusOne.copySign(buildScalar(+1.0)).getPartialDerivative(0).getReal(), 1.0e-15);
1561         assertEquals(-1.0, minusOne.copySign(buildScalar(+1.0)).getPartialDerivative(1).getReal(), 1.0e-15);
1562         assertEquals(-1.0, minusOne.copySign(buildScalar(-1.0)).getPartialDerivative(0).getReal(), 1.0e-15);
1563         assertEquals(+1.0, minusOne.copySign(buildScalar(-1.0)).getPartialDerivative(1).getReal(), 1.0e-15);
1564         assertEquals(+1.0, minusOne.copySign(buildScalar(+0.0)).getPartialDerivative(0).getReal(), 1.0e-15);
1565         assertEquals(-1.0, minusOne.copySign(buildScalar(+0.0)).getPartialDerivative(1).getReal(), 1.0e-15);
1566         assertEquals(-1.0, minusOne.copySign(buildScalar(-0.0)).getPartialDerivative(0).getReal(), 1.0e-15);
1567         assertEquals(+1.0, minusOne.copySign(buildScalar(-0.0)).getPartialDerivative(1).getReal(), 1.0e-15);
1568         assertEquals(+1.0, minusOne.copySign(buildScalar(Double.NaN)).getPartialDerivative(0).getReal(), 1.0e-15);
1569         assertEquals(-1.0, minusOne.copySign(buildScalar(Double.NaN)).getPartialDerivative(1).getReal(), 1.0e-15);
1570 
1571         FieldDerivativeStructure<T> plusOne = factory.variable(0, +1.0);
1572         assertEquals(+1.0, plusOne.copySign(+1.0).getPartialDerivative(0).getReal(), 1.0e-15);
1573         assertEquals(+1.0, plusOne.copySign(+1.0).getPartialDerivative(1).getReal(), 1.0e-15);
1574         assertEquals(-1.0, plusOne.copySign(-1.0).getPartialDerivative(0).getReal(), 1.0e-15);
1575         assertEquals(-1.0, plusOne.copySign(-1.0).getPartialDerivative(1).getReal(), 1.0e-15);
1576         assertEquals(+1.0, plusOne.copySign(+0.0).getPartialDerivative(0).getReal(), 1.0e-15);
1577         assertEquals(+1.0, plusOne.copySign(+0.0).getPartialDerivative(1).getReal(), 1.0e-15);
1578         assertEquals(-1.0, plusOne.copySign(-0.0).getPartialDerivative(0).getReal(), 1.0e-15);
1579         assertEquals(-1.0, plusOne.copySign(-0.0).getPartialDerivative(1).getReal(), 1.0e-15);
1580         assertEquals(+1.0, plusOne.copySign(Double.NaN).getPartialDerivative(0).getReal(), 1.0e-15);
1581         assertEquals(+1.0, plusOne.copySign(Double.NaN).getPartialDerivative(1).getReal(), 1.0e-15);
1582         assertEquals(+1.0, plusOne.copySign(buildScalar(+1.0)).getPartialDerivative(0).getReal(), 1.0e-15);
1583         assertEquals(+1.0, plusOne.copySign(buildScalar(+1.0)).getPartialDerivative(1).getReal(), 1.0e-15);
1584         assertEquals(-1.0, plusOne.copySign(buildScalar(-1.0)).getPartialDerivative(0).getReal(), 1.0e-15);
1585         assertEquals(-1.0, plusOne.copySign(buildScalar(-1.0)).getPartialDerivative(1).getReal(), 1.0e-15);
1586         assertEquals(+1.0, plusOne.copySign(buildScalar(+0.0)).getPartialDerivative(0).getReal(), 1.0e-15);
1587         assertEquals(+1.0, plusOne.copySign(buildScalar(+0.0)).getPartialDerivative(1).getReal(), 1.0e-15);
1588         assertEquals(-1.0, plusOne.copySign(buildScalar(-0.0)).getPartialDerivative(0).getReal(), 1.0e-15);
1589         assertEquals(-1.0, plusOne.copySign(buildScalar(-0.0)).getPartialDerivative(1).getReal(), 1.0e-15);
1590         assertEquals(+1.0, plusOne.copySign(buildScalar(Double.NaN)).getPartialDerivative(0).getReal(), 1.0e-15);
1591         assertEquals(+1.0, plusOne.copySign(buildScalar(Double.NaN)).getPartialDerivative(1).getReal(), 1.0e-15);
1592 
1593     }
1594 
1595     @Test
1596     public void testToDegreesDefinition() {
1597         double epsilon = 3.0e-16;
1598         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1599             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1600             for (double x = 0.1; x < 1.2; x += 0.001) {
1601                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1602                 assertEquals(FastMath.toDegrees(x), dsX.toDegrees().getReal(), epsilon * FastMath.toDegrees(x));
1603                 for (int n = 1; n <= maxOrder; ++n) {
1604                     if (n == 1) {
1605                         assertEquals(180 / FastMath.PI, dsX.toDegrees().getPartialDerivative(1).getReal(), epsilon);
1606                     } else {
1607                         assertEquals(0.0, dsX.toDegrees().getPartialDerivative(n).getReal(), epsilon);
1608                     }
1609                 }
1610             }
1611         }
1612     }
1613 
1614     @Test
1615     public void testToRadiansDefinition() {
1616         double epsilon = 3.0e-16;
1617         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1618             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1619             for (double x = 0.1; x < 1.2; x += 0.001) {
1620                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1621                 assertEquals(FastMath.toRadians(x), dsX.toRadians().getReal(), epsilon);
1622                 for (int n = 1; n <= maxOrder; ++n) {
1623                     if (n == 1) {
1624                         assertEquals(FastMath.PI / 180, dsX.toRadians().getPartialDerivative(1).getReal(), epsilon);
1625                     } else {
1626                         assertEquals(0.0, dsX.toRadians().getPartialDerivative(n).getReal(), epsilon);
1627                     }
1628                 }
1629             }
1630         }
1631     }
1632 
1633     @Test
1634     public void testDegRad() {
1635         double epsilon = 3.0e-16;
1636         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1637             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1638             for (double x = 0.1; x < 1.2; x += 0.001) {
1639                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1640                 FieldDerivativeStructure<T> rebuiltX = dsX.toDegrees().toRadians();
1641                 FieldDerivativeStructure<T> zero = rebuiltX.subtract(dsX);
1642                 for (int n = 0; n <= maxOrder; ++n) {
1643                     assertEquals(0.0, zero.getPartialDerivative(n).getReal(), epsilon);
1644                 }
1645             }
1646         }
1647     }
1648 
1649     @Test
1650     public void testComposeMismatchedDimensions() {
1651         assertThrows(MathIllegalArgumentException.class, () -> {
1652             final FDSFactory<T> factory = buildFactory(1, 3);
1653             factory.variable(0, 1.2).compose(new double[3]);
1654         });
1655     }
1656 
1657     @Test
1658     public abstract void testComposeField();
1659 
1660     protected void doTestComposeField(final double[] epsilon) {
1661         double[] maxError = new double[epsilon.length];
1662         for (int maxOrder = 0; maxOrder < epsilon.length; ++maxOrder) {
1663             @SuppressWarnings("unchecked")
1664             FieldPolynomialFunction<T>[] p = (FieldPolynomialFunction<T>[]) Array.newInstance(FieldPolynomialFunction.class,
1665                                                                                               maxOrder + 1);
1666             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1667             T[] coefficients = MathArrays.buildArray(factory.getValueField(), epsilon.length);
1668             for (int i = 0; i < coefficients.length; ++i) {
1669                 coefficients[i] = factory.getValueField().getZero().newInstance(i + 1);
1670             }
1671             p[0] = new FieldPolynomialFunction<>(coefficients);
1672             for (int i = 1; i <= maxOrder; ++i) {
1673                 p[i] = p[i - 1].polynomialDerivative();
1674             }
1675             for (double x = 0.1; x < 1.2; x += 0.001) {
1676                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1677                 FieldDerivativeStructure<T> dsY1 = dsX.getField().getZero();
1678                 for (int i = p[0].degree(); i >= 0; --i) {
1679                     dsY1 = dsY1.multiply(dsX).add(p[0].getCoefficients()[i]);
1680                 }
1681                 T[] f = MathArrays.buildArray(getField(), maxOrder + 1);
1682                 for (int i = 0; i < f.length; ++i) {
1683                     f[i] = p[i].value(x);
1684                 }
1685                 FieldDerivativeStructure<T> dsY2 = dsX.compose(f);
1686                 FieldDerivativeStructure<T> zero = dsY1.subtract(dsY2);
1687                 for (int n = 0; n <= maxOrder; ++n) {
1688                     maxError[n] = FastMath.max(maxError[n], FastMath.abs(zero.getPartialDerivative(n).getReal()));
1689                 }
1690             }
1691         }
1692         for (int n = 0; n < maxError.length; ++n) {
1693             assertEquals(0.0, maxError[n], epsilon[n]);
1694         }
1695     }
1696 
1697     @Test
1698     public abstract void testComposePrimitive();
1699 
1700     protected void doTestComposePrimitive(final double[] epsilon) {
1701         PolynomialFunction poly =
1702                 new PolynomialFunction(new double[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 });
1703         double[] maxError = new double[epsilon.length];
1704         for (int maxOrder = 0; maxOrder < epsilon.length; ++maxOrder) {
1705             final FDSFactory<T> factory = buildFactory(1, maxOrder);
1706             PolynomialFunction[] p = new PolynomialFunction[maxOrder + 1];
1707             p[0] = poly;
1708             for (int i = 1; i <= maxOrder; ++i) {
1709                 p[i] = p[i - 1].polynomialDerivative();
1710             }
1711             for (double x = 0.1; x < 1.2; x += 0.001) {
1712                 FieldDerivativeStructure<T> dsX = factory.variable(0, x);
1713                 FieldDerivativeStructure<T> dsY1 = dsX.getField().getZero();
1714                 for (int i = poly.degree(); i >= 0; --i) {
1715                     dsY1 = dsY1.multiply(dsX).add(poly.getCoefficients()[i]);
1716                 }
1717                 double[] f = new double[maxOrder + 1];
1718                 for (int i = 0; i < f.length; ++i) {
1719                     f[i] = p[i].value(x);
1720                 }
1721                 FieldDerivativeStructure<T> dsY2 = dsX.compose(f);
1722                 FieldDerivativeStructure<T> zero = dsY1.subtract(dsY2);
1723                 for (int n = 0; n <= maxOrder; ++n) {
1724                     maxError[n] = FastMath.max(maxError[n], FastMath.abs(zero.getPartialDerivative(n).getReal()));
1725                 }
1726             }
1727         }
1728         for (int n = 0; n < maxError.length; ++n) {
1729             assertEquals(0.0, maxError[n], epsilon[n]);
1730         }
1731     }
1732 
1733     @Test
1734     public void testIntegration() {
1735         // check that first-order integration on two variables does not depend on sequence of operations
1736         final RandomGenerator random = new Well19937a(0x87bb96d6e11557bdl);
1737         final FDSFactory<T> factory = buildFactory(3, 7);
1738         final int size = factory.getCompiler().getSize();
1739         for (int count = 0; count < 100; ++count) {
1740             final double[] data = new double[size];
1741             for (int i = 0; i < size; i++) {
1742                 data[i] = random.nextDouble();
1743             }
1744             final FieldDerivativeStructure<T> f       = factory.build(data);
1745             final FieldDerivativeStructure<T> i2fIxIy = f.integrate(0, 1).integrate(1, 1);
1746             final FieldDerivativeStructure<T> i2fIyIx = f.integrate(1, 1).integrate(0, 1);
1747             checkEquals(i2fIxIy, i2fIyIx, 0.);
1748         }
1749     }
1750 
1751     @Test
1752     public void testIntegrationGreaterThanOrder() {
1753         // check that integration to a too high order generates zero
1754         // as integration constants are set to zero
1755         final RandomGenerator random = new Well19937a(0x4744a847b11e4c6fl);
1756         final FDSFactory<T> factory = buildFactory(3, 7);
1757         final int size = factory.getCompiler().getSize();
1758         for (int count = 0; count < 100; ++count) {
1759             final double[] data = new double[size];
1760             for (int i = 0; i < size; i++) {
1761                 data[i] = random.nextDouble();
1762             }
1763             final FieldDerivativeStructure<T> f = factory.build(data);
1764             for (int index = 0; index < factory.getCompiler().getFreeParameters(); ++index) {
1765                 final FieldDerivativeStructure<T> integ = f.integrate(index, factory.getCompiler().getOrder() + 1);
1766                 factory.constant(0).equals(integ);
1767             }
1768         }
1769     }
1770 
1771     @Test
1772     public void testIntegrationNoOp() {
1773         // check that integration of order 0 is no-op
1774         final RandomGenerator random = new Well19937a(0x75a35152f30f644bl);
1775         final FDSFactory<T> factory = buildFactory(3, 7);
1776         final int size = factory.getCompiler().getSize();
1777         for (int count = 0; count < 100; ++count) {
1778             final double[] data = new double[size];
1779             for (int i = 0; i < size; i++) {
1780                 data[i] = random.nextDouble();
1781             }
1782             final FieldDerivativeStructure<T> f = factory.build(data);
1783             for (int index = 0; index < factory.getCompiler().getFreeParameters(); ++index) {
1784                 final FieldDerivativeStructure<T> integ = f.integrate(index, 0);
1785                 checkEquals(f, integ, 0.);
1786             }
1787         }
1788     }
1789 
1790     @Test
1791     public void testDifferentiationNoOp() {
1792         // check that differentiation of order 0 is no-op
1793         final RandomGenerator random = new Well19937a(0x3b6ae4c2f1282949l);
1794         final FDSFactory<T> factory = buildFactory(3, 7);
1795         final int size = factory.getCompiler().getSize();
1796         for (int count = 0; count < 100; ++count) {
1797             final double[] data = new double[size];
1798             for (int i = 0; i < size; i++) {
1799                 data[i] = random.nextDouble();
1800             }
1801             final FieldDerivativeStructure<T> f = factory.build(data);
1802             for (int index = 0; index < factory.getCompiler().getFreeParameters(); ++index) {
1803                 final FieldDerivativeStructure<T> integ = f.differentiate(index, 0);
1804                 f.equals(integ);
1805             }
1806         }
1807     }
1808 
1809     @Test
1810     public void testIntegrationDifferentiation() {
1811         // check that integration and differentiation for univariate functions are each other inverse except for constant
1812         // term and highest order one
1813         final RandomGenerator random = new Well19937a(0x67fe66c05e5ee222l);
1814         final FDSFactory<T> factory = buildFactory(1, 25);
1815         final int size = factory.getCompiler().getSize();
1816         for (int count = 0; count < 100; ++count) {
1817             final double[] data = new double[size];
1818             for (int i = 1; i < size - 1; i++) {
1819                 data[i] = random.nextDouble();
1820             }
1821             final int indexVar = 0;
1822             final FieldDerivativeStructure<T> f = factory.build(data);
1823             final FieldDerivativeStructure<T> f2 = f.integrate(indexVar, 1).differentiate(indexVar, 1);
1824             final FieldDerivativeStructure<T> f3 = f.differentiate(indexVar, 1).integrate(indexVar, 1);
1825             checkEquals(f2, f, 0.);
1826             checkEquals(f2, f3, 0.);
1827             // check special case when non-positive integration order actually returns differentiation
1828             final FieldDerivativeStructure<T> df = f.integrate(indexVar, -1);
1829             final FieldDerivativeStructure<T> df2 = f.differentiate(indexVar, 1);
1830             checkEquals(df, df2, 0.);
1831             // check special case when non-positive differentiation order actually returns integration
1832             final FieldDerivativeStructure<T> fi  = f.differentiate(indexVar, -1);
1833             final FieldDerivativeStructure<T> fi2 = f.integrate(indexVar, 1);
1834             checkEquals(fi, fi2, 0.);
1835         }
1836     }
1837 
1838     @Test
1839     public void testDifferentiation1() {
1840         // check differentiation operator with result obtained manually
1841         final int freeParam = 3;
1842         final int order = 5;
1843         final FDSFactory<T> factory = buildFactory(freeParam, order);
1844         final FieldDerivativeStructure<T> f = factory.variable(0, 1.0);
1845         final int[] orders = new int[freeParam];
1846         orders[0] = 2;
1847         orders[1] = 1;
1848         orders[2] = 1;
1849         final T value = factory.getValueField().getZero().newInstance(10.);
1850         f.setDerivativeComponent(factory.getCompiler().getPartialDerivativeIndex(orders), value);
1851         final FieldDerivativeStructure<T> dfDx = f.differentiate(0, 1);
1852         orders[0] -= 1;
1853         assertEquals(1., dfDx.getPartialDerivative(new int[freeParam]).getReal(), 0.);
1854         assertEquals(value.getReal(), dfDx.getPartialDerivative(orders).getReal(), 0.);
1855         factory.constant(0.0).equals(f.differentiate(0, order + 1));
1856     }
1857 
1858     @Test
1859     public void testDifferentiation2() {
1860         // check that first-order differentiation twice is same as second-order differentiation
1861         final RandomGenerator random = new Well19937a(0xec293aaee352de94l);
1862         final FDSFactory<T> factory = buildFactory(5, 4);
1863         final int size = factory.getCompiler().getSize();
1864         for (int count = 0; count < 100; ++count) {
1865             final double[] data = new double[size];
1866             for (int i = 0; i < size; i++) {
1867                 data[i] = random.nextDouble();
1868             }
1869             final FieldDerivativeStructure<T> f = factory.build(data);
1870             final FieldDerivativeStructure<T> d2fDx2 = f.differentiate(0, 1).differentiate(0, 1);
1871             final FieldDerivativeStructure<T> d2fDx2Bis = f.differentiate(0, 2);
1872             checkEquals(d2fDx2, d2fDx2Bis, 0.);
1873         }
1874     }
1875 
1876     @Test
1877     public void testDifferentiation3() {
1878         // check that first-order differentiation on two variables does not depend on sequence of operations
1879         final RandomGenerator random = new Well19937a(0x35409ecc1348e46cl);
1880         final FDSFactory<T> factory = buildFactory(3, 7);
1881         final int size = factory.getCompiler().getSize();
1882         for (int count = 0; count < 100; ++count) {
1883             final double[] data = new double[size];
1884             for (int i = 0; i < size; i++) {
1885                 data[i] = random.nextDouble();
1886             }
1887             final FieldDerivativeStructure<T> f = factory.build(data);
1888             final FieldDerivativeStructure<T> d2fDxDy = f.differentiate(0, 1).differentiate(1, 1);
1889             final FieldDerivativeStructure<T> d2fDyDx = f.differentiate(1, 1).differentiate(0, 1);
1890             checkEquals(d2fDxDy, d2fDyDx, 0.);
1891         }
1892     }
1893 
1894     @Test
1895     public void testField() {
1896         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
1897             final FDSFactory<T> factory = buildFactory(3, maxOrder);
1898             FieldDerivativeStructure<T> x = factory.variable(0, 1.0);
1899             checkF0F1(x.getField().getZero(), 0.0, 0.0, 0.0, 0.0);
1900             checkF0F1(x.getField().getOne(), 1.0, 0.0, 0.0, 0.0);
1901             assertEquals(maxOrder, x.getField().getZero().getOrder());
1902             assertEquals(3, x.getField().getZero().getFreeParameters());
1903             assertEquals(FieldDerivativeStructure.class, x.getField().getRuntimeClass());
1904         }
1905     }
1906 
1907     @Test
1908     public void testEquals() {
1909         final FDSFactory<T> factory33 = buildFactory(3, 3);
1910         final FDSFactory<T> factory35 = buildFactory(3, 5);
1911         final FDSFactory<T> factory55 = buildFactory(5, 5);
1912         Assertions.assertFalse(factory33.variable(0, 1.0).equals(factory35.variable(0, 1.0)));
1913         Assertions.assertFalse(factory33.variable(0, 1.0).equals(factory55.variable(0, 1.0)));
1914         Assertions.assertTrue(factory33.variable(0, 1.0).equals(factory33.variable(0, 1.0)));
1915         Assertions.assertFalse(factory33.variable(0, 1.0).equals(new UnivariateDerivative1(0, 0)));
1916         Assertions.assertFalse(factory33.variable(0, 1.0).equals(factory33.variable(1, 1.0)));
1917         final FieldDerivativeStructure<T> x = factory33.variable(0, 1.0);
1918         Assertions.assertEquals(x, x);
1919         Assertions.assertNotEquals(0, x.hashCode());
1920     }
1921 
1922     @Test
1923     public void testOneParameterConstructor() {
1924         double x = 1.2;
1925         double cos = FastMath.cos(x);
1926         double sin = FastMath.sin(x);
1927         final FDSFactory<T> factory = buildFactory(1, 4);
1928         FieldDerivativeStructure<T> yRef = factory.variable(0, x).cos();
1929         try {
1930             factory.build(0.0, 0.0);
1931             fail("an exception should have been thrown");
1932         } catch (MathIllegalArgumentException dme) {
1933             // expected
1934         } catch (Exception e) {
1935             fail("wrong exceptionc caught " + e.getClass().getName());
1936         }
1937         double[] derivatives = new double[] { cos, -sin, -cos, sin, cos };
1938         FieldDerivativeStructure<T> y = factory.build(derivatives);
1939         checkEquals(yRef, y, 1.0e-15);
1940         T[] all = y.getAllDerivatives();
1941         assertEquals(derivatives.length, all.length);
1942         for (int i = 0; i < all.length; ++i) {
1943             assertEquals(derivatives[i], all[i].getReal(), 1.0e-15);
1944         }
1945     }
1946 
1947     @Test
1948     public void testOneOrderConstructor() {
1949         double x =  1.2;
1950         double y =  2.4;
1951         double z = 12.5;
1952         final FDSFactory<T> factory = buildFactory(3, 1);
1953         FieldDerivativeStructure<T> xRef = factory.variable(0, x);
1954         FieldDerivativeStructure<T> yRef = factory.variable(1, y);
1955         FieldDerivativeStructure<T> zRef = factory.variable(2, z);
1956         try {
1957             factory.build(x + y - z, 1.0, 1.0);
1958             fail("an exception should have been thrown");
1959         } catch (MathIllegalArgumentException dme) {
1960             // expected
1961         } catch (Exception e) {
1962             fail("wrong exceptionc caught " + e.getClass().getName());
1963         }
1964         double[] derivatives = new double[] { x + y - z, 1.0, 1.0, -1.0 };
1965         FieldDerivativeStructure<T> t = factory.build(derivatives);
1966         checkEquals(xRef.add(yRef.subtract(zRef)), t, 1.0e-15);
1967         T[] all = xRef.add(yRef.subtract(zRef)).getAllDerivatives();
1968         assertEquals(derivatives.length, all.length);
1969         for (int i = 0; i < all.length; ++i) {
1970             assertEquals(derivatives[i], all[i].getReal(), 1.0e-15);
1971         }
1972     }
1973 
1974     @Test
1975     public void testLinearCombination1DSDS() {
1976         doTestLinearCombination1DSDS(1.0e-15);
1977     }
1978 
1979     protected void doTestLinearCombination1DSDS(final double tol) {
1980         final FDSFactory<T> factory = buildFactory(6, 1);
1981         final FieldDerivativeStructure<T>[] a = MathArrays.buildArray(factory.getDerivativeField(), 3);
1982         a[0] = factory.variable(0, -1321008684645961.0 / 268435456.0);
1983         a[1] = factory.variable(1, -5774608829631843.0 / 268435456.0);
1984         a[2] = factory.variable(2, -7645843051051357.0 / 8589934592.0);
1985         final FieldDerivativeStructure<T>[] b = MathArrays.buildArray(factory.getDerivativeField(), 3);
1986         b[0] = factory.variable(3, -5712344449280879.0 / 2097152.0);
1987         b[1] = factory.variable(4, -4550117129121957.0 / 2097152.0);
1988         b[2] = factory.variable(5, 8846951984510141.0 / 131072.0);
1989 
1990         final FieldDerivativeStructure<T> abSumInline = a[0].linearCombination(a[0], b[0], a[1], b[1], a[2], b[2]);
1991         final FieldDerivativeStructure<T> abSumArray = a[0].linearCombination(a, b);
1992 
1993         assertEquals(abSumInline.getReal(), abSumArray.getReal(), 0);
1994         assertEquals(-1.8551294182586248737720779899, abSumInline.getReal(), tol);
1995         assertEquals(b[0].getReal(), abSumInline.getPartialDerivative(1, 0, 0, 0, 0, 0).getReal(), 1.0e-15);
1996         assertEquals(b[1].getReal(), abSumInline.getPartialDerivative(0, 1, 0, 0, 0, 0).getReal(), 1.0e-15);
1997         assertEquals(b[2].getReal(), abSumInline.getPartialDerivative(0, 0, 1, 0, 0, 0).getReal(), 1.0e-15);
1998         assertEquals(a[0].getReal(), abSumInline.getPartialDerivative(0, 0, 0, 1, 0, 0).getReal(), 1.0e-15);
1999         assertEquals(a[1].getReal(), abSumInline.getPartialDerivative(0, 0, 0, 0, 1, 0).getReal(), 1.0e-15);
2000         assertEquals(a[2].getReal(), abSumInline.getPartialDerivative(0, 0, 0, 0, 0, 1).getReal(), 1.0e-15);
2001 
2002     }
2003 
2004     @Test
2005     public void testLinearCombination1FieldDS() {
2006         doTestLinearCombination1FieldDS(1.0e-15);
2007     }
2008 
2009     protected void doTestLinearCombination1FieldDS(final double tol) {
2010         final FDSFactory<T> factory = buildFactory(3, 1);
2011         final T[] a = MathArrays.buildArray(getField(), 3);
2012         a[0] = buildScalar(-1321008684645961.0 / 268435456.0);
2013         a[1] = buildScalar(-5774608829631843.0 / 268435456.0);
2014         a[2] = buildScalar(-7645843051051357.0 / 8589934592.0);
2015         final FieldDerivativeStructure<T>[] b = MathArrays.buildArray(factory.getDerivativeField(), 3);
2016         b[0] = factory.variable(0, -5712344449280879.0 / 2097152.0);
2017         b[1] = factory.variable(1, -4550117129121957.0 / 2097152.0);
2018         b[2] = factory.variable(2, 8846951984510141.0 / 131072.0);
2019 
2020         final FieldDerivativeStructure<T> abSumInline = b[0].linearCombination(a[0], b[0],
2021                                                                                a[1], b[1],
2022                                                                                a[2], b[2]);
2023         final FieldDerivativeStructure<T> abSumArray = b[0].linearCombination(a, b);
2024 
2025         assertEquals(abSumInline.getReal(), abSumArray.getReal(), 0);
2026         assertEquals(-1.8551294182586248737720779899, abSumInline.getReal(), tol);
2027         assertEquals(a[0].getReal(), abSumInline.getPartialDerivative(1, 0, 0).getReal(), 1.0e-15);
2028         assertEquals(a[1].getReal(), abSumInline.getPartialDerivative(0, 1, 0).getReal(), 1.0e-15);
2029         assertEquals(a[2].getReal(), abSumInline.getPartialDerivative(0, 0, 1).getReal(), 1.0e-15);
2030 
2031     }
2032 
2033     @Test
2034     public void testLinearCombination1DoubleDS() {
2035         doTestLinearCombination1DoubleDS(1.0e-15);
2036     }
2037 
2038     protected void doTestLinearCombination1DoubleDS(final double tol) {
2039         final FDSFactory<T> factory = buildFactory(3, 1);
2040         final double[] a = new double[] {
2041             -1321008684645961.0 / 268435456.0,
2042             -5774608829631843.0 / 268435456.0,
2043             -7645843051051357.0 / 8589934592.0
2044         };
2045         final FieldDerivativeStructure<T>[] b = MathArrays.buildArray(factory.getDerivativeField(), 3);
2046         b[0] = factory.variable(0, -5712344449280879.0 / 2097152.0);
2047         b[1] = factory.variable(1, -4550117129121957.0 / 2097152.0);
2048         b[2] = factory.variable(2, 8846951984510141.0 / 131072.0);
2049 
2050         final FieldDerivativeStructure<T> abSumInline = b[0].linearCombination(a[0], b[0],
2051                                                                        a[1], b[1],
2052                                                                        a[2], b[2]);
2053         final FieldDerivativeStructure<T> abSumArray = b[0].linearCombination(a, b);
2054 
2055         assertEquals(abSumInline.getReal(), abSumArray.getReal(), 0);
2056         assertEquals(-1.8551294182586248737720779899, abSumInline.getReal(), tol);
2057         assertEquals(a[0], abSumInline.getPartialDerivative(1, 0, 0).getReal(), 1.0e-15);
2058         assertEquals(a[1], abSumInline.getPartialDerivative(0, 1, 0).getReal(), 1.0e-15);
2059         assertEquals(a[2], abSumInline.getPartialDerivative(0, 0, 1).getReal(), 1.0e-15);
2060 
2061     }
2062 
2063     @Test
2064     public void testLinearCombination2DSDS() {
2065 
2066         final FDSFactory<T> factory = buildFactory(4, 1);
2067         final FieldDerivativeStructure<T>[] u = MathArrays.buildArray(factory.getDerivativeField(), 4);
2068         final FieldDerivativeStructure<T>[] v = MathArrays.buildArray(factory.getDerivativeField(), 4);
2069 
2070         // we compare accurate versus naive dot product implementations
2071         // on regular vectors (i.e. not extreme cases like in the previous test)
2072         Well1024a random = new Well1024a(0xc6af886975069f11l);
2073 
2074         for (int i = 0; i < 10000; ++i) {
2075             for (int j = 0; j < u.length; ++j) {
2076                 u[j] = factory.variable(j, 1e17 * random.nextDouble());
2077                 v[j] = factory.constant(1e17 * random.nextDouble());
2078             }
2079 
2080             FieldDerivativeStructure<T> lin = u[0].linearCombination(u[0], v[0], u[1], v[1]);
2081             double ref = u[0].getReal() * v[0].getReal() +
2082                          u[1].getReal() * v[1].getReal();
2083             assertEquals(ref, lin.getReal(), 1.0e-15 * FastMath.abs(ref));
2084             assertEquals(v[0].getReal(), lin.getPartialDerivative(1, 0, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[0].getReal()));
2085             assertEquals(v[1].getReal(), lin.getPartialDerivative(0, 1, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[1].getReal()));
2086 
2087             lin = u[0].linearCombination(u[0], v[0], u[1], v[1], u[2], v[2]);
2088             ref = u[0].getReal() * v[0].getReal() +
2089                   u[1].getReal() * v[1].getReal() +
2090                   u[2].getReal() * v[2].getReal();
2091             assertEquals(ref, lin.getReal(), 1.0e-15 * FastMath.abs(ref));
2092             assertEquals(v[0].getReal(), lin.getPartialDerivative(1, 0, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[0].getReal()));
2093             assertEquals(v[1].getReal(), lin.getPartialDerivative(0, 1, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[1].getReal()));
2094             assertEquals(v[2].getReal(), lin.getPartialDerivative(0, 0, 1, 0).getReal(), 1.0e-15 * FastMath.abs(v[2].getReal()));
2095 
2096             lin = u[0].linearCombination(u[0], v[0], u[1], v[1], u[2], v[2], u[3], v[3]);
2097             ref = u[0].getReal() * v[0].getReal() +
2098                   u[1].getReal() * v[1].getReal() +
2099                   u[2].getReal() * v[2].getReal() +
2100                   u[3].getReal() * v[3].getReal();
2101             assertEquals(ref, lin.getReal(), 1.0e-15 * FastMath.abs(ref));
2102             assertEquals(v[0].getReal(), lin.getPartialDerivative(1, 0, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[0].getReal()));
2103             assertEquals(v[1].getReal(), lin.getPartialDerivative(0, 1, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[1].getReal()));
2104             assertEquals(v[2].getReal(), lin.getPartialDerivative(0, 0, 1, 0).getReal(), 1.0e-15 * FastMath.abs(v[2].getReal()));
2105             assertEquals(v[3].getReal(), lin.getPartialDerivative(0, 0, 0, 1).getReal(), 1.0e-15 * FastMath.abs(v[3].getReal()));
2106 
2107         }
2108     }
2109 
2110     @Test
2111     public void testLinearCombination2DoubleDS() {
2112         final FDSFactory<T> factory = buildFactory(4, 1);
2113         final double[] u = new double[4];
2114         final FieldDerivativeStructure<T>[] v = MathArrays.buildArray(factory.getDerivativeField(), 4);
2115         // we compare accurate versus naive dot product implementations
2116         // on regular vectors (i.e. not extreme cases like in the previous test)
2117         Well1024a random = new Well1024a(0xc6af886975069f11l);
2118 
2119         for (int i = 0; i < 10000; ++i) {
2120             for (int j = 0; j < u.length; ++j) {
2121                 u[j] = 1e17 * random.nextDouble();
2122                 v[j] = factory.variable(j, 1e17 * random.nextDouble());
2123             }
2124 
2125             FieldDerivativeStructure<T> lin = v[0].linearCombination(u[0], v[0], u[1], v[1]);
2126             double ref = u[0] * v[0].getReal() +
2127                          u[1] * v[1].getReal();
2128             assertEquals(ref, lin.getReal(), 1.0e-15 * FastMath.abs(ref));
2129             assertEquals(u[0], lin.getPartialDerivative(1, 0, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[0].getReal()));
2130             assertEquals(u[1], lin.getPartialDerivative(0, 1, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[1].getReal()));
2131 
2132             lin = v[0].linearCombination(u[0], v[0], u[1], v[1], u[2], v[2]);
2133             ref = u[0] * v[0].getReal() +
2134                   u[1] * v[1].getReal() +
2135                   u[2] * v[2].getReal();
2136             assertEquals(ref, lin.getReal(), 1.0e-15 * FastMath.abs(ref));
2137             assertEquals(u[0], lin.getPartialDerivative(1, 0, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[0].getReal()));
2138             assertEquals(u[1], lin.getPartialDerivative(0, 1, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[1].getReal()));
2139             assertEquals(u[2], lin.getPartialDerivative(0, 0, 1, 0).getReal(), 1.0e-15 * FastMath.abs(v[2].getReal()));
2140 
2141             lin = v[0].linearCombination(u[0], v[0], u[1], v[1], u[2], v[2], u[3], v[3]);
2142             ref = u[0] * v[0].getReal() +
2143                   u[1] * v[1].getReal() +
2144                   u[2] * v[2].getReal() +
2145                   u[3] * v[3].getReal();
2146             assertEquals(ref, lin.getReal(), 1.0e-15 * FastMath.abs(ref));
2147             assertEquals(u[0], lin.getPartialDerivative(1, 0, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[0].getReal()));
2148             assertEquals(u[1], lin.getPartialDerivative(0, 1, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[1].getReal()));
2149             assertEquals(u[2], lin.getPartialDerivative(0, 0, 1, 0).getReal(), 1.0e-15 * FastMath.abs(v[2].getReal()));
2150             assertEquals(u[3], lin.getPartialDerivative(0, 0, 0, 1).getReal(), 1.0e-15 * FastMath.abs(v[3].getReal()));
2151 
2152         }
2153     }
2154 
2155     @Test
2156     public void testLinearCombination2FieldDS() {
2157         final FDSFactory<T> factory = buildFactory(4, 1);
2158         final T[] u = MathArrays.buildArray(getField(), 4);
2159         final FieldDerivativeStructure<T>[] v = MathArrays.buildArray(factory.getDerivativeField(), 4);
2160         // we compare accurate versus naive dot product implementations
2161         // on regular vectors (i.e. not extreme cases like in the previous test)
2162         Well1024a random = new Well1024a(0xc6af886975069f11l);
2163 
2164         for (int i = 0; i < 10000; ++i) {
2165             for (int j = 0; j < u.length; ++j) {
2166                 u[j] = buildScalar(1e17 * random.nextDouble());
2167                 v[j] = factory.variable(j, 1e17 * random.nextDouble());
2168             }
2169 
2170             FieldDerivativeStructure<T> lin = v[0].linearCombination(u[0], v[0], u[1], v[1]);
2171             double ref = u[0].getReal() * v[0].getReal() +
2172                          u[1].getReal() * v[1].getReal();
2173             assertEquals(ref, lin.getReal(), 1.0e-15 * FastMath.abs(ref));
2174             assertEquals(u[0].getReal(), lin.getPartialDerivative(1, 0, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[0].getReal()));
2175             assertEquals(u[1].getReal(), lin.getPartialDerivative(0, 1, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[1].getReal()));
2176 
2177             lin = v[0].linearCombination(u[0], v[0], u[1], v[1], u[2], v[2]);
2178             ref = u[0].getReal() * v[0].getReal() +
2179                   u[1].getReal() * v[1].getReal() +
2180                   u[2].getReal() * v[2].getReal();
2181             assertEquals(ref, lin.getReal(), 1.0e-15 * FastMath.abs(ref));
2182             assertEquals(u[0].getReal(), lin.getPartialDerivative(1, 0, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[0].getReal()));
2183             assertEquals(u[1].getReal(), lin.getPartialDerivative(0, 1, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[1].getReal()));
2184             assertEquals(u[2].getReal(), lin.getPartialDerivative(0, 0, 1, 0).getReal(), 1.0e-15 * FastMath.abs(v[2].getReal()));
2185 
2186             lin = v[0].linearCombination(u[0], v[0], u[1], v[1], u[2], v[2], u[3], v[3]);
2187             ref = u[0].getReal() * v[0].getReal() +
2188                   u[1].getReal() * v[1].getReal() +
2189                   u[2].getReal() * v[2].getReal() +
2190                   u[3].getReal() * v[3].getReal();
2191             assertEquals(ref, lin.getReal(), 1.0e-15 * FastMath.abs(ref));
2192             assertEquals(u[0].getReal(), lin.getPartialDerivative(1, 0, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[0].getReal()));
2193             assertEquals(u[1].getReal(), lin.getPartialDerivative(0, 1, 0, 0).getReal(), 1.0e-15 * FastMath.abs(v[1].getReal()));
2194             assertEquals(u[2].getReal(), lin.getPartialDerivative(0, 0, 1, 0).getReal(), 1.0e-15 * FastMath.abs(v[2].getReal()));
2195             assertEquals(u[3].getReal(), lin.getPartialDerivative(0, 0, 0, 1).getReal(), 1.0e-15 * FastMath.abs(v[3].getReal()));
2196 
2197         }
2198     }
2199 
2200     @Test
2201     public void testZero() {
2202         FDSFactory<T> factory = buildFactory(3, 2);
2203         FieldDerivativeStructure<T> zero = factory.constant(17).getField().getZero();
2204         T[] a = zero.getAllDerivatives();
2205         assertEquals(10, a.length);
2206         for (int i = 0; i < a.length; ++i) {
2207             assertEquals(buildScalar(0.0), a[i]);
2208         }
2209     }
2210 
2211     @Test
2212     public void testOne() {
2213         FDSFactory<T> factory = buildFactory(3, 2);
2214         FieldDerivativeStructure<T> one = factory.constant(17).getField().getOne();
2215         T[] a = one.getAllDerivatives();
2216         assertEquals(10, a.length);
2217         for (int i = 0; i < a.length; ++i) {
2218             assertEquals(i == 0 ? buildScalar(1.0) : buildScalar(0.0), a[i]);
2219         }
2220     }
2221 
2222     @Test
2223     public void testMap() {
2224         List<int[]> pairs = new ArrayList<>();
2225         for (int parameters = 1; parameters < 5; ++parameters) {
2226             for (int order = 0; order < 3; ++order) {
2227                 pairs.add(new int[] { parameters, order });
2228             }
2229         }
2230         Map<Field<?>, Integer> map = new HashMap<>();
2231         for (int i = 0; i < 1000; ++i) {
2232             // create a brand new factory for each derivative
2233             int parameters = pairs.get(i % pairs.size())[0];
2234             int order      = pairs.get(i % pairs.size())[1];
2235             FDSFactory<T> factory = buildFactory(parameters, order);
2236             map.put(factory.constant(buildScalar(17)).getField(), 0);
2237         }
2238 
2239         // despite we have created numerous factories,
2240         // there should be only one field for each pair parameters/order
2241         assertEquals(pairs.size(), map.size());
2242         @SuppressWarnings("unchecked")
2243         Field<FieldDerivativeStructure<T>> first = (Field<FieldDerivativeStructure<T>>) map.entrySet().iterator().next().getKey();
2244         assertEquals(first, first);
2245         assertNotEquals(first, getField());
2246 
2247         // even at same parameters and differentiation orders, different values generate different fields
2248         FieldDerivativeStructure<T> zero64 = buildFactory(3, 2).build();
2249         FieldDerivativeStructure<Dfp> zeroDFP = new FDSFactory<Dfp>(new DfpField(15), 3, 2).build();
2250         assertEquals(zero64.getFreeParameters(), zeroDFP.getFreeParameters());
2251         assertEquals(zero64.getOrder(), zeroDFP.getOrder());
2252         assertNotEquals(zero64.getField(), zeroDFP.getField());
2253     }
2254 
2255     @SuppressWarnings("unchecked")
2256     @Test
2257     public void testRebaseConditions() {
2258         final FDSFactory<T> f32 = buildFactory(3, 2);
2259         final FDSFactory<T> f22 = buildFactory(2, 2);
2260         final FDSFactory<T> f31 = buildFactory(3, 1);
2261         try {
2262             f32.variable(0, 0).rebase(f22.variable(0, 0), f22.variable(1, 1.0));
2263         } catch (MathIllegalArgumentException miae) {
2264             assertEquals(LocalizedCoreFormats.DIMENSIONS_MISMATCH, miae.getSpecifier());
2265             assertEquals(3, ((Integer) miae.getParts()[0]).intValue());
2266             assertEquals(2, ((Integer) miae.getParts()[1]).intValue());
2267         }
2268         try {
2269             f32.variable(0, 0).rebase(f31.variable(0, 0), f31.variable(1, 1.0), f31.variable(2, 2.0));
2270         } catch (MathIllegalArgumentException miae) {
2271             assertEquals(LocalizedCoreFormats.DIMENSIONS_MISMATCH, miae.getSpecifier());
2272             assertEquals(2, ((Integer) miae.getParts()[0]).intValue());
2273             assertEquals(1, ((Integer) miae.getParts()[1]).intValue());
2274         }
2275     }
2276 
2277     @SuppressWarnings("unchecked")
2278     @Test
2279     public void testRebaseNoVariables() {
2280         final FieldDerivativeStructure<T> x = buildFactory(0, 2).constant(1.0);
2281         assertSame(x, x.rebase());
2282     }
2283 
2284     @Test
2285     public void testRebaseValueMoreIntermediateThanBase() {
2286         doTestRebaseValue(createBaseVariables(buildFactory(2, 4), 1.5, -2.0),
2287                           q -> {
2288                               final FieldDerivativeStructure<T>[] a = MathArrays.buildArray(q[0].getFactory().getDerivativeField(), 3);
2289                               a[0] = q[0].add(q[1].multiply(3));
2290                               a[1] = q[0].log();
2291                               a[2] = q[1].divide(q[0].sin());
2292                               return a;
2293                           },
2294                           buildFactory(3, 4),
2295                           p -> p[0].add(p[1].divide(p[2])),
2296                           1.0e-15);
2297     }
2298 
2299     @Test
2300     public void testRebaseValueLessIntermediateThanBase() {
2301         doTestRebaseValue(createBaseVariables(buildFactory(3, 4), 1.5, -2.0, 0.5),
2302                           q -> {
2303                               final FieldDerivativeStructure<T>[] a = MathArrays.buildArray(q[0].getFactory().getDerivativeField(), 2);
2304                               a[0] = q[0].add(q[1].multiply(3));
2305                               a[1] = q[0].add(q[1]).subtract(q[2]);
2306                               return a;
2307                           },
2308                           buildFactory(2, 4),
2309                           p -> p[0].multiply(p[1]),
2310                           1.0e-15);
2311     }
2312 
2313     @Test
2314     public void testRebaseValueEqualIntermediateAndBase() {
2315         doTestRebaseValue(createBaseVariables(buildFactory(2, 4), 1.5, -2.0),
2316                           q -> {
2317                               final FieldDerivativeStructure<T>[] a = MathArrays.buildArray(q[0].getFactory().getDerivativeField(), 2);
2318                               a[0] = q[0].add(q[1].multiply(3));
2319                               a[1] = q[0].add(q[1]);
2320                               return a;
2321                           },
2322                           buildFactory(2, 4),
2323                           p -> p[0].multiply(p[1]),
2324                           1.0e-15);
2325     }
2326 
2327     private void doTestRebaseValue(final FieldDerivativeStructure<T>[] q,
2328                                    final CalculusFieldMultivariateVectorFunction<FieldDerivativeStructure<T>> qToP,
2329                                    final FDSFactory<T> factoryP,
2330                                    final CalculusFieldMultivariateFunction<FieldDerivativeStructure<T>> f,
2331                                    final double tol) {
2332 
2333         // intermediate variables as functions of base variables
2334         final FieldDerivativeStructure<T>[] pBase = qToP.value(q);
2335         
2336         // reference function
2337         final FieldDerivativeStructure<T> ref = f.value(pBase);
2338 
2339         // intermediate variables as independent variables
2340         final FieldDerivativeStructure<T>[] pIntermediate = creatIntermediateVariables(factoryP, pBase);
2341 
2342         // function of the intermediate variables
2343         final FieldDerivativeStructure<T> fI = f.value(pIntermediate);
2344 
2345         // function rebased to base variables
2346         final FieldDerivativeStructure<T> rebased = fI.rebase(pBase);
2347 
2348         assertEquals(q[0].getFreeParameters(),                   ref.getFreeParameters());
2349         assertEquals(q[0].getOrder(),                            ref.getOrder());
2350         assertEquals(factoryP.getCompiler().getFreeParameters(), fI.getFreeParameters());
2351         assertEquals(factoryP.getCompiler().getOrder(),          fI.getOrder());
2352         assertEquals(ref.getFreeParameters(),                    rebased.getFreeParameters());
2353         assertEquals(ref.getOrder(),                             rebased.getOrder());
2354 
2355         checkEquals(ref, rebased, tol);
2356 
2357     }
2358 
2359     final FieldDerivativeStructure<T>[] createBaseVariables(final FDSFactory<T> factory, double... q) {
2360         final FieldDerivativeStructure<T>[] qDS = MathArrays.buildArray(factory.getDerivativeField(), q.length);
2361         for (int i = 0; i < q.length; ++i) {
2362             qDS[i] = factory.variable(i, q[i]);
2363         }
2364         return qDS;
2365     }
2366 
2367     final FieldDerivativeStructure<T>[] creatIntermediateVariables(final FDSFactory<T> factory,
2368                                                                    @SuppressWarnings("unchecked") FieldDerivativeStructure<T>... pBase) {
2369         final FieldDerivativeStructure<T>[] pIntermediate = MathArrays.buildArray(factory.getDerivativeField(), pBase.length);
2370         for (int i = 0; i < pBase.length; ++i) {
2371             pIntermediate[i] = factory.variable(i, pBase[i].getValue());
2372         }
2373         return pIntermediate;
2374     }
2375 
2376     @Test
2377     public void testRunTimeClass() {
2378         FDSFactory<T> factory = buildFactory(3, 2);
2379         Field<FieldDerivativeStructure<T>> field = factory.getDerivativeField();
2380         assertEquals(FieldDerivativeStructure.class, field.getRuntimeClass());
2381         assertEquals(getField(), factory.getValueField());
2382         assertEquals("org.hipparchus.analysis.differentiation.FDSFactory$DerivativeField",
2383                             factory.getDerivativeField().getClass().getName());
2384     }
2385 
2386     private void checkF0F1(FieldDerivativeStructure<T> ds, double value, double...derivatives) {
2387 
2388         // check dimension
2389         assertEquals(derivatives.length, ds.getFreeParameters());
2390 
2391         // check value, directly and also as 0th order derivative
2392         assertEquals(value, ds.getReal(), 1.0e-15);
2393         assertEquals(value, ds.getPartialDerivative(new int[ds.getFreeParameters()]).getReal(), 1.0e-15);
2394 
2395         // check first order derivatives
2396         for (int i = 0; i < derivatives.length; ++i) {
2397             int[] orders = new int[derivatives.length];
2398             orders[i] = 1;
2399             assertEquals(derivatives[i], ds.getPartialDerivative(orders).getReal(), 1.0e-15);
2400         }
2401 
2402     }
2403 
2404     public static <T extends CalculusFieldElement<T>> void checkEquals(FieldDerivativeStructure<T> ds1,
2405                                                                        FieldDerivativeStructure<T> ds2,
2406                                                                        double epsilon) {
2407 
2408         // check dimension
2409         assertEquals(ds1.getFreeParameters(), ds2.getFreeParameters());
2410         assertEquals(ds1.getOrder(), ds2.getOrder());
2411 
2412         int[] derivatives = new int[ds1.getFreeParameters()];
2413         int sum = 0;
2414         while (true) {
2415 
2416             if (sum <= ds1.getOrder()) {
2417                 assertEquals(ds1.getPartialDerivative(derivatives).getReal(),
2418                                     ds2.getPartialDerivative(derivatives).getReal(),
2419                                     epsilon);
2420             }
2421 
2422             boolean increment = true;
2423             sum = 0;
2424             for (int i = derivatives.length - 1; i >= 0; --i) {
2425                 if (increment) {
2426                     if (derivatives[i] == ds1.getOrder()) {
2427                         derivatives[i] = 0;
2428                     } else {
2429                         derivatives[i]++;
2430                         increment = false;
2431                     }
2432                 }
2433                 sum += derivatives[i];
2434             }
2435             if (increment) {
2436                 return;
2437             }
2438 
2439         }
2440 
2441     }
2442 
2443 }