View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://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.apache.commons.math4.legacy.analysis.differentiation;
19  
20  import java.util.Arrays;
21  import java.util.List;
22  
23  import org.apache.commons.math4.legacy.field.ExtendedFieldElementAbstractTest;
24  import org.apache.commons.math4.legacy.TestUtils;
25  import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction;
26  import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
27  import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
28  import org.apache.commons.rng.UniformRandomProvider;
29  import org.apache.commons.rng.simple.RandomSource;
30  import org.apache.commons.numbers.core.ArithmeticUtils;
31  import org.apache.commons.numbers.combinatorics.Factorial;
32  import org.apache.commons.math4.core.jdkmath.JdkMath;
33  import org.apache.commons.numbers.core.Precision;
34  import org.junit.Assert;
35  import org.junit.Test;
36  
37  /**
38   * Test for class {@link DerivativeStructure}.
39   */
40  public class DerivativeStructureTest extends ExtendedFieldElementAbstractTest<DerivativeStructure> {
41  
42      @Override
43      protected DerivativeStructure build(final double x) {
44          return new DerivativeStructure(2, 1, 0, x);
45      }
46  
47      @Test(expected=NumberIsTooLargeException.class)
48      public void testWrongVariableIndex() {
49          new DerivativeStructure(3, 1, 3, 1.0);
50      }
51  
52      @Test(expected=DimensionMismatchException.class)
53      public void testMissingOrders() {
54          new DerivativeStructure(3, 1, 0, 1.0).getPartialDerivative(0, 1);
55      }
56  
57      @Test(expected=NumberIsTooLargeException.class)
58      public void testTooLargeOrder() {
59          new DerivativeStructure(3, 1, 0, 1.0).getPartialDerivative(1, 1, 2);
60      }
61  
62      @Test
63      public void testVariableWithoutDerivative0() {
64          DerivativeStructure v = new DerivativeStructure(1, 0, 0, 1.0);
65          Assert.assertEquals(1.0, v.getValue(), 1.0e-15);
66      }
67  
68      @Test(expected=NumberIsTooLargeException.class)
69      public void testVariableWithoutDerivative1() {
70          DerivativeStructure v = new DerivativeStructure(1, 0, 0, 1.0);
71          Assert.assertEquals(1.0, v.getPartialDerivative(1), 1.0e-15);
72      }
73  
74      @Test
75      public void testVariable() {
76          for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
77              checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0),
78                        1.0, 1.0, 0.0, 0.0);
79              checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0),
80                        2.0, 0.0, 1.0, 0.0);
81              checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0),
82                        3.0, 0.0, 0.0, 1.0);
83          }
84      }
85  
86      @Test
87      public void testConstant() {
88          for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
89              checkF0F1(new DerivativeStructure(3, maxOrder, JdkMath.PI),
90                        JdkMath.PI, 0.0, 0.0, 0.0);
91          }
92      }
93  
94      @Test
95      public void testCreateConstant() {
96          DerivativeStructure a = new DerivativeStructure(3, 2, 0, 1.3);
97          DerivativeStructure b = a.createConstant(2.5);
98          Assert.assertEquals(a.getFreeParameters(), b.getFreeParameters());
99          Assert.assertEquals(a.getOrder(), b.getOrder());
100         checkEquals(a.getField().getOne().multiply(2.5), b, 1.0e-15);
101     }
102 
103     @Test
104     public void testPrimitiveAdd() {
105         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
106             checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).add(5), 6.0, 1.0, 0.0, 0.0);
107             checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).add(5), 7.0, 0.0, 1.0, 0.0);
108             checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).add(5), 8.0, 0.0, 0.0, 1.0);
109         }
110     }
111 
112     @Test
113     public void testAdd() {
114         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
115             DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
116             DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
117             DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
118             DerivativeStructure xyz = x.add(y.add(z));
119             checkF0F1(xyz, x.getValue() + y.getValue() + z.getValue(), 1.0, 1.0, 1.0);
120         }
121     }
122 
123     @Test
124     public void testPrimitiveSubtract() {
125         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
126             checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).subtract(5), -4.0, 1.0, 0.0, 0.0);
127             checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).subtract(5), -3.0, 0.0, 1.0, 0.0);
128             checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).subtract(5), -2.0, 0.0, 0.0, 1.0);
129         }
130     }
131 
132     @Test
133     public void testSubtract() {
134         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
135             DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
136             DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
137             DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
138             DerivativeStructure xyz = x.subtract(y.subtract(z));
139             checkF0F1(xyz, x.getValue() - (y.getValue() - z.getValue()), 1.0, -1.0, 1.0);
140         }
141     }
142 
143     @Test
144     public void testPrimitiveMultiply() {
145         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
146             checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).multiply(5),  5.0, 5.0, 0.0, 0.0);
147             checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).multiply(5), 10.0, 0.0, 5.0, 0.0);
148             checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).multiply(5), 15.0, 0.0, 0.0, 5.0);
149         }
150     }
151 
152     @Test
153     public void testMultiply() {
154         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
155             DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
156             DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
157             DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
158             DerivativeStructure xyz = x.multiply(y.multiply(z));
159             for (int i = 0; i <= maxOrder; ++i) {
160                 for (int j = 0; j <= maxOrder; ++j) {
161                     for (int k = 0; k <= maxOrder; ++k) {
162                         if (i + j + k <= maxOrder) {
163                             Assert.assertEquals((i == 0 ? x.getValue() : (i == 1 ? 1.0 : 0.0)) *
164                                                 (j == 0 ? y.getValue() : (j == 1 ? 1.0 : 0.0)) *
165                                                 (k == 0 ? z.getValue() : (k == 1 ? 1.0 : 0.0)),
166                                                 xyz.getPartialDerivative(i, j, k),
167                                                 1.0e-15);
168                         }
169                     }
170                 }
171             }
172         }
173     }
174 
175     @Test
176     public void testNegate() {
177         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
178             checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).negate(), -1.0, -1.0, 0.0, 0.0);
179             checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).negate(), -2.0, 0.0, -1.0, 0.0);
180             checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).negate(), -3.0, 0.0, 0.0, -1.0);
181         }
182     }
183 
184     @Test
185     public void testReciprocal() {
186         for (double x = 0.1; x < 1.2; x += 0.1) {
187             DerivativeStructure r = new DerivativeStructure(1, 6, 0, x).reciprocal();
188             Assert.assertEquals(1 / x, r.getValue(), 1.0e-15);
189             for (int i = 1; i < r.getOrder(); ++i) {
190                 double expected = ArithmeticUtils.pow(-1, i) * Factorial.value(i) /
191                                   JdkMath.pow(x, i + 1);
192                 Assert.assertEquals(expected, r.getPartialDerivative(i), 1.0e-15 * JdkMath.abs(expected));
193             }
194         }
195     }
196 
197     @Test
198     public void testPow() {
199         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
200             for (int n = 0; n < 10; ++n) {
201 
202                 DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
203                 DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
204                 DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
205                 List<DerivativeStructure> list = Arrays.asList(x, y, z,
206                                                                x.add(y).add(z),
207                                                                x.multiply(y).multiply(z));
208 
209                 if (n == 0) {
210                     for (DerivativeStructure ds : list) {
211                         checkEquals(ds.getField().getOne(), ds.pow(n), 1.0e-15);
212                     }
213                 } else if (n == 1) {
214                     for (DerivativeStructure ds : list) {
215                         checkEquals(ds, ds.pow(n), 1.0e-15);
216                     }
217                 } else {
218                     for (DerivativeStructure ds : list) {
219                         DerivativeStructure p = ds.getField().getOne();
220                         for (int i = 0; i < n; ++i) {
221                             p = p.multiply(ds);
222                         }
223                         checkEquals(p, ds.pow(n), 1.0e-15);
224                     }
225                 }
226             }
227         }
228     }
229 
230     @Test
231     public void testPowDoubleDS() {
232         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
233 
234             DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 0.1);
235             DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 0.2);
236             DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 0.3);
237             List<DerivativeStructure> list = Arrays.asList(x, y, z,
238                                                            x.add(y).add(z),
239                                                            x.multiply(y).multiply(z));
240 
241             for (DerivativeStructure ds : list) {
242                 // the special case a = 0 is included here
243                 for (double a : new double[] { 0.0, 0.1, 1.0, 2.0, 5.0 }) {
244                     DerivativeStructure reference = (a == 0) ?
245                                                     x.getField().getZero() :
246                                                     new DerivativeStructure(3, maxOrder, a).pow(ds);
247                     DerivativeStructure result = DerivativeStructure.pow(a, ds);
248                     checkEquals(reference, result, 1.0e-15);
249                 }
250             }
251 
252             // negative base: -1^x can be evaluated for integers only, so value is sometimes OK, derivatives are always NaN
253             DerivativeStructure negEvenInteger = DerivativeStructure.pow(-2.0, new DerivativeStructure(3,  maxOrder, 0, 2.0));
254             Assert.assertEquals(4.0, negEvenInteger.getValue(), 1.0e-15);
255             Assert.assertTrue(Double.isNaN(negEvenInteger.getPartialDerivative(1, 0, 0)));
256             DerivativeStructure negOddInteger = DerivativeStructure.pow(-2.0, new DerivativeStructure(3,  maxOrder, 0, 3.0));
257             Assert.assertEquals(-8.0, negOddInteger.getValue(), 1.0e-15);
258             Assert.assertTrue(Double.isNaN(negOddInteger.getPartialDerivative(1, 0, 0)));
259             DerivativeStructure negNonInteger = DerivativeStructure.pow(-2.0, new DerivativeStructure(3,  maxOrder, 0, 2.001));
260             Assert.assertTrue(Double.isNaN(negNonInteger.getValue()));
261             Assert.assertTrue(Double.isNaN(negNonInteger.getPartialDerivative(1, 0, 0)));
262 
263             DerivativeStructure zeroNeg = DerivativeStructure.pow(0.0, new DerivativeStructure(3,  maxOrder, 0, -1.0));
264             Assert.assertTrue(Double.isNaN(zeroNeg.getValue()));
265             Assert.assertTrue(Double.isNaN(zeroNeg.getPartialDerivative(1, 0, 0)));
266             DerivativeStructure posNeg = DerivativeStructure.pow(2.0, new DerivativeStructure(3,  maxOrder, 0, -2.0));
267             Assert.assertEquals(1.0 / 4.0, posNeg.getValue(), 1.0e-15);
268             Assert.assertEquals(JdkMath.log(2.0) / 4.0, posNeg.getPartialDerivative(1, 0, 0), 1.0e-15);
269 
270             // very special case: a = 0 and power = 0
271             DerivativeStructure zeroZero = DerivativeStructure.pow(0.0, new DerivativeStructure(3,  maxOrder, 0, 0.0));
272 
273             // this should be OK for simple first derivative with one variable only ...
274             Assert.assertEquals(1.0, zeroZero.getValue(), 1.0e-15);
275             Assert.assertEquals(Double.NEGATIVE_INFINITY, zeroZero.getPartialDerivative(1, 0, 0), 1.0e-15);
276 
277             // the following checks show a LIMITATION of the current implementation
278             // we have no way to tell x is a pure linear variable x = 0
279             // we only say: "x is a structure with value = 0.0,
280             // first derivative with respect to x = 1.0, and all other derivatives
281             // (first order with respect to y and z and higher derivatives) all 0.0.
282             // We have function f(x) = a^x and x = 0 so we compute:
283             // f(0) = 1, f'(0) = ln(a), f''(0) = ln(a)^2. The limit of these values
284             // when a converges to 0 implies all derivatives keep switching between
285             // +infinity and -infinity.
286             //
287             // Function composition rule for first derivatives is:
288             // d[f(g(x,y,z))]/dy = f'(g(x,y,z)) * dg(x,y,z)/dy
289             // so given that in our case x represents g and does not depend
290             // on y or z, we have dg(x,y,z)/dy = 0
291             // applying the composition rules gives:
292             // d[f(g(x,y,z))]/dy = f'(g(x,y,z)) * dg(x,y,z)/dy
293             //                 = -infinity * 0
294             //                 = NaN
295             // if we knew x is really the x variable and not the identity
296             // function applied to x, we would not have computed f'(g(x,y,z)) * dg(x,y,z)/dy
297             // and we would have found that the result was 0 and not NaN
298             Assert.assertTrue(Double.isNaN(zeroZero.getPartialDerivative(0, 1, 0)));
299             Assert.assertTrue(Double.isNaN(zeroZero.getPartialDerivative(0, 0, 1)));
300 
301             // Function composition rule for second derivatives is:
302             // d2[f(g(x))]/dx2 = f''(g(x)) * [g'(x)]^2 + f'(g(x)) * g''(x)
303             // when function f is the a^x root and x = 0 we have:
304             // f(0) = 1, f'(0) = ln(a), f''(0) = ln(a)^2 which for a = 0 implies
305             // all derivatives keep switching between +infinity and -infinity
306             // so given that in our case x represents g, we have g(x) = 0,
307             // g'(x) = 1 and g''(x) = 0
308             // applying the composition rules gives:
309             // d2[f(g(x))]/dx2 = f''(g(x)) * [g'(x)]^2 + f'(g(x)) * g''(x)
310             //                 = +infinity * 1^2 + -infinity * 0
311             //                 = +infinity + NaN
312             //                 = NaN
313             // if we knew x is really the x variable and not the identity
314             // function applied to x, we would not have computed f'(g(x)) * g''(x)
315             // and we would have found that the result was +infinity and not NaN
316             if (maxOrder > 1) {
317                 Assert.assertTrue(Double.isNaN(zeroZero.getPartialDerivative(2, 0, 0)));
318                 Assert.assertTrue(Double.isNaN(zeroZero.getPartialDerivative(0, 2, 0)));
319                 Assert.assertTrue(Double.isNaN(zeroZero.getPartialDerivative(0, 0, 2)));
320                 Assert.assertTrue(Double.isNaN(zeroZero.getPartialDerivative(1, 1, 0)));
321                 Assert.assertTrue(Double.isNaN(zeroZero.getPartialDerivative(0, 1, 1)));
322                 Assert.assertTrue(Double.isNaN(zeroZero.getPartialDerivative(1, 1, 0)));
323             }
324 
325             // very special case: 0^0 where the power is a primitive
326             DerivativeStructure zeroDsZeroDouble = new DerivativeStructure(3,  maxOrder, 0, 0.0).pow(0.0);
327             boolean first = true;
328             for (final double d : zeroDsZeroDouble.getAllDerivatives()) {
329                 if (first) {
330                     Assert.assertEquals(1.0, d, Precision.EPSILON);
331                     first = false;
332                 } else {
333                     Assert.assertEquals(0.0, d, Precision.SAFE_MIN);
334                 }
335             }
336             DerivativeStructure zeroDsZeroInt = new DerivativeStructure(3,  maxOrder, 0, 0.0).pow(0);
337             first = true;
338             for (final double d : zeroDsZeroInt.getAllDerivatives()) {
339                 if (first) {
340                     Assert.assertEquals(1.0, d, Precision.EPSILON);
341                     first = false;
342                 } else {
343                     Assert.assertEquals(0.0, d, Precision.SAFE_MIN);
344                 }
345             }
346 
347             // 0^p with p smaller than 1.0
348             DerivativeStructure u = new DerivativeStructure(3, maxOrder, 1, -0.0).pow(0.25);
349             for (int i0 = 0; i0 <= maxOrder; ++i0) {
350                 for (int i1 = 0; i1 <= maxOrder; ++i1) {
351                     for (int i2 = 0; i2 <= maxOrder; ++i2) {
352                         if (i0 + i1 + i2 <= maxOrder) {
353                             Assert.assertEquals(0.0, u.getPartialDerivative(i0, i1, i2), 1.0e-10);
354                         }
355                     }
356                 }
357             }
358         }
359     }
360 
361     @Test
362     public void testExpression() {
363         double epsilon = 2.5e-13;
364         for (double x = 0; x < 2; x += 0.2) {
365             DerivativeStructure dsX = new DerivativeStructure(3, 5, 0, x);
366             for (double y = 0; y < 2; y += 0.2) {
367                 DerivativeStructure dsY = new DerivativeStructure(3, 5, 1, y);
368                 for (double z = 0; z >- 2; z -= 0.2) {
369                     DerivativeStructure dsZ = new DerivativeStructure(3, 5, 2, z);
370 
371                     // f(x, y, z) = x + 5 x y - 2 z + (8 z x - y)^3
372                     DerivativeStructure ds =
373                             new DerivativeStructure(1, dsX,
374                                                     5, dsX.multiply(dsY),
375                                                     -2, dsZ,
376                                                     1, new DerivativeStructure(8, dsZ.multiply(dsX),
377                                                                                -1, dsY).pow(3));
378                     DerivativeStructure dsOther =
379                             new DerivativeStructure(1, dsX,
380                                                     5, dsX.multiply(dsY),
381                                                     -2, dsZ).add(new DerivativeStructure(8, dsZ.multiply(dsX),
382                                                                                          -1, dsY).pow(3));
383                     double f = x + 5 * x * y - 2 * z + JdkMath.pow(8 * z * x - y, 3);
384                     Assert.assertEquals(f, ds.getValue(),
385                                         JdkMath.abs(epsilon * f));
386                     Assert.assertEquals(f, dsOther.getValue(),
387                                         JdkMath.abs(epsilon * f));
388 
389                     // df/dx = 1 + 5 y + 24 (8 z x - y)^2 z
390                     double dfdx = 1 + 5 * y + 24 * z * JdkMath.pow(8 * z * x - y, 2);
391                     Assert.assertEquals(dfdx, ds.getPartialDerivative(1, 0, 0),
392                                         JdkMath.abs(epsilon * dfdx));
393                     Assert.assertEquals(dfdx, dsOther.getPartialDerivative(1, 0, 0),
394                                         JdkMath.abs(epsilon * dfdx));
395 
396                     // df/dxdy = 5 + 48 z*(y - 8 z x)
397                     double dfdxdy = 5 + 48 * z * (y - 8 * z * x);
398                     Assert.assertEquals(dfdxdy, ds.getPartialDerivative(1, 1, 0),
399                                         JdkMath.abs(epsilon * dfdxdy));
400                     Assert.assertEquals(dfdxdy, dsOther.getPartialDerivative(1, 1, 0),
401                                         JdkMath.abs(epsilon * dfdxdy));
402 
403                     // df/dxdydz = 48 (y - 16 z x)
404                     double dfdxdydz = 48 * (y - 16 * z * x);
405                     Assert.assertEquals(dfdxdydz, ds.getPartialDerivative(1, 1, 1),
406                                         JdkMath.abs(epsilon * dfdxdydz));
407                     Assert.assertEquals(dfdxdydz, dsOther.getPartialDerivative(1, 1, 1),
408                                         JdkMath.abs(epsilon * dfdxdydz));
409                 }
410             }
411         }
412     }
413 
414     @Test
415     public void testCompositionOneVariableX() {
416         double epsilon = 1.0e-13;
417         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
418             for (double x = 0.1; x < 1.2; x += 0.1) {
419                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
420                 for (double y = 0.1; y < 1.2; y += 0.1) {
421                     DerivativeStructure dsY = new DerivativeStructure(1, maxOrder, y);
422                     DerivativeStructure f = dsX.divide(dsY).sqrt();
423                     double f0 = JdkMath.sqrt(x / y);
424                     Assert.assertEquals(f0, f.getValue(), JdkMath.abs(epsilon * f0));
425                     if (f.getOrder() > 0) {
426                         double f1 = 1 / (2 * JdkMath.sqrt(x * y));
427                         Assert.assertEquals(f1, f.getPartialDerivative(1), JdkMath.abs(epsilon * f1));
428                         if (f.getOrder() > 1) {
429                             double f2 = -f1 / (2 * x);
430                             Assert.assertEquals(f2, f.getPartialDerivative(2), JdkMath.abs(epsilon * f2));
431                             if (f.getOrder() > 2) {
432                                 double f3 = (f0 + x / (2 * y * f0)) / (4 * x * x * x);
433                                 Assert.assertEquals(f3, f.getPartialDerivative(3), JdkMath.abs(epsilon * f3));
434                             }
435                         }
436                     }
437                 }
438             }
439         }
440     }
441 
442     @Test
443     public void testTrigo() {
444         double epsilon = 2.0e-12;
445         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
446             for (double x = 0.1; x < 1.2; x += 0.1) {
447                 DerivativeStructure dsX = new DerivativeStructure(3, maxOrder, 0, x);
448                 for (double y = 0.1; y < 1.2; y += 0.1) {
449                     DerivativeStructure dsY = new DerivativeStructure(3, maxOrder, 1, y);
450                     for (double z = 0.1; z < 1.2; z += 0.1) {
451                         DerivativeStructure dsZ = new DerivativeStructure(3, maxOrder, 2, z);
452                         DerivativeStructure f = dsX.divide(dsY.cos().add(dsZ.tan())).sin();
453                         double a = JdkMath.cos(y) + JdkMath.tan(z);
454                         double f0 = JdkMath.sin(x / a);
455                         Assert.assertEquals(f0, f.getValue(), JdkMath.abs(epsilon * f0));
456                         if (f.getOrder() > 0) {
457                             double dfdx = JdkMath.cos(x / a) / a;
458                             Assert.assertEquals(dfdx, f.getPartialDerivative(1, 0, 0), JdkMath.abs(epsilon * dfdx));
459                             double dfdy =  x * JdkMath.sin(y) * dfdx / a;
460                             Assert.assertEquals(dfdy, f.getPartialDerivative(0, 1, 0), JdkMath.abs(epsilon * dfdy));
461                             double cz = JdkMath.cos(z);
462                             double cz2 = cz * cz;
463                             double dfdz = -x * dfdx / (a * cz2);
464                             Assert.assertEquals(dfdz, f.getPartialDerivative(0, 0, 1), JdkMath.abs(epsilon * dfdz));
465                             if (f.getOrder() > 1) {
466                                 double df2dx2 = -(f0 / (a * a));
467                                 Assert.assertEquals(df2dx2, f.getPartialDerivative(2, 0, 0), JdkMath.abs(epsilon * df2dx2));
468                                 double df2dy2 = x * JdkMath.cos(y) * dfdx / a -
469                                                 x * x * JdkMath.sin(y) * JdkMath.sin(y) * f0 / (a * a * a * a) +
470                                                 2 * JdkMath.sin(y) * dfdy / a;
471                                 Assert.assertEquals(df2dy2, f.getPartialDerivative(0, 2, 0), JdkMath.abs(epsilon * df2dy2));
472                                 double c4 = cz2 * cz2;
473                                 double df2dz2 = x * (2 * a * (1 - a * cz * JdkMath.sin(z)) * dfdx - x * f0 / a ) / (a * a * a * c4);
474                                 Assert.assertEquals(df2dz2, f.getPartialDerivative(0, 0, 2), JdkMath.abs(epsilon * df2dz2));
475                                 double df2dxdy = dfdy / x  - x * JdkMath.sin(y) * f0 / (a * a * a);
476                                 Assert.assertEquals(df2dxdy, f.getPartialDerivative(1, 1, 0), JdkMath.abs(epsilon * df2dxdy));
477                             }
478                         }
479                     }
480                 }
481             }
482         }
483     }
484 
485     @Test
486     public void testSqrtDefinition() {
487         double[] epsilon = new double[] { 5.0e-16, 5.0e-16, 2.0e-15, 5.0e-14, 2.0e-12 };
488         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
489             for (double x = 0.1; x < 1.2; x += 0.001) {
490                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
491                 DerivativeStructure sqrt1 = dsX.pow(0.5);
492                 DerivativeStructure sqrt2 = dsX.sqrt();
493                 DerivativeStructure zero = sqrt1.subtract(sqrt2);
494                 for (int n = 0; n <= maxOrder; ++n) {
495                     Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
496                 }
497             }
498         }
499     }
500 
501     @Test
502     public void testRootNSingularity() {
503         for (int n = 2; n < 10; ++n) {
504             for (int maxOrder = 0; maxOrder < 12; ++maxOrder) {
505                 DerivativeStructure dsZero = new DerivativeStructure(1, maxOrder, 0, 0.0);
506                 DerivativeStructure rootN  = dsZero.rootN(n);
507                 Assert.assertEquals(0.0, rootN.getValue(), 1.0e-20);
508                 if (maxOrder > 0) {
509                     Assert.assertTrue(Double.isInfinite(rootN.getPartialDerivative(1)));
510                     Assert.assertTrue(rootN.getPartialDerivative(1) > 0);
511                     for (int order = 2; order <= maxOrder; ++order) {
512                         // the following checks shows a LIMITATION of the current implementation
513                         // we have no way to tell dsZero is a pure linear variable x = 0
514                         // we only say: "dsZero is a structure with value = 0.0,
515                         // first derivative = 1.0, second and higher derivatives = 0.0".
516                         // Function composition rule for second derivatives is:
517                         // d2[f(g(x))]/dx2 = f''(g(x)) * [g'(x)]^2 + f'(g(x)) * g''(x)
518                         // when function f is the nth root and x = 0 we have:
519                         // f(0) = 0, f'(0) = +infinity, f''(0) = -infinity (and higher
520                         // derivatives keep switching between +infinity and -infinity)
521                         // so given that in our case dsZero represents g, we have g(x) = 0,
522                         // g'(x) = 1 and g''(x) = 0
523                         // applying the composition rules gives:
524                         // d2[f(g(x))]/dx2 = f''(g(x)) * [g'(x)]^2 + f'(g(x)) * g''(x)
525                         //                 = -infinity * 1^2 + +infinity * 0
526                         //                 = -infinity + NaN
527                         //                 = NaN
528                         // if we knew dsZero is really the x variable and not the identity
529                         // function applied to x, we would not have computed f'(g(x)) * g''(x)
530                         // and we would have found that the result was -infinity and not NaN
531                         Assert.assertTrue(Double.isNaN(rootN.getPartialDerivative(order)));
532                     }
533                 }
534 
535                 // the following shows that the limitation explained above is NOT a bug...
536                 // if we set up the higher order derivatives for g appropriately, we do
537                 // compute the higher order derivatives of the composition correctly
538                 double[] gDerivatives = new double[ 1 + maxOrder];
539                 gDerivatives[0] = 0.0;
540                 for (int k = 1; k <= maxOrder; ++k) {
541                     gDerivatives[k] = JdkMath.pow(-1.0, k + 1);
542                 }
543                 DerivativeStructure correctRoot = new DerivativeStructure(1, maxOrder, gDerivatives).rootN(n);
544                 Assert.assertEquals(0.0, correctRoot.getValue(), 1.0e-20);
545                 if (maxOrder > 0) {
546                     Assert.assertTrue(Double.isInfinite(correctRoot.getPartialDerivative(1)));
547                     Assert.assertTrue(correctRoot.getPartialDerivative(1) > 0);
548                     for (int order = 2; order <= maxOrder; ++order) {
549                         Assert.assertTrue(Double.isInfinite(correctRoot.getPartialDerivative(order)));
550                         if ((order & 1) == 0) {
551                             Assert.assertTrue(correctRoot.getPartialDerivative(order) < 0);
552                         } else {
553                             Assert.assertTrue(correctRoot.getPartialDerivative(order) > 0);
554                         }
555                     }
556                 }
557             }
558         }
559     }
560 
561     @Test
562     public void testSqrtPow2() {
563         double[] epsilon = new double[] { 1.0e-16, 3.0e-16, 2.0e-15, 6.0e-14, 6.0e-12 };
564         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
565             for (double x = 0.1; x < 1.2; x += 0.001) {
566                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
567                 DerivativeStructure rebuiltX = dsX.multiply(dsX).sqrt();
568                 DerivativeStructure zero = rebuiltX.subtract(dsX);
569                 for (int n = 0; n <= maxOrder; ++n) {
570                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
571                 }
572             }
573         }
574     }
575 
576     @Test
577     public void testCbrtDefinition() {
578         double[] epsilon = new double[] { 4.0e-16, 9.0e-16, 6.0e-15, 2.0e-13, 4.0e-12 };
579         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
580             for (double x = 0.1; x < 1.2; x += 0.001) {
581                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
582                 DerivativeStructure cbrt1 = dsX.pow(1.0 / 3.0);
583                 DerivativeStructure cbrt2 = dsX.cbrt();
584                 DerivativeStructure zero = cbrt1.subtract(cbrt2);
585                 for (int n = 0; n <= maxOrder; ++n) {
586                     Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
587                 }
588             }
589         }
590     }
591 
592     @Test
593     public void testCbrtPow3() {
594         double[] epsilon = new double[] { 1.0e-16, 5.0e-16, 8.0e-15, 3.0e-13, 4.0e-11 };
595         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
596             for (double x = 0.1; x < 1.2; x += 0.001) {
597                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
598                 DerivativeStructure rebuiltX = dsX.multiply(dsX.multiply(dsX)).cbrt();
599                 DerivativeStructure zero = rebuiltX.subtract(dsX);
600                 for (int n = 0; n <= maxOrder; ++n) {
601                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
602                 }
603             }
604         }
605     }
606 
607     @Test
608     public void testPowReciprocalPow() {
609         double[] epsilon = new double[] { 2.0e-15, 2.0e-14, 3.0e-13, 8.0e-12, 3.0e-10 };
610         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
611             for (double x = 0.1; x < 1.2; x += 0.01) {
612                 DerivativeStructure dsX = new DerivativeStructure(2, maxOrder, 0, x);
613                 for (double y = 0.1; y < 1.2; y += 0.01) {
614                     DerivativeStructure dsY = new DerivativeStructure(2, maxOrder, 1, y);
615                     DerivativeStructure rebuiltX = dsX.pow(dsY).pow(dsY.reciprocal());
616                     DerivativeStructure zero = rebuiltX.subtract(dsX);
617                     for (int n = 0; n <= maxOrder; ++n) {
618                         for (int m = 0; m <= maxOrder; ++m) {
619                             if (n + m <= maxOrder) {
620                                 Assert.assertEquals(0.0, zero.getPartialDerivative(n, m), epsilon[n + m]);
621                             }
622                         }
623                     }
624                 }
625             }
626         }
627     }
628 
629     @Test
630     public void testHypotDefinition() {
631         double epsilon = 1.0e-20;
632         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
633             for (double x = -1.7; x < 2; x += 0.2) {
634                 DerivativeStructure dsX = new DerivativeStructure(2, maxOrder, 0, x);
635                 for (double y = -1.7; y < 2; y += 0.2) {
636                     DerivativeStructure dsY = new DerivativeStructure(2, maxOrder, 1, y);
637                     DerivativeStructure hypot = DerivativeStructure.hypot(dsY, dsX);
638                     DerivativeStructure ref = dsX.multiply(dsX).add(dsY.multiply(dsY)).sqrt();
639                     DerivativeStructure zero = hypot.subtract(ref);
640                     for (int n = 0; n <= maxOrder; ++n) {
641                         for (int m = 0; m <= maxOrder; ++m) {
642                             if (n + m <= maxOrder) {
643                                 Assert.assertEquals(0, zero.getPartialDerivative(n, m), epsilon);
644                             }
645                         }
646                     }
647                 }
648             }
649         }
650     }
651 
652     @Test
653     public void testHypotNoOverflow() {
654 
655         DerivativeStructure dsX = new DerivativeStructure(2, 5, 0, +3.0e250);
656         DerivativeStructure dsY = new DerivativeStructure(2, 5, 1, -4.0e250);
657         DerivativeStructure hypot = DerivativeStructure.hypot(dsX, dsY);
658         Assert.assertEquals(5.0e250, hypot.getValue(), 1.0e235);
659         Assert.assertEquals(dsX.getValue() / hypot.getValue(), hypot.getPartialDerivative(1, 0), 1.0e-10);
660         Assert.assertEquals(dsY.getValue() / hypot.getValue(), hypot.getPartialDerivative(0, 1), 1.0e-10);
661 
662         DerivativeStructure sqrt  = dsX.multiply(dsX).add(dsY.multiply(dsY)).sqrt();
663         Assert.assertTrue(Double.isInfinite(sqrt.getValue()));
664     }
665 
666     @Test
667     public void testHypotNeglectible() {
668 
669         DerivativeStructure dsSmall = new DerivativeStructure(2, 5, 0, +3.0e-10);
670         DerivativeStructure dsLarge = new DerivativeStructure(2, 5, 1, -4.0e25);
671 
672         Assert.assertEquals(dsLarge.abs().getValue(),
673                             DerivativeStructure.hypot(dsSmall, dsLarge).getValue(),
674                             1.0e-10);
675         Assert.assertEquals(0,
676                             DerivativeStructure.hypot(dsSmall, dsLarge).getPartialDerivative(1, 0),
677                             1.0e-10);
678         Assert.assertEquals(-1,
679                             DerivativeStructure.hypot(dsSmall, dsLarge).getPartialDerivative(0, 1),
680                             1.0e-10);
681 
682         Assert.assertEquals(dsLarge.abs().getValue(),
683                             DerivativeStructure.hypot(dsLarge, dsSmall).getValue(),
684                             1.0e-10);
685         Assert.assertEquals(0,
686                             DerivativeStructure.hypot(dsLarge, dsSmall).getPartialDerivative(1, 0),
687                             1.0e-10);
688         Assert.assertEquals(-1,
689                             DerivativeStructure.hypot(dsLarge, dsSmall).getPartialDerivative(0, 1),
690                             1.0e-10);
691     }
692 
693     @Test
694     public void testHypotSpecial() {
695         Assert.assertTrue(Double.isNaN(DerivativeStructure.hypot(new DerivativeStructure(2, 5, 0, Double.NaN),
696                                                                  new DerivativeStructure(2, 5, 0, +3.0e250)).getValue()));
697         Assert.assertTrue(Double.isNaN(DerivativeStructure.hypot(new DerivativeStructure(2, 5, 0, +3.0e250),
698                                                                  new DerivativeStructure(2, 5, 0, Double.NaN)).getValue()));
699         Assert.assertTrue(Double.isInfinite(DerivativeStructure.hypot(new DerivativeStructure(2, 5, 0, Double.POSITIVE_INFINITY),
700                                                                       new DerivativeStructure(2, 5, 0, +3.0e250)).getValue()));
701         Assert.assertTrue(Double.isInfinite(DerivativeStructure.hypot(new DerivativeStructure(2, 5, 0, +3.0e250),
702                                                                       new DerivativeStructure(2, 5, 0, Double.POSITIVE_INFINITY)).getValue()));
703     }
704 
705     @Test
706     public void testPrimitiveRemainder() {
707         double epsilon = 1.0e-15;
708         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
709             for (double x = -1.7; x < 2; x += 0.2) {
710                 DerivativeStructure dsX = new DerivativeStructure(2, maxOrder, 0, x);
711                 for (double y = -1.7; y < 2; y += 0.2) {
712                     DerivativeStructure remainder = dsX.remainder(y);
713                     DerivativeStructure ref = dsX.subtract(x - JdkMath.IEEEremainder(x, y));
714                     DerivativeStructure zero = remainder.subtract(ref);
715                     for (int n = 0; n <= maxOrder; ++n) {
716                         for (int m = 0; m <= maxOrder; ++m) {
717                             if (n + m <= maxOrder) {
718                                 Assert.assertEquals(0, zero.getPartialDerivative(n, m), epsilon);
719                             }
720                         }
721                     }
722                 }
723             }
724         }
725     }
726 
727     @Test
728     public void testRemainder() {
729         double epsilon = 2.0e-15;
730         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
731             for (double x = -1.7; x < 2; x += 0.2) {
732                 DerivativeStructure dsX = new DerivativeStructure(2, maxOrder, 0, x);
733                 for (double y = -1.7; y < 2; y += 0.2) {
734                     DerivativeStructure dsY = new DerivativeStructure(2, maxOrder, 1, y);
735                     DerivativeStructure remainder = dsX.remainder(dsY);
736                     DerivativeStructure ref = dsX.subtract(dsY.multiply((x - JdkMath.IEEEremainder(x, y)) / y));
737                     DerivativeStructure zero = remainder.subtract(ref);
738                     for (int n = 0; n <= maxOrder; ++n) {
739                         for (int m = 0; m <= maxOrder; ++m) {
740                             if (n + m <= maxOrder) {
741                                 Assert.assertEquals(0, zero.getPartialDerivative(n, m), epsilon);
742                             }
743                         }
744                     }
745                 }
746             }
747         }
748     }
749 
750     @Override
751     @Test
752     public void testExp() {
753         double[] epsilon = new double[] { 1.0e-16, 1.0e-16, 1.0e-16, 1.0e-16, 1.0e-16 };
754         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
755             for (double x = 0.1; x < 1.2; x += 0.001) {
756                 double refExp = JdkMath.exp(x);
757                 DerivativeStructure exp = new DerivativeStructure(1, maxOrder, 0, x).exp();
758                 for (int n = 0; n <= maxOrder; ++n) {
759                     Assert.assertEquals(refExp, exp.getPartialDerivative(n), epsilon[n]);
760                 }
761             }
762         }
763     }
764 
765     @Test
766     public void testExpm1Definition() {
767         double epsilon = 3.0e-16;
768         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
769             for (double x = 0.1; x < 1.2; x += 0.001) {
770                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
771                 DerivativeStructure expm11 = dsX.expm1();
772                 DerivativeStructure expm12 = dsX.exp().subtract(dsX.getField().getOne());
773                 DerivativeStructure zero = expm11.subtract(expm12);
774                 for (int n = 0; n <= maxOrder; ++n) {
775                     Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon);
776                 }
777             }
778         }
779     }
780 
781     @Override
782     @Test
783     public void testLog() {
784         double[] epsilon = new double[] { 1.0e-16, 1.0e-16, 3.0e-14, 7.0e-13, 3.0e-11 };
785         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
786             for (double x = 0.1; x < 1.2; x += 0.001) {
787                 DerivativeStructure log = new DerivativeStructure(1, maxOrder, 0, x).log();
788                 Assert.assertEquals(JdkMath.log(x), log.getValue(), epsilon[0]);
789                 for (int n = 1; n <= maxOrder; ++n) {
790                     double refDer = -Factorial.value(n - 1) / JdkMath.pow(-x, n);
791                     Assert.assertEquals(refDer, log.getPartialDerivative(n), epsilon[n]);
792                 }
793             }
794         }
795     }
796 
797     @Test
798     public void testLog1pDefinition() {
799         double epsilon = 3.0e-16;
800         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
801             for (double x = 0.1; x < 1.2; x += 0.001) {
802                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
803                 DerivativeStructure log1p1 = dsX.log1p();
804                 DerivativeStructure log1p2 = dsX.add(dsX.getField().getOne()).log();
805                 DerivativeStructure zero = log1p1.subtract(log1p2);
806                 for (int n = 0; n <= maxOrder; ++n) {
807                     Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon);
808                 }
809             }
810         }
811     }
812 
813     @Test
814     public void testLog10Definition() {
815         double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 8.0e-15, 3.0e-13, 8.0e-12 };
816         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
817             for (double x = 0.1; x < 1.2; x += 0.001) {
818                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
819                 DerivativeStructure log101 = dsX.log10();
820                 DerivativeStructure log102 = dsX.log().divide(JdkMath.log(10.0));
821                 DerivativeStructure zero = log101.subtract(log102);
822                 for (int n = 0; n <= maxOrder; ++n) {
823                     Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
824                 }
825             }
826         }
827     }
828 
829     @Test
830     public void testLogExp() {
831         double[] epsilon = new double[] { 2.0e-16, 2.0e-16, 3.0e-16, 2.0e-15, 6.0e-15 };
832         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
833             for (double x = 0.1; x < 1.2; x += 0.001) {
834                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
835                 DerivativeStructure rebuiltX = dsX.exp().log();
836                 DerivativeStructure zero = rebuiltX.subtract(dsX);
837                 for (int n = 0; n <= maxOrder; ++n) {
838                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
839                 }
840             }
841         }
842     }
843 
844     @Test
845     public void testLog1pExpm1() {
846         double[] epsilon = new double[] { 6.0e-17, 3.0e-16, 5.0e-16, 9.0e-16, 6.0e-15 };
847         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
848             for (double x = 0.1; x < 1.2; x += 0.001) {
849                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
850                 DerivativeStructure rebuiltX = dsX.expm1().log1p();
851                 DerivativeStructure zero = rebuiltX.subtract(dsX);
852                 for (int n = 0; n <= maxOrder; ++n) {
853                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
854                 }
855             }
856         }
857     }
858 
859     @Test
860     public void testLog10Power() {
861         double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 9.0e-16, 6.0e-15, 6.0e-14 };
862         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
863             for (double x = 0.1; x < 1.2; x += 0.001) {
864                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
865                 DerivativeStructure rebuiltX = new DerivativeStructure(1, maxOrder, 10.0).pow(dsX).log10();
866                 DerivativeStructure zero = rebuiltX.subtract(dsX);
867                 for (int n = 0; n <= maxOrder; ++n) {
868                     Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
869                 }
870             }
871         }
872     }
873 
874     @Test
875     public void testSinCos() {
876         double epsilon = 5.0e-16;
877         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
878             for (double x = 0.1; x < 1.2; x += 0.001) {
879                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
880                 DerivativeStructure sin = dsX.sin();
881                 DerivativeStructure cos = dsX.cos();
882                 double s = JdkMath.sin(x);
883                 double c = JdkMath.cos(x);
884                 for (int n = 0; n <= maxOrder; ++n) {
885                     switch (n % 4) {
886                     case 0 :
887                         Assert.assertEquals( s, sin.getPartialDerivative(n), epsilon);
888                         Assert.assertEquals( c, cos.getPartialDerivative(n), epsilon);
889                         break;
890                     case 1 :
891                         Assert.assertEquals( c, sin.getPartialDerivative(n), epsilon);
892                         Assert.assertEquals(-s, cos.getPartialDerivative(n), epsilon);
893                         break;
894                     case 2 :
895                         Assert.assertEquals(-s, sin.getPartialDerivative(n), epsilon);
896                         Assert.assertEquals(-c, cos.getPartialDerivative(n), epsilon);
897                         break;
898                     default :
899                         Assert.assertEquals(-c, sin.getPartialDerivative(n), epsilon);
900                         Assert.assertEquals( s, cos.getPartialDerivative(n), epsilon);
901                         break;
902                     }
903                 }
904             }
905         }
906     }
907 
908     @Test
909     public void testSinAsin() {
910         double[] epsilon = new double[] { 3.0e-16, 5.0e-16, 3.0e-15, 2.0e-14, 4.0e-13 };
911         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
912             for (double x = 0.1; x < 1.2; x += 0.001) {
913                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
914                 DerivativeStructure rebuiltX = dsX.sin().asin();
915                 DerivativeStructure zero = rebuiltX.subtract(dsX);
916                 for (int n = 0; n <= maxOrder; ++n) {
917                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
918                 }
919             }
920         }
921     }
922 
923     @Test
924     public void testCosAcos() {
925         double[] epsilon = new double[] { 6.0e-16, 6.0e-15, 2.0e-13, 4.0e-12, 2.0e-10 };
926         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
927             for (double x = 0.1; x < 1.2; x += 0.001) {
928                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
929                 DerivativeStructure rebuiltX = dsX.cos().acos();
930                 DerivativeStructure zero = rebuiltX.subtract(dsX);
931                 for (int n = 0; n <= maxOrder; ++n) {
932                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
933                 }
934             }
935         }
936     }
937 
938     @Test
939     public void testTanAtan() {
940         double[] epsilon = new double[] { 6.0e-17, 2.0e-16, 2.0e-15, 4.0e-14, 2.0e-12 };
941         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
942             for (double x = 0.1; x < 1.2; x += 0.001) {
943                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
944                 DerivativeStructure rebuiltX = dsX.tan().atan();
945                 DerivativeStructure zero = rebuiltX.subtract(dsX);
946                 for (int n = 0; n <= maxOrder; ++n) {
947                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
948                 }
949             }
950         }
951     }
952 
953     @Test
954     public void testTangentDefinition() {
955         double[] epsilon = new double[] { 5.0e-16, 2.0e-15, 3.0e-14, 5.0e-13, 2.0e-11 };
956         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
957             for (double x = 0.1; x < 1.2; x += 0.001) {
958                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
959                 DerivativeStructure tan1 = dsX.sin().divide(dsX.cos());
960                 DerivativeStructure tan2 = dsX.tan();
961                 DerivativeStructure zero = tan1.subtract(tan2);
962                 for (int n = 0; n <= maxOrder; ++n) {
963                     Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
964                 }
965             }
966         }
967     }
968 
969     @Override
970     @Test
971     public void testAtan2() {
972         double[] epsilon = new double[] { 5.0e-16, 3.0e-15, 2.2e-14, 1.0e-12, 8.0e-11 };
973         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
974             for (double x = -1.7; x < 2; x += 0.2) {
975                 DerivativeStructure dsX = new DerivativeStructure(2, maxOrder, 0, x);
976                 for (double y = -1.7; y < 2; y += 0.2) {
977                     DerivativeStructure dsY = new DerivativeStructure(2, maxOrder, 1, y);
978                     DerivativeStructure atan2 = DerivativeStructure.atan2(dsY, dsX);
979                     DerivativeStructure ref = dsY.divide(dsX).atan();
980                     if (x < 0) {
981                         ref = (y < 0) ? ref.subtract(JdkMath.PI) : ref.add(JdkMath.PI);
982                     }
983                     DerivativeStructure zero = atan2.subtract(ref);
984                     for (int n = 0; n <= maxOrder; ++n) {
985                         for (int m = 0; m <= maxOrder; ++m) {
986                             if (n + m <= maxOrder) {
987                                 Assert.assertEquals(0, zero.getPartialDerivative(n, m), epsilon[n + m]);
988                             }
989                         }
990                     }
991                 }
992             }
993         }
994     }
995 
996     @Test
997     public void testAtan2SpecialCases() {
998 
999         DerivativeStructure pp =
1000                 DerivativeStructure.atan2(new DerivativeStructure(2, 2, 1, +0.0),
1001                                           new DerivativeStructure(2, 2, 1, +0.0));
1002         Assert.assertEquals(0, pp.getValue(), 1.0e-15);
1003         Assert.assertEquals(+1, JdkMath.copySign(1, pp.getValue()), 1.0e-15);
1004 
1005         DerivativeStructure pn =
1006                 DerivativeStructure.atan2(new DerivativeStructure(2, 2, 1, +0.0),
1007                                           new DerivativeStructure(2, 2, 1, -0.0));
1008         Assert.assertEquals(JdkMath.PI, pn.getValue(), 1.0e-15);
1009 
1010         DerivativeStructure np =
1011                 DerivativeStructure.atan2(new DerivativeStructure(2, 2, 1, -0.0),
1012                                           new DerivativeStructure(2, 2, 1, +0.0));
1013         Assert.assertEquals(0, np.getValue(), 1.0e-15);
1014         Assert.assertEquals(-1, JdkMath.copySign(1, np.getValue()), 1.0e-15);
1015 
1016         DerivativeStructure nn =
1017                 DerivativeStructure.atan2(new DerivativeStructure(2, 2, 1, -0.0),
1018                                           new DerivativeStructure(2, 2, 1, -0.0));
1019         Assert.assertEquals(-JdkMath.PI, nn.getValue(), 1.0e-15);
1020     }
1021 
1022     @Test
1023     public void testSinhDefinition() {
1024         double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 5.0e-16, 2.0e-15, 6.0e-15 };
1025         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1026             for (double x = 0.1; x < 1.2; x += 0.001) {
1027                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
1028                 DerivativeStructure sinh1 = dsX.exp().subtract(dsX.exp().reciprocal()).multiply(0.5);
1029                 DerivativeStructure sinh2 = dsX.sinh();
1030                 DerivativeStructure zero = sinh1.subtract(sinh2);
1031                 for (int n = 0; n <= maxOrder; ++n) {
1032                     Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
1033                 }
1034             }
1035         }
1036     }
1037 
1038     @Test
1039     public void testCoshDefinition() {
1040         double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 5.0e-16, 2.0e-15, 6.0e-15 };
1041         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1042             for (double x = 0.1; x < 1.2; x += 0.001) {
1043                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
1044                 DerivativeStructure cosh1 = dsX.exp().add(dsX.exp().reciprocal()).multiply(0.5);
1045                 DerivativeStructure cosh2 = dsX.cosh();
1046                 DerivativeStructure zero = cosh1.subtract(cosh2);
1047                 for (int n = 0; n <= maxOrder; ++n) {
1048                     Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
1049                 }
1050             }
1051         }
1052     }
1053 
1054     @Test
1055     public void testTanhDefinition() {
1056         double[] epsilon = new double[] { 3.0e-16, 5.0e-16, 7.0e-16, 3.0e-15, 2.0e-14 };
1057         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1058             for (double x = 0.1; x < 1.2; x += 0.001) {
1059                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
1060                 DerivativeStructure tanh1 = dsX.exp().subtract(dsX.exp().reciprocal()).divide(dsX.exp().add(dsX.exp().reciprocal()));
1061                 DerivativeStructure tanh2 = dsX.tanh();
1062                 DerivativeStructure zero = tanh1.subtract(tanh2);
1063                 for (int n = 0; n <= maxOrder; ++n) {
1064                     Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
1065                 }
1066             }
1067         }
1068     }
1069 
1070     @Test
1071     public void testSinhAsinh() {
1072         double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 4.0e-16, 7.0e-16, 3.0e-15, 8.0e-15 };
1073         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1074             for (double x = 0.1; x < 1.2; x += 0.001) {
1075                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
1076                 DerivativeStructure rebuiltX = dsX.sinh().asinh();
1077                 DerivativeStructure zero = rebuiltX.subtract(dsX);
1078                 for (int n = 0; n <= maxOrder; ++n) {
1079                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
1080                 }
1081             }
1082         }
1083     }
1084 
1085     @Test
1086     public void testCoshAcosh() {
1087         double[] epsilon = new double[] { 2.0e-15, 1.0e-14, 2.0e-13, 6.0e-12, 3.0e-10, 2.0e-8 };
1088         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1089             for (double x = 0.1; x < 1.2; x += 0.001) {
1090                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
1091                 DerivativeStructure rebuiltX = dsX.cosh().acosh();
1092                 DerivativeStructure zero = rebuiltX.subtract(dsX);
1093                 for (int n = 0; n <= maxOrder; ++n) {
1094                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
1095                 }
1096             }
1097         }
1098     }
1099 
1100     @Test
1101     public void testTanhAtanh() {
1102         double[] epsilon = new double[] { 3.0e-16, 2.0e-16, 7.0e-16, 4.0e-15, 3.0e-14, 4.0e-13 };
1103         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1104             for (double x = 0.1; x < 1.2; x += 0.001) {
1105                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
1106                 DerivativeStructure rebuiltX = dsX.tanh().atanh();
1107                 DerivativeStructure zero = rebuiltX.subtract(dsX);
1108                 for (int n = 0; n <= maxOrder; ++n) {
1109                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
1110                 }
1111             }
1112         }
1113     }
1114 
1115     @Test
1116     public void testCompositionOneVariableY() {
1117         double epsilon = 1.0e-13;
1118         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1119             for (double x = 0.1; x < 1.2; x += 0.1) {
1120                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, x);
1121                 for (double y = 0.1; y < 1.2; y += 0.1) {
1122                     DerivativeStructure dsY = new DerivativeStructure(1, maxOrder, 0, y);
1123                     DerivativeStructure f = dsX.divide(dsY).sqrt();
1124                     double f0 = JdkMath.sqrt(x / y);
1125                     Assert.assertEquals(f0, f.getValue(), JdkMath.abs(epsilon * f0));
1126                     if (f.getOrder() > 0) {
1127                         double f1 = -x / (2 * y * y * f0);
1128                         Assert.assertEquals(f1, f.getPartialDerivative(1), JdkMath.abs(epsilon * f1));
1129                         if (f.getOrder() > 1) {
1130                             double f2 = (f0 - x / (4 * y * f0)) / (y * y);
1131                             Assert.assertEquals(f2, f.getPartialDerivative(2), JdkMath.abs(epsilon * f2));
1132                             if (f.getOrder() > 2) {
1133                                 double f3 = (x / (8 * y * f0) - 2 * f0) / (y * y * y);
1134                                 Assert.assertEquals(f3, f.getPartialDerivative(3), JdkMath.abs(epsilon * f3));
1135                             }
1136                         }
1137                     }
1138                 }
1139             }
1140         }
1141     }
1142 
1143     @Test
1144     public void testTaylorPolynomial() {
1145         for (double x = 0; x < 1.2; x += 0.1) {
1146             DerivativeStructure dsX = new DerivativeStructure(3, 4, 0, x);
1147             for (double y = 0; y < 1.2; y += 0.2) {
1148                 DerivativeStructure dsY = new DerivativeStructure(3, 4, 1, y);
1149                 for (double z = 0; z < 1.2; z += 0.2) {
1150                     DerivativeStructure dsZ = new DerivativeStructure(3, 4, 2, z);
1151                     DerivativeStructure f = dsX.multiply(dsY).add(dsZ).multiply(dsX).multiply(dsY);
1152                     for (double dx = -0.2; dx < 0.2; dx += 0.2) {
1153                         for (double dy = -0.2; dy < 0.2; dy += 0.1) {
1154                             for (double dz = -0.2; dz < 0.2; dz += 0.1) {
1155                                 double ref = (x + dx) * (y + dy) * ((x + dx) * (y + dy) + (z + dz));
1156                                 Assert.assertEquals(ref, f.taylor(dx, dy, dz), 2.0e-15);
1157                             }
1158                         }
1159                     }
1160                 }
1161             }
1162         }
1163     }
1164 
1165     @Test
1166     public void testTaylorAtan2() {
1167         double[] expected = new double[] { 0.214, 0.0241, 0.00422, 6.48e-4, 8.04e-5 };
1168         double x0 =  0.1;
1169         double y0 = -0.3;
1170         for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
1171             DerivativeStructure dsX   = new DerivativeStructure(2, maxOrder, 0, x0);
1172             DerivativeStructure dsY   = new DerivativeStructure(2, maxOrder, 1, y0);
1173             DerivativeStructure atan2 = DerivativeStructure.atan2(dsY, dsX);
1174             double maxError = 0;
1175             for (double dx = -0.05; dx < 0.05; dx += 0.001) {
1176                 for (double dy = -0.05; dy < 0.05; dy += 0.001) {
1177                     double ref = JdkMath.atan2(y0 + dy, x0 + dx);
1178                     maxError = JdkMath.max(maxError, JdkMath.abs(ref - atan2.taylor(dx, dy)));
1179                 }
1180             }
1181             Assert.assertEquals(0.0, expected[maxOrder] - maxError, 0.01 * expected[maxOrder]);
1182         }
1183     }
1184 
1185     @Override
1186     @Test
1187     public void testAbs() {
1188 
1189         DerivativeStructure minusOne = new DerivativeStructure(1, 1, 0, -1.0);
1190         Assert.assertEquals(+1.0, minusOne.abs().getPartialDerivative(0), 1.0e-15);
1191         Assert.assertEquals(-1.0, minusOne.abs().getPartialDerivative(1), 1.0e-15);
1192 
1193         DerivativeStructure plusOne = new DerivativeStructure(1, 1, 0, +1.0);
1194         Assert.assertEquals(+1.0, plusOne.abs().getPartialDerivative(0), 1.0e-15);
1195         Assert.assertEquals(+1.0, plusOne.abs().getPartialDerivative(1), 1.0e-15);
1196 
1197         DerivativeStructure minusZero = new DerivativeStructure(1, 1, 0, -0.0);
1198         Assert.assertEquals(+0.0, minusZero.abs().getPartialDerivative(0), 1.0e-15);
1199         Assert.assertEquals(-1.0, minusZero.abs().getPartialDerivative(1), 1.0e-15);
1200 
1201         DerivativeStructure plusZero = new DerivativeStructure(1, 1, 0, +0.0);
1202         Assert.assertEquals(+0.0, plusZero.abs().getPartialDerivative(0), 1.0e-15);
1203         Assert.assertEquals(+1.0, plusZero.abs().getPartialDerivative(1), 1.0e-15);
1204     }
1205 
1206     @Override
1207     @Test
1208     public void testSignum() {
1209 
1210         DerivativeStructure minusOne = new DerivativeStructure(1, 1, 0, -1.0);
1211         Assert.assertEquals(-1.0, minusOne.signum().getPartialDerivative(0), 1.0e-15);
1212         Assert.assertEquals( 0.0, minusOne.signum().getPartialDerivative(1), 1.0e-15);
1213 
1214         DerivativeStructure plusOne = new DerivativeStructure(1, 1, 0, +1.0);
1215         Assert.assertEquals(+1.0, plusOne.signum().getPartialDerivative(0), 1.0e-15);
1216         Assert.assertEquals( 0.0, plusOne.signum().getPartialDerivative(1), 1.0e-15);
1217 
1218         DerivativeStructure minusZero = new DerivativeStructure(1, 1, 0, -0.0);
1219         Assert.assertEquals(-0.0, minusZero.signum().getPartialDerivative(0), 1.0e-15);
1220         Assert.assertTrue(Double.doubleToLongBits(minusZero.signum().getValue()) < 0);
1221         Assert.assertEquals( 0.0, minusZero.signum().getPartialDerivative(1), 1.0e-15);
1222 
1223         DerivativeStructure plusZero = new DerivativeStructure(1, 1, 0, +0.0);
1224         Assert.assertEquals(+0.0, plusZero.signum().getPartialDerivative(0), 1.0e-15);
1225         Assert.assertEquals(0, Double.doubleToLongBits(plusZero.signum().getValue()));
1226         Assert.assertEquals( 0.0, plusZero.signum().getPartialDerivative(1), 1.0e-15);
1227     }
1228 
1229     @Test
1230     public void testCeilFloorRintLong() {
1231 
1232         DerivativeStructure x = new DerivativeStructure(1, 1, 0, -1.5);
1233         Assert.assertEquals(-1.5, x.getPartialDerivative(0), 1.0e-15);
1234         Assert.assertEquals(+1.0, x.getPartialDerivative(1), 1.0e-15);
1235         Assert.assertEquals(-1.0, x.ceil().getPartialDerivative(0), 1.0e-15);
1236         Assert.assertEquals(+0.0, x.ceil().getPartialDerivative(1), 1.0e-15);
1237         Assert.assertEquals(-2.0, x.floor().getPartialDerivative(0), 1.0e-15);
1238         Assert.assertEquals(+0.0, x.floor().getPartialDerivative(1), 1.0e-15);
1239         Assert.assertEquals(-2.0, x.rint().getPartialDerivative(0), 1.0e-15);
1240         Assert.assertEquals(+0.0, x.rint().getPartialDerivative(1), 1.0e-15);
1241         Assert.assertEquals(-2.0, x.subtract(x.getField().getOne()).rint().getPartialDerivative(0), 1.0e-15);
1242         Assert.assertEquals(-1L, x.round());
1243     }
1244 
1245     @Test
1246     public void testCopySign() {
1247 
1248         DerivativeStructure minusOne = new DerivativeStructure(1, 1, 0, -1.0);
1249         Assert.assertEquals(+1.0, minusOne.copySign(+1.0).getPartialDerivative(0), 1.0e-15);
1250         Assert.assertEquals(-1.0, minusOne.copySign(+1.0).getPartialDerivative(1), 1.0e-15);
1251         Assert.assertEquals(-1.0, minusOne.copySign(-1.0).getPartialDerivative(0), 1.0e-15);
1252         Assert.assertEquals(+1.0, minusOne.copySign(-1.0).getPartialDerivative(1), 1.0e-15);
1253         Assert.assertEquals(+1.0, minusOne.copySign(+0.0).getPartialDerivative(0), 1.0e-15);
1254         Assert.assertEquals(-1.0, minusOne.copySign(+0.0).getPartialDerivative(1), 1.0e-15);
1255         Assert.assertEquals(-1.0, minusOne.copySign(-0.0).getPartialDerivative(0), 1.0e-15);
1256         Assert.assertEquals(+1.0, minusOne.copySign(-0.0).getPartialDerivative(1), 1.0e-15);
1257         Assert.assertEquals(+1.0, minusOne.copySign(Double.NaN).getPartialDerivative(0), 1.0e-15);
1258         Assert.assertEquals(-1.0, minusOne.copySign(Double.NaN).getPartialDerivative(1), 1.0e-15);
1259 
1260         DerivativeStructure plusOne = new DerivativeStructure(1, 1, 0, +1.0);
1261         Assert.assertEquals(+1.0, plusOne.copySign(+1.0).getPartialDerivative(0), 1.0e-15);
1262         Assert.assertEquals(+1.0, plusOne.copySign(+1.0).getPartialDerivative(1), 1.0e-15);
1263         Assert.assertEquals(-1.0, plusOne.copySign(-1.0).getPartialDerivative(0), 1.0e-15);
1264         Assert.assertEquals(-1.0, plusOne.copySign(-1.0).getPartialDerivative(1), 1.0e-15);
1265         Assert.assertEquals(+1.0, plusOne.copySign(+0.0).getPartialDerivative(0), 1.0e-15);
1266         Assert.assertEquals(+1.0, plusOne.copySign(+0.0).getPartialDerivative(1), 1.0e-15);
1267         Assert.assertEquals(-1.0, plusOne.copySign(-0.0).getPartialDerivative(0), 1.0e-15);
1268         Assert.assertEquals(-1.0, plusOne.copySign(-0.0).getPartialDerivative(1), 1.0e-15);
1269         Assert.assertEquals(+1.0, plusOne.copySign(Double.NaN).getPartialDerivative(0), 1.0e-15);
1270         Assert.assertEquals(+1.0, plusOne.copySign(Double.NaN).getPartialDerivative(1), 1.0e-15);
1271     }
1272 
1273     @Test
1274     public void testToDegreesDefinition() {
1275         double epsilon = 3.0e-16;
1276         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1277             for (double x = 0.1; x < 1.2; x += 0.001) {
1278                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
1279                 Assert.assertEquals(JdkMath.toDegrees(x), dsX.toDegrees().getValue(), epsilon);
1280                 for (int n = 1; n <= maxOrder; ++n) {
1281                     if (n == 1) {
1282                         Assert.assertEquals(180 / JdkMath.PI, dsX.toDegrees().getPartialDerivative(1), epsilon);
1283                     } else {
1284                         Assert.assertEquals(0.0, dsX.toDegrees().getPartialDerivative(n), epsilon);
1285                     }
1286                 }
1287             }
1288         }
1289     }
1290 
1291     @Test
1292     public void testToRadiansDefinition() {
1293         double epsilon = 3.0e-16;
1294         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1295             for (double x = 0.1; x < 1.2; x += 0.001) {
1296                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
1297                 Assert.assertEquals(JdkMath.toRadians(x), dsX.toRadians().getValue(), epsilon);
1298                 for (int n = 1; n <= maxOrder; ++n) {
1299                     if (n == 1) {
1300                         Assert.assertEquals(JdkMath.PI / 180, dsX.toRadians().getPartialDerivative(1), epsilon);
1301                     } else {
1302                         Assert.assertEquals(0.0, dsX.toRadians().getPartialDerivative(n), epsilon);
1303                     }
1304                 }
1305             }
1306         }
1307     }
1308 
1309     @Test
1310     public void testDegRad() {
1311         double epsilon = 3.0e-16;
1312         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1313             for (double x = 0.1; x < 1.2; x += 0.001) {
1314                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
1315                 DerivativeStructure rebuiltX = dsX.toDegrees().toRadians();
1316                 DerivativeStructure zero = rebuiltX.subtract(dsX);
1317                 for (int n = 0; n <= maxOrder; ++n) {
1318                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon);
1319                 }
1320             }
1321         }
1322     }
1323 
1324     @Test(expected=DimensionMismatchException.class)
1325     public void testComposeMismatchedDimensions() {
1326         new DerivativeStructure(1, 3, 0, 1.2).compose(new double[3]);
1327     }
1328 
1329     @Test
1330     public void testCompose() {
1331         double[] epsilon = new double[] { 1.0e-20, 5.0e-14, 2.0e-13, 3.0e-13, 2.0e-13, 1.0e-20 };
1332         PolynomialFunction poly =
1333                 new PolynomialFunction(new double[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 });
1334         for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
1335             PolynomialFunction[] p = new PolynomialFunction[maxOrder + 1];
1336             p[0] = poly;
1337             for (int i = 1; i <= maxOrder; ++i) {
1338                 p[i] = p[i - 1].polynomialDerivative();
1339             }
1340             for (double x = 0.1; x < 1.2; x += 0.001) {
1341                 DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
1342                 DerivativeStructure dsY1 = dsX.getField().getZero();
1343                 for (int i = poly.degree(); i >= 0; --i) {
1344                     dsY1 = dsY1.multiply(dsX).add(poly.getCoefficients()[i]);
1345                 }
1346                 double[] f = new double[maxOrder + 1];
1347                 for (int i = 0; i < f.length; ++i) {
1348                     f[i] = p[i].value(x);
1349                 }
1350                 DerivativeStructure dsY2 = dsX.compose(f);
1351                 DerivativeStructure zero = dsY1.subtract(dsY2);
1352                 for (int n = 0; n <= maxOrder; ++n) {
1353                     Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
1354                 }
1355             }
1356         }
1357     }
1358 
1359     @Test
1360     public void testField() {
1361         for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
1362             DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
1363             checkF0F1(x.getField().getZero(), 0.0, 0.0, 0.0, 0.0);
1364             checkF0F1(x.getField().getOne(), 1.0, 0.0, 0.0, 0.0);
1365             Assert.assertEquals(maxOrder, x.getField().getZero().getOrder());
1366             Assert.assertEquals(3, x.getField().getZero().getFreeParameters());
1367             Assert.assertEquals(DerivativeStructure.class, x.getField().getRuntimeClass());
1368         }
1369     }
1370 
1371     @Test
1372     public void testOneParameterConstructor() {
1373         double x = 1.2;
1374         double cos = JdkMath.cos(x);
1375         double sin = JdkMath.sin(x);
1376         DerivativeStructure yRef = new DerivativeStructure(1, 4, 0, x).cos();
1377         try {
1378             new DerivativeStructure(1, 4, 0.0, 0.0);
1379             Assert.fail("an exception should have been thrown");
1380         } catch (DimensionMismatchException dme) {
1381             // expected
1382         } catch (Exception e) {
1383             Assert.fail("wrong exceptionc caught " + e.getClass().getName());
1384         }
1385         double[] derivatives = new double[] { cos, -sin, -cos, sin, cos };
1386         DerivativeStructure y = new DerivativeStructure(1,  4, derivatives);
1387         checkEquals(yRef, y, 1.0e-15);
1388         TestUtils.assertEquals(derivatives, y.getAllDerivatives(), 1.0e-15);
1389     }
1390 
1391     @Test
1392     public void testOneOrderConstructor() {
1393         double x =  1.2;
1394         double y =  2.4;
1395         double z = 12.5;
1396         DerivativeStructure xRef = new DerivativeStructure(3, 1, 0, x);
1397         DerivativeStructure yRef = new DerivativeStructure(3, 1, 1, y);
1398         DerivativeStructure zRef = new DerivativeStructure(3, 1, 2, z);
1399         try {
1400             new DerivativeStructure(3, 1, x + y - z, 1.0, 1.0);
1401             Assert.fail("an exception should have been thrown");
1402         } catch (DimensionMismatchException dme) {
1403             // expected
1404         } catch (Exception e) {
1405             Assert.fail("wrong exceptionc caught " + e.getClass().getName());
1406         }
1407         double[] derivatives = new double[] { x + y - z, 1.0, 1.0, -1.0 };
1408         DerivativeStructure t = new DerivativeStructure(3, 1, derivatives);
1409         checkEquals(xRef.add(yRef.subtract(zRef)), t, 1.0e-15);
1410         TestUtils.assertEquals(derivatives, xRef.add(yRef.subtract(zRef)).getAllDerivatives(), 1.0e-15);
1411     }
1412 
1413     @Test
1414     public void testLinearCombination1DSDS() {
1415         final DerivativeStructure[] a = new DerivativeStructure[] {
1416             new DerivativeStructure(6, 1, 0, -1321008684645961.0 / 268435456.0),
1417             new DerivativeStructure(6, 1, 1, -5774608829631843.0 / 268435456.0),
1418             new DerivativeStructure(6, 1, 2, -7645843051051357.0 / 8589934592.0)
1419         };
1420         final DerivativeStructure[] b = new DerivativeStructure[] {
1421             new DerivativeStructure(6, 1, 3, -5712344449280879.0 / 2097152.0),
1422             new DerivativeStructure(6, 1, 4, -4550117129121957.0 / 2097152.0),
1423             new DerivativeStructure(6, 1, 5, 8846951984510141.0 / 131072.0)
1424         };
1425 
1426         final DerivativeStructure abSumInline = a[0].linearCombination(a[0], b[0], a[1], b[1], a[2], b[2]);
1427         final DerivativeStructure abSumArray = a[0].linearCombination(a, b);
1428 
1429         Assert.assertEquals(abSumInline.getValue(), abSumArray.getValue(), 0);
1430         Assert.assertEquals(-1.8551294182586248737720779899, abSumInline.getValue(), 1.0e-15);
1431         Assert.assertEquals(b[0].getValue(), abSumInline.getPartialDerivative(1, 0, 0, 0, 0, 0), 1.0e-15);
1432         Assert.assertEquals(b[1].getValue(), abSumInline.getPartialDerivative(0, 1, 0, 0, 0, 0), 1.0e-15);
1433         Assert.assertEquals(b[2].getValue(), abSumInline.getPartialDerivative(0, 0, 1, 0, 0, 0), 1.0e-15);
1434         Assert.assertEquals(a[0].getValue(), abSumInline.getPartialDerivative(0, 0, 0, 1, 0, 0), 1.0e-15);
1435         Assert.assertEquals(a[1].getValue(), abSumInline.getPartialDerivative(0, 0, 0, 0, 1, 0), 1.0e-15);
1436         Assert.assertEquals(a[2].getValue(), abSumInline.getPartialDerivative(0, 0, 0, 0, 0, 1), 1.0e-15);
1437     }
1438 
1439     @Test
1440     public void testLinearCombination1DoubleDS() {
1441         final double[] a = new double[] {
1442             -1321008684645961.0 / 268435456.0,
1443             -5774608829631843.0 / 268435456.0,
1444             -7645843051051357.0 / 8589934592.0
1445         };
1446         final DerivativeStructure[] b = new DerivativeStructure[] {
1447             new DerivativeStructure(3, 1, 0, -5712344449280879.0 / 2097152.0),
1448             new DerivativeStructure(3, 1, 1, -4550117129121957.0 / 2097152.0),
1449             new DerivativeStructure(3, 1, 2, 8846951984510141.0 / 131072.0)
1450         };
1451 
1452         final DerivativeStructure abSumInline = b[0].linearCombination(a[0], b[0],
1453                                                                        a[1], b[1],
1454                                                                        a[2], b[2]);
1455         final DerivativeStructure abSumArray = b[0].linearCombination(a, b);
1456 
1457         Assert.assertEquals(abSumInline.getValue(), abSumArray.getValue(), 0);
1458         Assert.assertEquals(-1.8551294182586248737720779899, abSumInline.getValue(), 1.0e-15);
1459         Assert.assertEquals(a[0], abSumInline.getPartialDerivative(1, 0, 0), 1.0e-15);
1460         Assert.assertEquals(a[1], abSumInline.getPartialDerivative(0, 1, 0), 1.0e-15);
1461         Assert.assertEquals(a[2], abSumInline.getPartialDerivative(0, 0, 1), 1.0e-15);
1462     }
1463 
1464     @Test
1465     public void testLinearCombination2DSDS() {
1466         // we compare accurate versus naive dot product implementations
1467         // on regular vectors (i.e. not extreme cases like in the previous test)
1468         UniformRandomProvider random = RandomSource.WELL_1024_A.create(0xc6af886975069f11L);
1469 
1470         for (int i = 0; i < 10000; ++i) {
1471             final DerivativeStructure[] u = new DerivativeStructure[4];
1472             final DerivativeStructure[] v = new DerivativeStructure[4];
1473             for (int j = 0; j < u.length; ++j) {
1474                 u[j] = new DerivativeStructure(u.length, 1, j, 1e17 * random.nextDouble());
1475                 v[j] = new DerivativeStructure(u.length, 1, 1e17 * random.nextDouble());
1476             }
1477 
1478             DerivativeStructure lin = u[0].linearCombination(u[0], v[0], u[1], v[1]);
1479             double ref = u[0].getValue() * v[0].getValue() +
1480                          u[1].getValue() * v[1].getValue();
1481             Assert.assertEquals(ref, lin.getValue(), 1.0e-15 * JdkMath.abs(ref));
1482             Assert.assertEquals(v[0].getValue(), lin.getPartialDerivative(1, 0, 0, 0), 1.0e-15 * JdkMath.abs(v[0].getValue()));
1483             Assert.assertEquals(v[1].getValue(), lin.getPartialDerivative(0, 1, 0, 0), 1.0e-15 * JdkMath.abs(v[1].getValue()));
1484 
1485             lin = u[0].linearCombination(u[0], v[0], u[1], v[1], u[2], v[2]);
1486             ref = u[0].getValue() * v[0].getValue() +
1487                   u[1].getValue() * v[1].getValue() +
1488                   u[2].getValue() * v[2].getValue();
1489             Assert.assertEquals(ref, lin.getValue(), 1.0e-15 * JdkMath.abs(ref));
1490             Assert.assertEquals(v[0].getValue(), lin.getPartialDerivative(1, 0, 0, 0), 1.0e-15 * JdkMath.abs(v[0].getValue()));
1491             Assert.assertEquals(v[1].getValue(), lin.getPartialDerivative(0, 1, 0, 0), 1.0e-15 * JdkMath.abs(v[1].getValue()));
1492             Assert.assertEquals(v[2].getValue(), lin.getPartialDerivative(0, 0, 1, 0), 1.0e-15 * JdkMath.abs(v[2].getValue()));
1493 
1494             lin = u[0].linearCombination(u[0], v[0], u[1], v[1], u[2], v[2], u[3], v[3]);
1495             ref = u[0].getValue() * v[0].getValue() +
1496                   u[1].getValue() * v[1].getValue() +
1497                   u[2].getValue() * v[2].getValue() +
1498                   u[3].getValue() * v[3].getValue();
1499             Assert.assertEquals(ref, lin.getValue(), 1.0e-15 * JdkMath.abs(ref));
1500             Assert.assertEquals(v[0].getValue(), lin.getPartialDerivative(1, 0, 0, 0), 1.0e-15 * JdkMath.abs(v[0].getValue()));
1501             Assert.assertEquals(v[1].getValue(), lin.getPartialDerivative(0, 1, 0, 0), 1.0e-15 * JdkMath.abs(v[1].getValue()));
1502             Assert.assertEquals(v[2].getValue(), lin.getPartialDerivative(0, 0, 1, 0), 1.0e-15 * JdkMath.abs(v[2].getValue()));
1503             Assert.assertEquals(v[3].getValue(), lin.getPartialDerivative(0, 0, 0, 1), 1.0e-15 * JdkMath.abs(v[3].getValue()));
1504         }
1505     }
1506 
1507     @Test
1508     public void testLinearCombination2DoubleDS() {
1509         // we compare accurate versus naive dot product implementations
1510         // on regular vectors (i.e. not extreme cases like in the previous test)
1511         UniformRandomProvider random = RandomSource.WELL_1024_A.create(0xc6af886975069f11L);
1512 
1513         for (int i = 0; i < 10000; ++i) {
1514             final double[] u = new double[4];
1515             final DerivativeStructure[] v = new DerivativeStructure[4];
1516             for (int j = 0; j < u.length; ++j) {
1517                 u[j] = 1e17 * random.nextDouble();
1518                 v[j] = new DerivativeStructure(u.length, 1, j, 1e17 * random.nextDouble());
1519             }
1520 
1521             DerivativeStructure lin = v[0].linearCombination(u[0], v[0], u[1], v[1]);
1522             double ref = u[0] * v[0].getValue() +
1523                          u[1] * v[1].getValue();
1524             Assert.assertEquals(ref, lin.getValue(), 1.0e-15 * JdkMath.abs(ref));
1525             Assert.assertEquals(u[0], lin.getPartialDerivative(1, 0, 0, 0), 1.0e-15 * JdkMath.abs(v[0].getValue()));
1526             Assert.assertEquals(u[1], lin.getPartialDerivative(0, 1, 0, 0), 1.0e-15 * JdkMath.abs(v[1].getValue()));
1527 
1528             lin = v[0].linearCombination(u[0], v[0], u[1], v[1], u[2], v[2]);
1529             ref = u[0] * v[0].getValue() +
1530                   u[1] * v[1].getValue() +
1531                   u[2] * v[2].getValue();
1532             Assert.assertEquals(ref, lin.getValue(), 1.0e-15 * JdkMath.abs(ref));
1533             Assert.assertEquals(u[0], lin.getPartialDerivative(1, 0, 0, 0), 1.0e-15 * JdkMath.abs(v[0].getValue()));
1534             Assert.assertEquals(u[1], lin.getPartialDerivative(0, 1, 0, 0), 1.0e-15 * JdkMath.abs(v[1].getValue()));
1535             Assert.assertEquals(u[2], lin.getPartialDerivative(0, 0, 1, 0), 1.0e-15 * JdkMath.abs(v[2].getValue()));
1536 
1537             lin = v[0].linearCombination(u[0], v[0], u[1], v[1], u[2], v[2], u[3], v[3]);
1538             ref = u[0] * v[0].getValue() +
1539                   u[1] * v[1].getValue() +
1540                   u[2] * v[2].getValue() +
1541                   u[3] * v[3].getValue();
1542             Assert.assertEquals(ref, lin.getValue(), 1.0e-15 * JdkMath.abs(ref));
1543             Assert.assertEquals(u[0], lin.getPartialDerivative(1, 0, 0, 0), 1.0e-15 * JdkMath.abs(v[0].getValue()));
1544             Assert.assertEquals(u[1], lin.getPartialDerivative(0, 1, 0, 0), 1.0e-15 * JdkMath.abs(v[1].getValue()));
1545             Assert.assertEquals(u[2], lin.getPartialDerivative(0, 0, 1, 0), 1.0e-15 * JdkMath.abs(v[2].getValue()));
1546             Assert.assertEquals(u[3], lin.getPartialDerivative(0, 0, 0, 1), 1.0e-15 * JdkMath.abs(v[3].getValue()));
1547         }
1548     }
1549 
1550     private void checkF0F1(DerivativeStructure ds, double value, double...derivatives) {
1551 
1552         // check dimension
1553         Assert.assertEquals(derivatives.length, ds.getFreeParameters());
1554 
1555         // check value, directly and also as 0th order derivative
1556         Assert.assertEquals(value, ds.getValue(), 1.0e-15);
1557         Assert.assertEquals(value, ds.getPartialDerivative(new int[ds.getFreeParameters()]), 1.0e-15);
1558 
1559         // check first order derivatives
1560         for (int i = 0; i < derivatives.length; ++i) {
1561             int[] orders = new int[derivatives.length];
1562             orders[i] = 1;
1563             Assert.assertEquals(derivatives[i], ds.getPartialDerivative(orders), 1.0e-15);
1564         }
1565     }
1566 
1567     private void checkEquals(DerivativeStructure ds1, DerivativeStructure ds2, double epsilon) {
1568 
1569         // check dimension
1570         Assert.assertEquals(ds1.getFreeParameters(), ds2.getFreeParameters());
1571         Assert.assertEquals(ds1.getOrder(), ds2.getOrder());
1572 
1573         int[] derivatives = new int[ds1.getFreeParameters()];
1574         int sum = 0;
1575         while (true) {
1576 
1577             if (sum <= ds1.getOrder()) {
1578                 Assert.assertEquals(ds1.getPartialDerivative(derivatives),
1579                                     ds2.getPartialDerivative(derivatives),
1580                                     epsilon);
1581             }
1582 
1583             boolean increment = true;
1584             sum = 0;
1585             for (int i = derivatives.length - 1; i >= 0; --i) {
1586                 if (increment) {
1587                     if (derivatives[i] == ds1.getOrder()) {
1588                         derivatives[i] = 0;
1589                     } else {
1590                         derivatives[i]++;
1591                         increment = false;
1592                     }
1593                 }
1594                 sum += derivatives[i];
1595             }
1596             if (increment) {
1597                 return;
1598             }
1599         }
1600     }
1601 }