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