Java tutorial
/* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.apache.commons.math3.analysis.function; import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction; import org.apache.commons.math3.analysis.FunctionUtils; import org.apache.commons.math3.analysis.UnivariateFunction; import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; import org.apache.commons.math3.util.FastMath; /** * <a href="http://en.wikipedia.org/wiki/Sinc_function">Sinc</a> function, * defined by * <pre><code> * sinc(x) = 1 if x = 0, * sin(x) / x otherwise. * </code></pre> * * @since 3.0 * @version $Id: Sinc.java 1383441 2012-09-11 14:56:39Z luc $ */ public class Sinc implements UnivariateDifferentiableFunction, DifferentiableUnivariateFunction { /** * Value below which the computations are done using Taylor series. * <p> * The Taylor series for sinc even order derivatives are: * <pre> * d^(2n)sinc/dx^(2n) = Sum_(k>=0) (-1)^(n+k) / ((2k)!(2n+2k+1)) x^(2k) * = (-1)^n [ 1/(2n+1) - x^2/(4n+6) + x^4/(48n+120) - x^6/(1440n+5040) + O(x^8) ] * </pre> * </p> * <p> * The Taylor series for sinc odd order derivatives are: * <pre> * d^(2n+1)sinc/dx^(2n+1) = Sum_(k>=0) (-1)^(n+k+1) / ((2k+1)!(2n+2k+3)) x^(2k+1) * = (-1)^(n+1) [ x/(2n+3) - x^3/(12n+30) + x^5/(240n+840) - x^7/(10080n+45360) + O(x^9) ] * </pre> * </p> * <p> * So the ratio of the fourth term with respect to the first term * is always smaller than x^6/720, for all derivative orders. * This implies that neglecting this term and using only the first three terms induces * a relative error bounded by x^6/720. The SHORTCUT value is chosen such that this * relative error is below double precision accuracy when |x| <= SHORTCUT. * </p> */ private static final double SHORTCUT = 6.0e-3; /** For normalized sinc function. */ private final boolean normalized; /** * The sinc function, {@code sin(x) / x}. */ public Sinc() { this(false); } /** * Instantiates the sinc function. * * @param normalized If {@code true}, the function is * <code> sin(πx) / πx</code>, otherwise {@code sin(x) / x}. */ public Sinc(boolean normalized) { this.normalized = normalized; } /** {@inheritDoc} */ public double value(final double x) { final double scaledX = normalized ? FastMath.PI * x : x; if (FastMath.abs(scaledX) <= SHORTCUT) { // use Taylor series final double scaledX2 = scaledX * scaledX; return ((scaledX2 - 20) * scaledX2 + 120) / 120; } else { // use definition expression return FastMath.sin(scaledX) / scaledX; } } /** {@inheritDoc} * @deprecated as of 3.1, replaced by {@link #value(DerivativeStructure)} */ @Deprecated public UnivariateFunction derivative() { return FunctionUtils.toDifferentiableUnivariateFunction(this).derivative(); } /** {@inheritDoc} * @since 3.1 */ public DerivativeStructure value(final DerivativeStructure t) { final double scaledX = (normalized ? FastMath.PI : 1) * t.getValue(); final double scaledX2 = scaledX * scaledX; double[] f = new double[t.getOrder() + 1]; if (FastMath.abs(scaledX) <= SHORTCUT) { for (int i = 0; i < f.length; ++i) { final int k = i / 2; if ((i & 0x1) == 0) { // even derivation order f[i] = (((k & 0x1) == 0) ? 1 : -1) * (1.0 / (i + 1) - scaledX2 * (1.0 / (2 * i + 6) - scaledX2 / (24 * i + 120))); } else { // odd derivation order f[i] = (((k & 0x1) == 0) ? -scaledX : scaledX) * (1.0 / (i + 2) - scaledX2 * (1.0 / (6 * i + 24) - scaledX2 / (120 * i + 720))); } } } else { final double inv = 1 / scaledX; final double cos = FastMath.cos(scaledX); final double sin = FastMath.sin(scaledX); f[0] = inv * sin; // the nth order derivative of sinc has the form: // dn(sinc(x)/dxn = [S_n(x) sin(x) + C_n(x) cos(x)] / x^(n+1) // where S_n(x) is an even polynomial with degree n-1 or n (depending on parity) // and C_n(x) is an odd polynomial with degree n-1 or n (depending on parity) // S_0(x) = 1, S_1(x) = -1, S_2(x) = -x^2 + 2, S_3(x) = 3x^2 - 6... // C_0(x) = 0, C_1(x) = x, C_2(x) = -2x, C_3(x) = -x^3 + 6x... // the general recurrence relations for S_n and C_n are: // S_n(x) = x S_(n-1)'(x) - n S_(n-1)(x) - x C_(n-1)(x) // C_n(x) = x C_(n-1)'(x) - n C_(n-1)(x) + x S_(n-1)(x) // as per polynomials parity, we can store both S_n and C_n in the same array final double[] sc = new double[f.length]; sc[0] = 1; double coeff = inv; for (int n = 1; n < f.length; ++n) { double s = 0; double c = 0; // update and evaluate polynomials S_n(x) and C_n(x) final int kStart; if ((n & 0x1) == 0) { // even derivation order, S_n is degree n and C_n is degree n-1 sc[n] = 0; kStart = n; } else { // odd derivation order, S_n is degree n-1 and C_n is degree n sc[n] = sc[n - 1]; c = sc[n]; kStart = n - 1; } // in this loop, k is always even for (int k = kStart; k > 1; k -= 2) { // sine part sc[k] = (k - n) * sc[k] - sc[k - 1]; s = s * scaledX2 + sc[k]; // cosine part sc[k - 1] = (k - 1 - n) * sc[k - 1] + sc[k - 2]; c = c * scaledX2 + sc[k - 1]; } sc[0] *= -n; s = s * scaledX2 + sc[0]; coeff *= inv; f[n] = coeff * (s * sin + c * scaledX * cos); } } if (normalized) { double scale = FastMath.PI; for (int i = 1; i < f.length; ++i) { f[i] *= scale; scale *= FastMath.PI; } } return t.compose(f); } }