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.function;
24  
25  import org.hipparchus.analysis.differentiation.Derivative;
26  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
27  import org.hipparchus.exception.MathIllegalArgumentException;
28  import org.hipparchus.util.FastMath;
29  import org.hipparchus.util.SinCos;
30  
31  /**
32   * <a href="http://en.wikipedia.org/wiki/Sinc_function">Sinc</a> function,
33   * defined by
34   * <pre><code>
35   *   sinc(x) = 1            if x = 0,
36   *             sin(x) / x   otherwise.
37   * </code></pre>
38   *
39   */
40  public class Sinc implements UnivariateDifferentiableFunction {
41      /**
42       * Value below which the computations are done using Taylor series.
43       * <p>
44       * The Taylor series for sinc even order derivatives are:
45       * <pre>
46       * d^(2n)sinc/dx^(2n)     = Sum_(k>=0) (-1)^(n+k) / ((2k)!(2n+2k+1)) x^(2k)
47       *                        = (-1)^n     [ 1/(2n+1) - x^2/(4n+6) + x^4/(48n+120) - x^6/(1440n+5040) + O(x^8) ]
48       * </pre>
49       * </p>
50       * <p>
51       * The Taylor series for sinc odd order derivatives are:
52       * <pre>
53       * d^(2n+1)sinc/dx^(2n+1) = Sum_(k>=0) (-1)^(n+k+1) / ((2k+1)!(2n+2k+3)) x^(2k+1)
54       *                        = (-1)^(n+1) [ x/(2n+3) - x^3/(12n+30) + x^5/(240n+840) - x^7/(10080n+45360) + O(x^9) ]
55       * </pre>
56       * </p>
57       * <p>
58       * So the ratio of the fourth term with respect to the first term
59       * is always smaller than x^6/720, for all derivative orders.
60       * This implies that neglecting this term and using only the first three terms induces
61       * a relative error bounded by x^6/720. The SHORTCUT value is chosen such that this
62       * relative error is below double precision accuracy when |x| <= SHORTCUT.
63       * </p>
64       */
65      private static final double SHORTCUT = 6.0e-3;
66      /** For normalized sinc function. */
67      private final boolean normalized;
68  
69      /**
70       * The sinc function, {@code sin(x) / x}.
71       */
72      public Sinc() {
73          this(false);
74      }
75  
76      /**
77       * Instantiates the sinc function.
78       *
79       * @param normalized If {@code true}, the function is
80       * <code> sin(&pi;x) / &pi;x</code>, otherwise {@code sin(x) / x}.
81       */
82      public Sinc(boolean normalized) {
83          this.normalized = normalized;
84      }
85  
86      /** {@inheritDoc} */
87      @Override
88      public double value(final double x) {
89          final double scaledX = normalized ? FastMath.PI * x : x;
90          if (FastMath.abs(scaledX) <= SHORTCUT) {
91              // use Taylor series
92              final double scaledX2 = scaledX * scaledX;
93              return ((scaledX2 - 20) * scaledX2 + 120) / 120;
94          } else {
95              // use definition expression
96              return FastMath.sin(scaledX) / scaledX;
97          }
98      }
99  
100     /** {@inheritDoc}
101      */
102     @Override
103     public <T extends Derivative<T>> T value(T t)
104         throws MathIllegalArgumentException {
105 
106         final double scaledX   = (normalized ? FastMath.PI : 1) * t.getValue();
107         final double scaledX2  = scaledX * scaledX;
108 
109         double[] f = new double[t.getOrder() + 1];
110 
111         if (FastMath.abs(scaledX) <= SHORTCUT) {
112 
113             for (int i = 0; i < f.length; ++i) {
114                 final int k = i / 2;
115                 if ((i & 0x1) == 0) {
116                     // even derivation order
117                     f[i] = (((k & 0x1) == 0) ? 1 : -1) *
118                            (1.0 / (i + 1) - scaledX2 * (1.0 / (2 * i + 6) - scaledX2 / (24 * i + 120)));
119                 } else {
120                     // odd derivation order
121                     f[i] = (((k & 0x1) == 0) ? -scaledX : scaledX) *
122                            (1.0 / (i + 2) - scaledX2 * (1.0 / (6 * i + 24) - scaledX2 / (120 * i + 720)));
123                 }
124             }
125 
126         } else {
127 
128             final double inv    = 1 / scaledX;
129             final SinCos sinCos = FastMath.sinCos(scaledX);
130 
131             f[0] = inv * sinCos.sin();
132 
133             // the nth order derivative of sinc has the form:
134             // dn(sinc(x)/dxn = [S_n(x) sin(x) + C_n(x) cos(x)] / x^(n+1)
135             // where S_n(x) is an even polynomial with degree n-1 or n (depending on parity)
136             // and C_n(x) is an odd polynomial with degree n-1 or n (depending on parity)
137             // S_0(x) = 1, S_1(x) = -1, S_2(x) = -x^2 + 2, S_3(x) = 3x^2 - 6...
138             // C_0(x) = 0, C_1(x) = x, C_2(x) = -2x, C_3(x) = -x^3 + 6x...
139             // the general recurrence relations for S_n and C_n are:
140             // S_n(x) = x S_(n-1)'(x) - n S_(n-1)(x) - x C_(n-1)(x)
141             // C_n(x) = x C_(n-1)'(x) - n C_(n-1)(x) + x S_(n-1)(x)
142             // as per polynomials parity, we can store both S_n and C_n in the same array
143             final double[] sc = new double[f.length];
144             sc[0] = 1;
145 
146             double coeff = inv;
147             for (int n = 1; n < f.length; ++n) {
148 
149                 double s = 0;
150                 double c = 0;
151 
152                 // update and evaluate polynomials S_n(x) and C_n(x)
153                 final int kStart;
154                 if ((n & 0x1) == 0) {
155                     // even derivation order, S_n is degree n and C_n is degree n-1
156                     sc[n] = 0;
157                     kStart = n;
158                 } else {
159                     // odd derivation order, S_n is degree n-1 and C_n is degree n
160                     sc[n] = sc[n - 1];
161                     c = sc[n];
162                     kStart = n - 1;
163                 }
164 
165                 // in this loop, k is always even
166                 for (int k = kStart; k > 1; k -= 2) {
167 
168                     // sine part
169                     sc[k]     = (k - n) * sc[k] - sc[k - 1];
170                     s         = s * scaledX2 + sc[k];
171 
172                     // cosine part
173                     sc[k - 1] = (k - 1 - n) * sc[k - 1] + sc[k -2];
174                     c         = c * scaledX2 + sc[k - 1];
175 
176                 }
177                 sc[0] *= -n;
178                 s      = s * scaledX2 + sc[0];
179 
180                 coeff *= inv;
181                 f[n]   = coeff * (s * sinCos.sin() + c * scaledX * sinCos.cos());
182 
183             }
184 
185         }
186 
187         if (normalized) {
188             double scale = FastMath.PI;
189             for (int i = 1; i < f.length; ++i) {
190                 f[i]  *= scale;
191                 scale *= FastMath.PI;
192             }
193         }
194 
195         return t.compose(f);
196 
197     }
198 
199 }