Completed support fo asin, acos and atan in DSCompiler.

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1371805 13f79535-47bb-0310-9956-ffa450edef68

src/main/java/org/apache/commons/math3/analysis/differentiation/DSCompiler.java

@@ -1097,10 +1097,39 @@ public class DSCompiler {
         final double x = operand[operandOffset];
         function[0] = FastMath.acos(x);
         if (order > 0) {
-            function[1] = -1.0 / FastMath.sqrt(1 - x * x);
-            for (int i = 2; i <= order; ++i) {
-                // TODO compute higher order derivatives
-                function[i] = Double.NaN;
+            // the nth order derivative of acos has the form:
+            // dn(acos(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
+            // where P_n(x) is a degree n-1 polynomial with same parity as n-1
+            // P_1(x) = -1, P_2(x) = -x, P_3(x) = -2x^2 - 1 ...
+            // the general recurrence relation for P_n is:
+            // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
+            // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
+            final double[] p = new double[order];
+            p[0] = -1;
+            final double x2    = x * x;
+            final double f     = 1.0 / (1 - x2);
+            double coeff = FastMath.sqrt(f);
+            function[1] = coeff * p[0];
+            for (int n = 2; n <= order; ++n) {
+
+                // update and evaluate polynomial P_n(x)
+                double v = 0;
+                p[n - 1] = (n - 1) * p[n - 2];
+                for (int k = n - 1; k >= 0; k -= 2) {
+                    v = v * x2 + p[k];
+                    if (k > 2) {
+                        p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
+                    } else if (k == 2) {
+                        p[0] = p[1];
+                    }
+                }
+                if ((n & 0x1) == 0) {
+                    v *= x;
+                }
+
+                coeff *= f;
+                function[n] = coeff * v;
+
             }
         }
 
@@ -1125,10 +1154,39 @@ public class DSCompiler {
         final double x = operand[operandOffset];
         function[0] = FastMath.asin(x);
         if (order > 0) {
-            function[1] = 1.0 / FastMath.sqrt(1 - x * x);
-            for (int i = 2; i <= order; ++i) {
-                // TODO compute higher order derivatives
-                function[i] = Double.NaN;
+            // the nth order derivative of asin has the form:
+            // dn(asin(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
+            // where P_n(x) is a degree n-1 polynomial with same parity as n-1
+            // P_1(x) = 1, P_2(x) = x, P_3(x) = 2x^2 + 1 ...
+            // the general recurrence relation for P_n is:
+            // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
+            // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
+            final double[] p = new double[order];
+            p[0] = 1;
+            final double x2    = x * x;
+            final double f     = 1.0 / (1 - x2);
+            double coeff = FastMath.sqrt(f);
+            function[1] = coeff * p[0];
+            for (int n = 2; n <= order; ++n) {
+
+                // update and evaluate polynomial P_n(x)
+                double v = 0;
+                p[n - 1] = (n - 1) * p[n - 2];
+                for (int k = n - 1; k >= 0; k -= 2) {
+                    v = v * x2 + p[k];
+                    if (k > 2) {
+                        p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
+                    } else if (k == 2) {
+                        p[0] = p[1];
+                    }
+                }
+                if ((n & 0x1) == 0) {
+                    v *= x;
+                }
+
+                coeff *= f;
+                function[n] = coeff * v;
+
             }
         }
 
@@ -1153,10 +1211,39 @@ public class DSCompiler {
         final double x = operand[operandOffset];
         function[0] = FastMath.atan(x);
         if (order > 0) {
-            function[1] = 1.0 / (1 + x * x);
-            for (int i = 2; i <= order; ++i) {
-                // TODO compute higher order derivatives
-                function[i] = Double.NaN;
+            // the nth order derivative of atan has the form:
+            // dn(atan(x)/dxn = Q_n(x) / (1 + x^2)^n
+            // where Q_n(x) is a degree n-1 polynomial with same parity as n-1
+            // Q_1(x) = 1, Q_2(x) = -2x, Q_3(x) = 6x^2 - 2 ...
+            // the general recurrence relation for Q_n is:
+            // Q_n(x) = (1+x^2) Q_(n-1)'(x) - 2(n-1) x Q_(n-1)(x)
+            // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array
+            final double[] q = new double[order];
+            q[0] = 1;
+            final double x2    = x * x;
+            final double f     = 1.0 / (1 + x2);
+            double coeff = f;
+            function[1] = coeff * q[0];
+            for (int n = 2; n <= order; ++n) {
+
+                // update and evaluate polynomial Q_n(x)
+                double v = 0;
+                q[n - 1] = -n * q[n - 2];
+                for (int k = n - 1; k >= 0; k -= 2) {
+                    v = v * x2 + q[k];
+                    if (k > 2) {
+                        q[k - 2] = (k - 1) * q[k - 1] + (k - 1 - 2 * n) * q[k - 3];
+                    } else if (k == 2) {
+                        q[0] = q[1];
+                    }
+                }
+                if ((n & 0x1) == 0) {
+                    v *= x;
+                }
+
+                coeff *= f;
+                function[n] = coeff * v;
+
             }
         }
 

src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java

@@ -459,6 +459,51 @@ public class DerivativeStructureTest {
         }
     }
 
+    @Test
+    public void testSinAsin() {
+        double[] epsilon = new double[] { 3.0e-16, 5.0e-16, 3.0e-15, 2.0e-14, 4.0e-13 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
+                DerivativeStructure rebuiltX = dsX.sin().asin();
+                DerivativeStructure zero = rebuiltX.subtract(dsX);
+                for (int n = 0; n <= maxOrder; ++n) {
+                    Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testCosAcos() {
+        double[] epsilon = new double[] { 6.0e-16, 6.0e-15, 2.0e-13, 4.0e-12, 2.0e-10 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
+                DerivativeStructure rebuiltX = dsX.cos().acos();
+                DerivativeStructure zero = rebuiltX.subtract(dsX);
+                for (int n = 0; n <= maxOrder; ++n) {
+                    Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testTanAtan() {
+        double[] epsilon = new double[] { 6.0e-17, 2.0e-16, 2.0e-15, 4.0e-14, 2.0e-12 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
+                DerivativeStructure rebuiltX = dsX.tan().atan();
+                DerivativeStructure zero = rebuiltX.subtract(dsX);
+                for (int n = 0; n <= maxOrder; ++n) {
+                    Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
     @Test
     public void testTangentDefinition() {
         double[] epsilon = new double[] { 5.0e-16, 2.0e-15, 3.0e-14, 5.0e-13, 2.0e-11 };