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 }