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.math3.analysis.function; 019 020import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction; 021import org.apache.commons.math3.analysis.FunctionUtils; 022import org.apache.commons.math3.analysis.UnivariateFunction; 023import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; 024import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; 025import org.apache.commons.math3.exception.DimensionMismatchException; 026import org.apache.commons.math3.util.FastMath; 027 028/** 029 * <a href="http://en.wikipedia.org/wiki/Sinc_function">Sinc</a> function, 030 * defined by 031 * <pre><code> 032 * sinc(x) = 1 if x = 0, 033 * sin(x) / x otherwise. 034 * </code></pre> 035 * 036 * @since 3.0 037 */ 038public 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 throws DimensionMismatchException { 110 111 final double scaledX = (normalized ? FastMath.PI : 1) * t.getValue(); 112 final double scaledX2 = scaledX * scaledX; 113 114 double[] f = new double[t.getOrder() + 1]; 115 116 if (FastMath.abs(scaledX) <= SHORTCUT) { 117 118 for (int i = 0; i < f.length; ++i) { 119 final int k = i / 2; 120 if ((i & 0x1) == 0) { 121 // even derivation order 122 f[i] = (((k & 0x1) == 0) ? 1 : -1) * 123 (1.0 / (i + 1) - scaledX2 * (1.0 / (2 * i + 6) - scaledX2 / (24 * i + 120))); 124 } else { 125 // odd derivation order 126 f[i] = (((k & 0x1) == 0) ? -scaledX : scaledX) * 127 (1.0 / (i + 2) - scaledX2 * (1.0 / (6 * i + 24) - scaledX2 / (120 * i + 720))); 128 } 129 } 130 131 } else { 132 133 final double inv = 1 / scaledX; 134 final double cos = FastMath.cos(scaledX); 135 final double sin = FastMath.sin(scaledX); 136 137 f[0] = inv * sin; 138 139 // the nth order derivative of sinc has the form: 140 // dn(sinc(x)/dxn = [S_n(x) sin(x) + C_n(x) cos(x)] / x^(n+1) 141 // where S_n(x) is an even polynomial with degree n-1 or n (depending on parity) 142 // and C_n(x) is an odd polynomial with degree n-1 or n (depending on parity) 143 // S_0(x) = 1, S_1(x) = -1, S_2(x) = -x^2 + 2, S_3(x) = 3x^2 - 6... 144 // C_0(x) = 0, C_1(x) = x, C_2(x) = -2x, C_3(x) = -x^3 + 6x... 145 // the general recurrence relations for S_n and C_n are: 146 // S_n(x) = x S_(n-1)'(x) - n S_(n-1)(x) - x C_(n-1)(x) 147 // C_n(x) = x C_(n-1)'(x) - n C_(n-1)(x) + x S_(n-1)(x) 148 // as per polynomials parity, we can store both S_n and C_n in the same array 149 final double[] sc = new double[f.length]; 150 sc[0] = 1; 151 152 double coeff = inv; 153 for (int n = 1; n < f.length; ++n) { 154 155 double s = 0; 156 double c = 0; 157 158 // update and evaluate polynomials S_n(x) and C_n(x) 159 final int kStart; 160 if ((n & 0x1) == 0) { 161 // even derivation order, S_n is degree n and C_n is degree n-1 162 sc[n] = 0; 163 kStart = n; 164 } else { 165 // odd derivation order, S_n is degree n-1 and C_n is degree n 166 sc[n] = sc[n - 1]; 167 c = sc[n]; 168 kStart = n - 1; 169 } 170 171 // in this loop, k is always even 172 for (int k = kStart; k > 1; k -= 2) { 173 174 // sine part 175 sc[k] = (k - n) * sc[k] - sc[k - 1]; 176 s = s * scaledX2 + sc[k]; 177 178 // cosine part 179 sc[k - 1] = (k - 1 - n) * sc[k - 1] + sc[k -2]; 180 c = c * scaledX2 + sc[k - 1]; 181 182 } 183 sc[0] *= -n; 184 s = s * scaledX2 + sc[0]; 185 186 coeff *= inv; 187 f[n] = coeff * (s * sin + c * scaledX * cos); 188 189 } 190 191 } 192 193 if (normalized) { 194 double scale = FastMath.PI; 195 for (int i = 1; i < f.length; ++i) { 196 f[i] *= scale; 197 scale *= FastMath.PI; 198 } 199 } 200 201 return t.compose(f); 202 203 } 204 205}