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