    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(&pi;x) / &pi;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    }