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.legacy.analysis.function; 019 020import org.apache.commons.math4.legacy.analysis.differentiation.DerivativeStructure; 021import org.apache.commons.math4.legacy.analysis.differentiation.UnivariateDifferentiableFunction; 022import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 023import org.apache.commons.math4.core.jdkmath.JdkMath; 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(πx) / π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 ? JdkMath.PI * x : x; 085 if (JdkMath.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 JdkMath.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 ? JdkMath.PI : 1) * t.getValue(); 103 final double scaledX2 = scaledX * scaledX; 104 105 double[] f = new double[t.getOrder() + 1]; 106 107 if (JdkMath.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 } else { 122 123 final double inv = 1 / scaledX; 124 final double cos = JdkMath.cos(scaledX); 125 final double sin = JdkMath.sin(scaledX); 126 127 f[0] = inv * sin; 128 129 // the nth order derivative of sinc has the form: 130 // dn(sinc(x)/dxn = [S_n(x) sin(x) + C_n(x) cos(x)] / x^(n+1) 131 // where S_n(x) is an even polynomial with degree n-1 or n (depending on parity) 132 // and C_n(x) is an odd polynomial with degree n-1 or n (depending on parity) 133 // S_0(x) = 1, S_1(x) = -1, S_2(x) = -x^2 + 2, S_3(x) = 3x^2 - 6... 134 // C_0(x) = 0, C_1(x) = x, C_2(x) = -2x, C_3(x) = -x^3 + 6x... 135 // the general recurrence relations for S_n and C_n are: 136 // S_n(x) = x S_(n-1)'(x) - n S_(n-1)(x) - x C_(n-1)(x) 137 // C_n(x) = x C_(n-1)'(x) - n C_(n-1)(x) + x S_(n-1)(x) 138 // as per polynomials parity, we can store both S_n and C_n in the same array 139 final double[] sc = new double[f.length]; 140 sc[0] = 1; 141 142 double coeff = inv; 143 for (int n = 1; n < f.length; ++n) { 144 145 double s = 0; 146 double c = 0; 147 148 // update and evaluate polynomials S_n(x) and C_n(x) 149 final int kStart; 150 if ((n & 0x1) == 0) { 151 // even derivation order, S_n is degree n and C_n is degree n-1 152 sc[n] = 0; 153 kStart = n; 154 } else { 155 // odd derivation order, S_n is degree n-1 and C_n is degree n 156 sc[n] = sc[n - 1]; 157 c = sc[n]; 158 kStart = n - 1; 159 } 160 161 // in this loop, k is always even 162 for (int k = kStart; k > 1; k -= 2) { 163 164 // sine part 165 sc[k] = (k - n) * sc[k] - sc[k - 1]; 166 s = s * scaledX2 + sc[k]; 167 168 // cosine part 169 sc[k - 1] = (k - 1 - n) * sc[k - 1] + sc[k -2]; 170 c = c * scaledX2 + sc[k - 1]; 171 } 172 sc[0] *= -n; 173 s = s * scaledX2 + sc[0]; 174 175 coeff *= inv; 176 f[n] = coeff * (s * sin + c * scaledX * cos); 177 } 178 } 179 180 if (normalized) { 181 double scale = JdkMath.PI; 182 for (int i = 1; i < f.length; ++i) { 183 f[i] *= scale; 184 scale *= JdkMath.PI; 185 } 186 } 187 188 return t.compose(f); 189 } 190}