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(πx) / π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 }