001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math3.analysis.function;
019
020 import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction;
021 import org.apache.commons.math3.analysis.FunctionUtils;
022 import org.apache.commons.math3.analysis.UnivariateFunction;
023 import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
024 import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction;
025 import org.apache.commons.math3.util.FastMath;
026
027 /**
028 * <a href="http://en.wikipedia.org/wiki/Sinc_function">Sinc</a> function,
029 * defined by
030 * <pre><code>
031 * sinc(x) = 1 if x = 0,
032 * sin(x) / x otherwise.
033 * </code></pre>
034 *
035 * @since 3.0
036 * @version $Id: Sinc.java 1383441 2012-09-11 14:56:39Z luc $
037 */
038 public class Sinc implements UnivariateDifferentiableFunction, DifferentiableUnivariateFunction {
039 /**
040 * Value below which the computations are done using Taylor series.
041 * <p>
042 * The Taylor series for sinc even order derivatives are:
043 * <pre>
044 * d^(2n)sinc/dx^(2n) = Sum_(k>=0) (-1)^(n+k) / ((2k)!(2n+2k+1)) x^(2k)
045 * = (-1)^n [ 1/(2n+1) - x^2/(4n+6) + x^4/(48n+120) - x^6/(1440n+5040) + O(x^8) ]
046 * </pre>
047 * </p>
048 * <p>
049 * The Taylor series for sinc odd order derivatives are:
050 * <pre>
051 * d^(2n+1)sinc/dx^(2n+1) = Sum_(k>=0) (-1)^(n+k+1) / ((2k+1)!(2n+2k+3)) x^(2k+1)
052 * = (-1)^(n+1) [ x/(2n+3) - x^3/(12n+30) + x^5/(240n+840) - x^7/(10080n+45360) + O(x^9) ]
053 * </pre>
054 * </p>
055 * <p>
056 * So the ratio of the fourth term with respect to the first term
057 * is always smaller than x^6/720, for all derivative orders.
058 * This implies that neglecting this term and using only the first three terms induces
059 * a relative error bounded by x^6/720. The SHORTCUT value is chosen such that this
060 * relative error is below double precision accuracy when |x| <= SHORTCUT.
061 * </p>
062 */
063 private static final double SHORTCUT = 6.0e-3;
064 /** For normalized sinc function. */
065 private final boolean normalized;
066
067 /**
068 * The sinc function, {@code sin(x) / x}.
069 */
070 public Sinc() {
071 this(false);
072 }
073
074 /**
075 * Instantiates the sinc function.
076 *
077 * @param normalized If {@code true}, the function is
078 * <code> sin(πx) / πx</code>, otherwise {@code sin(x) / x}.
079 */
080 public Sinc(boolean normalized) {
081 this.normalized = normalized;
082 }
083
084 /** {@inheritDoc} */
085 public double value(final double x) {
086 final double scaledX = normalized ? FastMath.PI * x : x;
087 if (FastMath.abs(scaledX) <= SHORTCUT) {
088 // use Taylor series
089 final double scaledX2 = scaledX * scaledX;
090 return ((scaledX2 - 20) * scaledX2 + 120) / 120;
091 } else {
092 // use definition expression
093 return FastMath.sin(scaledX) / scaledX;
094 }
095 }
096
097 /** {@inheritDoc}
098 * @deprecated as of 3.1, replaced by {@link #value(DerivativeStructure)}
099 */
100 @Deprecated
101 public UnivariateFunction derivative() {
102 return FunctionUtils.toDifferentiableUnivariateFunction(this).derivative();
103 }
104
105 /** {@inheritDoc}
106 * @since 3.1
107 */
108 public DerivativeStructure value(final DerivativeStructure t) {
109
110 final double scaledX = (normalized ? FastMath.PI : 1) * t.getValue();
111 final double scaledX2 = scaledX * scaledX;
112
113 double[] f = new double[t.getOrder() + 1];
114
115 if (FastMath.abs(scaledX) <= SHORTCUT) {
116
117 for (int i = 0; i < f.length; ++i) {
118 final int k = i / 2;
119 if ((i & 0x1) == 0) {
120 // even derivation order
121 f[i] = (((k & 0x1) == 0) ? 1 : -1) *
122 (1.0 / (i + 1) - scaledX2 * (1.0 / (2 * i + 6) - scaledX2 / (24 * i + 120)));
123 } else {
124 // odd derivation order
125 f[i] = (((k & 0x1) == 0) ? -scaledX : scaledX) *
126 (1.0 / (i + 2) - scaledX2 * (1.0 / (6 * i + 24) - scaledX2 / (120 * i + 720)));
127 }
128 }
129
130 } else {
131
132 final double inv = 1 / scaledX;
133 final double cos = FastMath.cos(scaledX);
134 final double sin = FastMath.sin(scaledX);
135
136 f[0] = inv * sin;
137
138 // the nth order derivative of sinc has the form:
139 // dn(sinc(x)/dxn = [S_n(x) sin(x) + C_n(x) cos(x)] / x^(n+1)
140 // where S_n(x) is an even polynomial with degree n-1 or n (depending on parity)
141 // and C_n(x) is an odd polynomial with degree n-1 or n (depending on parity)
142 // S_0(x) = 1, S_1(x) = -1, S_2(x) = -x^2 + 2, S_3(x) = 3x^2 - 6...
143 // C_0(x) = 0, C_1(x) = x, C_2(x) = -2x, C_3(x) = -x^3 + 6x...
144 // the general recurrence relations for S_n and C_n are:
145 // S_n(x) = x S_(n-1)'(x) - n S_(n-1)(x) - x C_(n-1)(x)
146 // C_n(x) = x C_(n-1)'(x) - n C_(n-1)(x) + x S_(n-1)(x)
147 // as per polynomials parity, we can store both S_n and C_n in the same array
148 final double[] sc = new double[f.length];
149 sc[0] = 1;
150
151 double coeff = inv;
152 for (int n = 1; n < f.length; ++n) {
153
154 double s = 0;
155 double c = 0;
156
157 // update and evaluate polynomials S_n(x) and C_n(x)
158 final int kStart;
159 if ((n & 0x1) == 0) {
160 // even derivation order, S_n is degree n and C_n is degree n-1
161 sc[n] = 0;
162 kStart = n;
163 } else {
164 // odd derivation order, S_n is degree n-1 and C_n is degree n
165 sc[n] = sc[n - 1];
166 c = sc[n];
167 kStart = n - 1;
168 }
169
170 // in this loop, k is always even
171 for (int k = kStart; k > 1; k -= 2) {
172
173 // sine part
174 sc[k] = (k - n) * sc[k] - sc[k - 1];
175 s = s * scaledX2 + sc[k];
176
177 // cosine part
178 sc[k - 1] = (k - 1 - n) * sc[k - 1] + sc[k -2];
179 c = c * scaledX2 + sc[k - 1];
180
181 }
182 sc[0] *= -n;
183 s = s * scaledX2 + sc[0];
184
185 coeff *= inv;
186 f[n] = coeff * (s * sin + c * scaledX * cos);
187
188 }
189
190 }
191
192 if (normalized) {
193 double scale = FastMath.PI;
194 for (int i = 1; i < f.length; ++i) {
195 f[i] *= scale;
196 scale *= FastMath.PI;
197 }
198 }
199
200 return t.compose(f);
201
202 }
203
204 }