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
018package org.apache.commons.math4.analysis.function;
019
020import org.apache.commons.math4.analysis.differentiation.DerivativeStructure;
021import org.apache.commons.math4.analysis.differentiation.UnivariateDifferentiableFunction;
022import org.apache.commons.math4.exception.DimensionMismatchException;
023import org.apache.commons.math4.util.FastMath;
024
025/**
026 * <a href="http://en.wikipedia.org/wiki/Sinc_function">Sinc</a> function,
027 * defined by
028 * <pre><code>
029 *   sinc(x) = 1            if x = 0,
030 *             sin(x) / x   otherwise.
031 * </code></pre>
032 *
033 * @since 3.0
034 */
035public class Sinc implements UnivariateDifferentiableFunction {
036    /**
037     * Value below which the computations are done using Taylor series.
038     * <p>
039     * The Taylor series for sinc even order derivatives are:
040     * <pre>
041     * d^(2n)sinc/dx^(2n)     = Sum_(k>=0) (-1)^(n+k) / ((2k)!(2n+2k+1)) x^(2k)
042     *                        = (-1)^n     [ 1/(2n+1) - x^2/(4n+6) + x^4/(48n+120) - x^6/(1440n+5040) + O(x^8) ]
043     * </pre>
044     * </p>
045     * <p>
046     * The Taylor series for sinc odd order derivatives are:
047     * <pre>
048     * d^(2n+1)sinc/dx^(2n+1) = Sum_(k>=0) (-1)^(n+k+1) / ((2k+1)!(2n+2k+3)) x^(2k+1)
049     *                        = (-1)^(n+1) [ x/(2n+3) - x^3/(12n+30) + x^5/(240n+840) - x^7/(10080n+45360) + O(x^9) ]
050     * </pre>
051     * </p>
052     * <p>
053     * So the ratio of the fourth term with respect to the first term
054     * is always smaller than x^6/720, for all derivative orders.
055     * This implies that neglecting this term and using only the first three terms induces
056     * a relative error bounded by x^6/720. The SHORTCUT value is chosen such that this
057     * relative error is below double precision accuracy when |x| <= SHORTCUT.
058     * </p>
059     */
060    private static final double SHORTCUT = 6.0e-3;
061    /** For normalized sinc function. */
062    private final boolean normalized;
063
064    /**
065     * The sinc function, {@code sin(x) / x}.
066     */
067    public Sinc() {
068        this(false);
069    }
070
071    /**
072     * Instantiates the sinc function.
073     *
074     * @param normalized If {@code true}, the function is
075     * <code> sin(&pi;x) / &pi;x</code>, otherwise {@code sin(x) / x}.
076     */
077    public Sinc(boolean normalized) {
078        this.normalized = normalized;
079    }
080
081    /** {@inheritDoc} */
082    @Override
083    public double value(final double x) {
084        final double scaledX = normalized ? FastMath.PI * x : x;
085        if (FastMath.abs(scaledX) <= SHORTCUT) {
086            // use Taylor series
087            final double scaledX2 = scaledX * scaledX;
088            return ((scaledX2 - 20) * scaledX2 + 120) / 120;
089        } else {
090            // use definition expression
091            return FastMath.sin(scaledX) / scaledX;
092        }
093    }
094
095    /** {@inheritDoc}
096     * @since 3.1
097     */
098    @Override
099    public DerivativeStructure value(final DerivativeStructure t)
100        throws DimensionMismatchException {
101
102        final double scaledX  = (normalized ? FastMath.PI : 1) * t.getValue();
103        final double scaledX2 = scaledX * scaledX;
104
105        double[] f = new double[t.getOrder() + 1];
106
107        if (FastMath.abs(scaledX) <= SHORTCUT) {
108
109            for (int i = 0; i < f.length; ++i) {
110                final int k = i / 2;
111                if ((i & 0x1) == 0) {
112                    // even derivation order
113                    f[i] = (((k & 0x1) == 0) ? 1 : -1) *
114                           (1.0 / (i + 1) - scaledX2 * (1.0 / (2 * i + 6) - scaledX2 / (24 * i + 120)));
115                } else {
116                    // odd derivation order
117                    f[i] = (((k & 0x1) == 0) ? -scaledX : scaledX) *
118                           (1.0 / (i + 2) - scaledX2 * (1.0 / (6 * i + 24) - scaledX2 / (120 * i + 720)));
119                }
120            }
121
122        } else {
123
124            final double inv = 1 / scaledX;
125            final double cos = FastMath.cos(scaledX);
126            final double sin = FastMath.sin(scaledX);
127
128            f[0] = inv * sin;
129
130            // the nth order derivative of sinc has the form:
131            // dn(sinc(x)/dxn = [S_n(x) sin(x) + C_n(x) cos(x)] / x^(n+1)
132            // where S_n(x) is an even polynomial with degree n-1 or n (depending on parity)
133            // and C_n(x) is an odd polynomial with degree n-1 or n (depending on parity)
134            // S_0(x) = 1, S_1(x) = -1, S_2(x) = -x^2 + 2, S_3(x) = 3x^2 - 6...
135            // C_0(x) = 0, C_1(x) = x, C_2(x) = -2x, C_3(x) = -x^3 + 6x...
136            // the general recurrence relations for S_n and C_n are:
137            // S_n(x) = x S_(n-1)'(x) - n S_(n-1)(x) - x C_(n-1)(x)
138            // C_n(x) = x C_(n-1)'(x) - n C_(n-1)(x) + x S_(n-1)(x)
139            // as per polynomials parity, we can store both S_n and C_n in the same array
140            final double[] sc = new double[f.length];
141            sc[0] = 1;
142
143            double coeff = inv;
144            for (int n = 1; n < f.length; ++n) {
145
146                double s = 0;
147                double c = 0;
148
149                // update and evaluate polynomials S_n(x) and C_n(x)
150                final int kStart;
151                if ((n & 0x1) == 0) {
152                    // even derivation order, S_n is degree n and C_n is degree n-1
153                    sc[n] = 0;
154                    kStart = n;
155                } else {
156                    // odd derivation order, S_n is degree n-1 and C_n is degree n
157                    sc[n] = sc[n - 1];
158                    c = sc[n];
159                    kStart = n - 1;
160                }
161
162                // in this loop, k is always even
163                for (int k = kStart; k > 1; k -= 2) {
164
165                    // sine part
166                    sc[k]     = (k - n) * sc[k] - sc[k - 1];
167                    s         = s * scaledX2 + sc[k];
168
169                    // cosine part
170                    sc[k - 1] = (k - 1 - n) * sc[k - 1] + sc[k -2];
171                    c         = c * scaledX2 + sc[k - 1];
172
173                }
174                sc[0] *= -n;
175                s      = s * scaledX2 + sc[0];
176
177                coeff *= inv;
178                f[n]   = coeff * (s * sin + c * scaledX * cos);
179
180            }
181
182        }
183
184        if (normalized) {
185            double scale = FastMath.PI;
186            for (int i = 1; i < f.length; ++i) {
187                f[i]  *= scale;
188                scale *= FastMath.PI;
189            }
190        }
191
192        return t.compose(f);
193
194    }
195
196}