NFFT  3.3.0
kernels.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id$ */
20 
24 #include "config.h"
25 
26 #include <stdio.h>
27 #include <math.h>
28 #include <float.h>
29 #ifdef HAVE_COMPLEX_H
30 #include <complex.h>
31 #endif
32 
33 #include "kernels.h"
34 
40 C gaussian(R x, int der, const R *param) /* K(x)=EXP(-x^2/c^2) */
41 {
42  R c = param[0];
43  R value = K(0.0);
44 
45  switch (der)
46  {
47  case 0 : value = EXP(-x*x/(c*c)); break;
48  case 1 : value = -K(2.0)*x/(c*c)*EXP(-x*x/(c*c)); break;
49  case 2 : value = K(2.0)*EXP(-x*x/(c*c))*(-c*c+K(2.0)*x*x)/(c*c*c*c); break;
50  case 3 : value = -K(4.0)*x*EXP(-x*x/(c*c))*(-K(3.0)*c*c+K(2.0)*x*x)/(c*c*c*c*c*c); break;
51  case 4 : value = K(4.0)*EXP(-x*x/(c*c))*(K(3.0)*c*c*c*c-K(12.0)*c*c*x*x+K(4.0)*x*x*x*x)/(c*c*c*c*c*c*c*c); break;
52  case 5 : value = -K(8.0)*x*EXP(-x*x/(c*c))*(K(15.0)*c*c*c*c-K(20.0)*c*c*x*x+K(4.0)*x*x*x*x)/POW(c,K(10.0)); break;
53  case 6 : value = K(8.0)*EXP(-x*x/(c*c))*(-K(15.0)*c*c*c*c*c*c+K(90.0)*x*x*c*c*c*c-K(60.0)*x*x*x*x*c*c+K(8.0)*x*x*x*x*x*x)/POW(c,K(12.0)); break;
54  case 7 : value = -K(16.0)*x*EXP(-x*x/(c*c))*(-K(105.0)*c*c*c*c*c*c+K(210.0)*x*x*c*c*c*c-K(84.0)*x*x*x*x*c*c+K(8.0)*x*x*x*x*x*x)/POW(c,K(14.0)); break;
55  case 8 : value = K(16.0)*EXP(-x*x/(c*c))*(K(105.0)*c*c*c*c*c*c*c*c-K(840.0)*x*x*c*c*c*c*c*c+K(840.0)*x*x*x*x*c*c*c*c-K(224.0)*x*x*x*x*x*x*c*c+K(16.0)*x*x*x*x*x*x*x*x)/POW(c,K(16.0)); break;
56  case 9 : value = -K(32.0)*x*EXP(-x*x/(c*c))*(K(945.0)*c*c*c*c*c*c*c*c-K(2520.0)*x*x*c*c*c*c*c*c+K(1512.0)*x*x*x*x*c*c*c*c-K(288.0)*x*x*x*x*x*x*c*c+K(16.0)*x*x*x*x*x*x*x*x)/POW(c,K(18.0)); break;
57  case 10 : value = K(32.0)*EXP(-x*x/(c*c))*(-K(945.0)*POW(c,K(10.0))+K(9450.0)*x*x*c*c*c*c*c*c*c*c-K(12600.0)*x*x*x*x*c*c*c*c*c*c+K(5040.0)*x*x*x*x*x*x*c*c*c*c-K(720.0)*x*x*x*x*x*x*x*x*c*c+K(32.0)*POW(x,K(10.0)))/POW(c,K(20.0)); break;
58  case 11 : value = -K(64.0)*x*EXP(-x*x/(c*c))*(-K(10395.0)*POW(c,K(10.0))+K(34650.0)*x*x*c*c*c*c*c*c*c*c-K(27720.0)*x*x*x*x*c*c*c*c*c*c+K(7920.0)*x*x*x*x*x*x*c*c*c*c-K(880.0)*x*x*x*x*x*x*x*x*c*c+K(32.0)*POW(x,K(10.0)))/POW(c,K(22.0)); break;
59  case 12 : value = K(64.0)*EXP(-x*x/(c*c))*(K(10395.0)*POW(c,K(12.0))-K(124740.0)*x*x*POW(c,K(10.0))+K(207900.0)*x*x*x*x*c*c*c*c*c*c*c*c-K(110880.0)*x*x*x*x*x*x*c*c*c*c*c*c+K(23760.0)*x*x*x*x*x*x*x*x*c*c*c*c-K(2112.0)*POW(x,K(10.0))*c*c+K(64.0)*POW(x,K(12.0)))/POW(c,K(24.0)); break;
60  default : value = K(0.0);
61  }
62 
63  return value;
64 }
65 
66 C multiquadric(R x, int der, const R *param) /* K(x)=SQRT(x^2+c^2) */
67 {
68  R c=param[0];
69  R value=K(0.0);
70 
71  switch (der)
72  {
73  case 0 : value=SQRT(x*x+c*c); break;
74  case 1 : value=K(1.0)/(SQRT(x*x+c*c))*x; break;
75  case 2 : value=c*c/SQRT(POW(x*x+c*c,K(3.0))); break;
76  case 3 : value=-K(3.0)*x*c*c/SQRT(POW(x*x+c*c,K(5.0))); break;
77  case 4 : value=K(3.0)*c*c*(K(4.0)*x*x-c*c)/SQRT(POW(x*x+c*c,K(7.0))); break;
78  case 5 : value=-K(15.0)*x*c*c*(K(4.0)*x*x-K(3.0)*c*c)/SQRT(POW(x*x+c*c,K(9.0))); break;
79  case 6 : value=K(45.0)*c*c*(K(8.0)*x*x*x*x-K(12.0)*x*x*c*c+c*c*c*c)/SQRT(POW(x*x+c*c,K(11.0))); break;
80  case 7 : value=-K(315.0)*x*c*c*(K(8.0)*x*x*x*x-K(20.0)*x*x*c*c+K(5.0)*c*c*c*c)/SQRT(POW(x*x+c*c,K(13.0))); break;
81  case 8 : value=K(315.0)*c*c*(K(64.0)*x*x*x*x*x*x-K(240.0)*x*x*x*x*c*c+K(120.0)*x*x*c*c*c*c-K(5.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(15.0))); break;
82  case 9 : value=-K(2835.0)*x*c*c*(K(64.0)*x*x*x*x*x*x-K(336.0)*x*x*x*x*c*c+K(280.0)*x*x*c*c*c*c-K(35.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(17.0))); break;
83  case 10 : value=K(14175.0)*c*c*(K(128.0)*x*x*x*x*x*x*x*x-K(896.0)*x*x*x*x*x*x*c*c+K(1120.0)*x*x*x*x*c*c*c*c-K(280.0)*x*x*c*c*c*c*c*c+K(7.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(19.0))); break;
84  case 11 : value=-K(155925.0)*x*c*c*(K(128.0)*x*x*x*x*x*x*x*x-K(1152.0)*x*x*x*x*x*x*c*c+K(2016.0)*x*x*x*x*c*c*c*c-K(840.0)*x*x*c*c*c*c*c*c+K(63.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(21.0))); break;
85  case 12 : value=K(467775.0)*c*c*(K(1260.0)*x*x*c*c*c*c*c*c*c*c-K(21.0)*POW(c,K(10.0))+K(512.0)*POW(x,K(10.0))-K(5760.0)*x*x*x*x*x*x*x*x*c*c+K(13440.0)*x*x*x*x*x*x*c*c*c*c-K(8400.0)*x*x*x*x*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(23.0))); break;
86  default : value=K(0.0);
87  }
88 
89  return value;
90 }
91 
92 C inverse_multiquadric(R x, int der, const R *param) /* K(x)=1/SQRT(x^2+c^2) */
93 {
94  R c=param[0];
95  R value=K(0.0);
96 
97  switch (der)
98  {
99  case 0 : value=K(1.0)/SQRT(x*x+c*c); break;
100  case 1 : value=-K(1.0)/(SQRT(POW(x*x+c*c,K(3.0))))*x; break;
101  case 2 : value=(K(2.0)*x*x-c*c)/SQRT(POW(x*x+c*c,K(5.0))); break;
102  case 3 : value=-K(3.0)*x*(K(2.0)*x*x-K(3.0)*c*c)/SQRT(POW(x*x+c*c,K(7.0))); break;
103  case 4 : value=K(3.0)*(K(8.0)*x*x*x*x-K(24.0)*x*x*c*c+K(3.0)*c*c*c*c)/SQRT(POW(x*x+c*c,K(9.0))); break;
104  case 5 : value=-K(15.0)*x*(K(8.0)*x*x*x*x-K(40.0)*x*x*c*c+K(15.0)*c*c*c*c)/SQRT(POW(x*x+c*c,K(11.0))); break;
105  case 6 : value=K(45.0)*(K(16.0)*x*x*x*x*x*x-K(120.0)*x*x*x*x*c*c+K(90.0)*x*x*c*c*c*c-K(5.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(13.0))); break;
106  case 7 : value=-K(315.0)*x*(K(16.0)*x*x*x*x*x*x-K(168.0)*x*x*x*x*c*c+K(210.0)*x*x*c*c*c*c-K(35.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(15.0))); break;
107  case 8 : value=K(315.0)*(K(128.0)*x*x*x*x*x*x*x*x-K(1792.0)*x*x*x*x*x*x*c*c+K(3360.0)*x*x*x*x*c*c*c*c-K(1120.0)*x*x*c*c*c*c*c*c+K(35.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(17.0))); break;
108  case 9 : value=-K(2835.0)*x*(K(128.0)*x*x*x*x*x*x*x*x-K(2304.0)*x*x*x*x*x*x*c*c+K(6048.0)*x*x*x*x*c*c*c*c-K(3360.0)*x*x*c*c*c*c*c*c+K(315.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(19.0))); break;
109  case 10 : value=K(14175.0)*(K(256.0)*POW(x,K(10.0))-K(5760.0)*x*x*x*x*x*x*x*x*c*c+K(20160.0)*x*x*x*x*x*x*c*c*c*c-K(16800.0)*x*x*x*x*c*c*c*c*c*c+K(3150.0)*x*x*c*c*c*c*c*c*c*c-K(63.0)*POW(c,K(10.0)))/SQRT(POW(x*x+c*c,K(21.0))); break;
110  case 11 : value=-K(155925.0)*x*(K(256.0)*POW(x,K(10.0))-K(7040.0)*x*x*x*x*x*x*x*x*c*c+K(31680.0)*x*x*x*x*x*x*c*c*c*c-K(36960.0)*x*x*x*x*c*c*c*c*c*c+K(11550.0)*x*x*c*c*c*c*c*c*c*c-K(693.0)*POW(c,K(10.0)))/SQRT(POW(x*x+c*c,K(23.0))); break;
111  case 12 : value=K(467775.0)*(K(231.0)*POW(c,K(12.0))+K(190080.0)*x*x*x*x*x*x*x*x*c*c*c*c-K(16632.0)*x*x*POW(c,K(10.0))-K(295680.0)*x*x*x*x*x*x*c*c*c*c*c*c+K(138600.0)*x*x*x*x*c*c*c*c*c*c*c*c+K(1024.0)*POW(x,K(12.0))-K(33792.0)*POW(x,K(10.0))*c*c)/SQRT(POW(x*x+c*c,K(25.0))); break;
112  default : value=K(0.0);
113  }
114 
115  return value;
116 }
117 
118 C logarithm(R x, int der, const R *param) /* K(x)=LOG |x| */
119 {
120  R value=K(0.0);
121 
122  (void)param;
123 
124  if (FABS(x)<DBL_EPSILON) value=K(0.0);
125  else switch (der)
126  {
127  case 0 : value=LOG(FABS(x)); break;
128  case 1 : value=(x<0 ? -1 : 1)/FABS(x); break;
129  case 2 : value=-1/(x*x); break;
130  case 3 : value=K(2.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(3.0)); break;
131  case 4 : value=-K(6.0)/(x*x*x*x); break;
132  case 5 : value=K(24.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(5.0)); break;
133  case 6 : value=-K(120.0)/(x*x*x*x*x*x); break;
134  case 7 : value=K(720.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(7.0)); break;
135  case 8 : value=-K(5040.0)/(x*x*x*x*x*x*x*x); break;
136  case 9 : value=K(40320.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(9.0)); break;
137  case 10 : value=-K(362880.0)/POW(x,K(10.0)); break;
138  case 11 : value=K(3628800.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(11.0)); break;
139  case 12 : value=-K(39916800.0)/POW(x,K(12.0)); break;
140  case 13 : value=K(479001600.0)/POW(x,K(13.0)); break;
141  case 14 : value=-K(6227020800.0)/POW(x,K(14.0)); break;
142  case 15 : value=K(87178291200.0)/POW(x,K(15.0)); break;
143  case 16 : value=-K(1307674368000.0)/POW(x,K(16.0)); break;
144  case 17 : value=K(20922789888000.0)/POW(x,K(17.0)); break;
145  default : value=K(0.0);
146  }
147 
148  return value;
149 }
150 
151 C thinplate_spline(R x, int der, const R *param) /* K(x) = x^2 LOG |x| */
152 {
153  R value=K(0.0);
154 
155  (void)param;
156 
157  if (FABS(x)<DBL_EPSILON) value=K(0.0);
158  else switch (der)
159  {
160  case 0 : value=x*x*LOG(FABS(x)); break;
161  case 1 : value=K(2.0)*x*LOG(FABS(x))+x; break;
162  case 2 : value=K(2.0)*LOG(FABS(x))+K(3.0); break;
163  case 3 : value=K(2.0)/x; break;
164  case 4 : value=-K(2.0)/(x*x); break;
165  case 5 : value=K(4.0)/(x*x*x); break;
166  case 6 : value=-K(12.0)/(x*x*x*x); break;
167  case 7 : value=K(48.0)/(x*x*x*x*x); break;
168  case 8 : value=-K(240.0)/(x*x*x*x*x*x); break;
169  case 9 : value=K(1440.0)/(x*x*x*x*x*x*x); break;
170  case 10 : value=-K(10080.0)/(x*x*x*x*x*x*x*x); break;
171  case 11 : value=K(80640.0)/(x*x*x*x*x*x*x*x*x); break;
172  case 12 : value=-K(725760.0)/POW(x,K(10.0)); break;
173  default : value=K(0.0);
174  }
175 
176  return value;
177 }
178 
179 C one_over_square(R x, int der, const R *param) /* K(x) = 1/x^2 */
180 {
181  R value=K(0.0);
182 
183  (void)param;
184 
185  if (FABS(x)<DBL_EPSILON) value=K(0.0);
186  else switch (der)
187  {
188  case 0 : value=K(1.0)/(x*x); break;
189  case 1 : value=-K(2.0)/(x*x*x); break;
190  case 2 : value=K(6.0)/(x*x*x*x); break;
191  case 3 : value=-K(24.0)/(x*x*x*x*x); break;
192  case 4 : value=K(120.0)/(x*x*x*x*x*x); break;
193  case 5 : value=-K(720.0)/(x*x*x*x*x*x*x); break;
194  case 6 : value=K(5040.0)/(x*x*x*x*x*x*x*x); break;
195  case 7 : value=-K(40320.0)/(x*x*x*x*x*x*x*x*x); break;
196  case 8 : value=K(362880.0)/POW(x,K(10.0)); break;
197  case 9 : value=-K(3628800.0)/POW(x,K(11.0)); break;
198  case 10 : value=K(39916800.0)/POW(x,K(12.0)); break;
199  case 11 : value=-K(479001600.0)/POW(x,K(13.0)); break;
200  case 12 : value=K(6227020800.0)/POW(x,K(14.0)); break;
201  default : value=K(0.0);
202  }
203 
204  return value;
205 }
206 
207 C one_over_modulus(R x, int der, const R *param) /* K(x) = 1/|x| */
208 {
209  R value=K(0.0);
210 
211  (void)param;
212 
213  if (FABS(x)<DBL_EPSILON) value=K(0.0);
214  else switch (der)
215  {
216  case 0 : value=K(1.0)/FABS(x); break;
217  case 1 : value=-1/x/FABS(x); break;
218  case 2 : value=K(2.0)/POW(FABS(x),K(3.0)); break;
219  case 3 : value=-K(6.0)/(x*x*x)/FABS(x); break;
220  case 4 : value=K(24.0)/POW(FABS(x),K(5.0)); break;
221  case 5 : value=-K(120.0)/(x*x*x*x*x)/FABS(x); break;
222  case 6 : value=K(720.0)/POW(FABS(x),K(7.0)); break;
223  case 7 : value=-K(5040.0)/(x*x*x*x*x*x*x)/FABS(x); break;
224  case 8 : value=K(40320.0)/POW(FABS(x),K(9.0)); break;
225  case 9 : value=-K(362880.0)/(x*x*x*x*x*x*x*x*x)/FABS(x); break;
226  case 10 : value=K(3628800.0)/POW(FABS(x),K(11.0)); break;
227  case 11 : value=-K(39916800.0)/POW(x,K(11.0))/FABS(x); break;
228  case 12 : value=K(479001600.0)/POW(FABS(x),K(13.0)); break;
229  default : value=K(0.0);
230  }
231 
232  return value;
233 }
234 
235 C one_over_x(R x, int der, const R *param) /* K(x) = 1/x */
236 {
237  R value=K(0.0);
238 
239  (void)param;
240 
241  if (FABS(x)<DBL_EPSILON) value=K(0.0);
242  else switch (der)
243  {
244  case 0 : value=K(1.0)/x; break;
245  case 1 : value=-K(1.0)/(x*x); break;
246  case 2 : value=K(2.0)/(x*x*x); break;
247  case 3 : value=-K(6.0)/(x*x*x*x); break;
248  case 4 : value=K(24.0)/(x*x*x*x*x); break;
249  case 5 : value=-K(120.0)/(x*x*x*x*x*x); break;
250  case 6 : value=K(720.0)/(x*x*x*x*x*x*x); break;
251  case 7 : value=-K(5040.0)/(x*x*x*x*x*x*x*x); break;
252  case 8 : value=K(40320.0)/(x*x*x*x*x*x*x*x*x); break;
253  case 9 : value=-K(362880.0)/POW(x,K(10.0)); break;
254  case 10 : value=K(3628800.0)/POW(x,K(11.0)); break;
255  case 11 : value=-K(39916800.0)/POW(x,K(12.0)); break;
256  case 12 : value=K(479001600.0)/POW(x,K(13.0)); break;
257  default : value=K(0.0);
258  }
259 
260  return value;
261 }
262 
263 C inverse_multiquadric3(R x, int der, const R *param) /* K(x) = 1/SQRT(x^2+c^2)^3 */
264 {
265  R c=param[0];
266  R value=K(0.0);
267 
268  switch (der)
269  {
270  case 0 : value=K(1.0)/(SQRT(POW(x*x+c*c,K(3.0)))); break;
271  case 1 : value=-K(3.0)/SQRT(POW(x*x+c*c,K(5.0)))*x; break;
272  case 2 : value=K(3.0)*(K(4.0)*x*x-c*c)/SQRT(POW(x*x+c*c,K(7.0))); break;
273  case 3 : value=-K(15.0)*x*(K(4.0)*x*x-K(3.0)*c*c)/SQRT(POW(x*x+c*c,K(9.0))); break;
274  case 4 : value=K(45.0)*(K(8.0)*x*x*x*x-K(12.0)*x*x*c*c+c*c*c*c)/SQRT(POW(x*x+c*c,K(11.0))); break;
275  case 5 : value=-K(315.0)*x*(K(8.0)*x*x*x*x-K(20.0)*x*x*c*c+K(5.0)*c*c*c*c)/SQRT(POW(x*x+c*c,K(13.0))); break;
276  case 6 : value=K(315.0)*(K(64.0)*x*x*x*x*x*x-K(240.0)*x*x*x*x*c*c+K(120.0)*x*x*c*c*c*c-K(5.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(15.0))); break;
277  case 7 : value=-K(2835.0)*x*(K(64.0)*x*x*x*x*x*x-K(336.0)*x*x*x*x*c*c+K(280.0)*x*x*c*c*c*c-K(35.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(17.0))); break;
278  case 8 : value=K(14175.0)*(K(128.0)*x*x*x*x*x*x*x*x-K(896.0)*x*x*x*x*x*x*c*c+K(1120.0)*x*x*x*x*c*c*c*c-K(280.0)*x*x*c*c*c*c*c*c+K(7.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(19.0))); break;
279  case 9 : value=-K(155925.0)*x*(K(128.0)*x*x*x*x*x*x*x*x-K(1152.0)*x*x*x*x*x*x*c*c+K(2016.0)*x*x*x*x*c*c*c*c-K(840.0)*x*x*c*c*c*c*c*c+K(63.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(21.0))); break;
280  case 10 : value=K(467775.0)*(K(512.0)*POW(x,K(10.0))-K(5760.0)*x*x*x*x*x*x*x*x*c*c+K(13440.0)*x*x*x*x*x*x*c*c*c*c-K(8400.0)*x*x*x*x*c*c*c*c*c*c+K(1260.0)*x*x*c*c*c*c*c*c*c*c-K(21.0)*POW(c,K(10.0)))/SQRT(POW(x*x+c*c,K(23.0))); break;
281  case 11 : value=-K(6081075.0)*x*(K(512.0)*POW(x,K(10.0))-K(7040.0)*x*x*x*x*x*x*x*x*c*c+K(21120.0)*x*x*x*x*x*x*c*c*c*c-K(18480.0)*x*x*x*x*c*c*c*c*c*c+K(4620.0)*x*x*c*c*c*c*c*c*c*c-K(231.0)*POW(c,K(10.0)))/SQRT(POW(x*x+c*c,K(25.0))); break;
282  case 12 : value=K(42567525.0)*(K(1024.0)*POW(x,K(12.0))+K(27720.0)*x*x*x*x*c*c*c*c*c*c*c*c+K(33.0)*POW(c,K(12.0))-K(2772.0)*x*x*POW(c,K(10.0))-K(73920.0)*x*x*x*x*x*x*c*c*c*c*c*c+K(63360.0)*x*x*x*x*x*x*x*x*c*c*c*c-K(16896.0)*POW(x,K(10.0))*c*c)/SQRT(POW(x*x+c*c,K(27.0))); break;
283  default : value=K(0.0);
284  }
285 
286  return value;
287 }
288 
289 C sinc_kernel(R x, int der, const R *param) /* K(x) = SIN(cx)/x */
290 {
291  R c=param[0];
292  R value=K(0.0);
293 
294  if (FABS(x)<DBL_EPSILON) value=K(0.0);
295  else switch (der)
296  {
297  case 0 : value=SIN(c*x)/x; break;
298  case 1 : value=(COS(c*x)*c*x-SIN(c*x))/(x*x); break;
299  case 2 : value=-(SIN(c*x)*c*c*x*x+K(2.0)*COS(c*x)*c*x-K(2.0)*SIN(c*x))/(x*x*x); break;
300  case 3 : value=-(COS(c*x)*c*c*c*x*x*x-K(3.0)*SIN(c*x)*c*c*x*x-K(6.0)*COS(c*x)*c*x+K(6.0)*SIN(c*x))/(x*x*x*x); break;
301  case 4 : value=(SIN(c*x)*c*c*c*c*x*x*x*x+K(4.0)*COS(c*x)*c*c*c*x*x*x-K(12.0)*SIN(c*x)*c*c*x*x-K(24.0)*COS(c*x)*c*x+K(24.0)*SIN(c*x))/(x*x*x*x*x); break;
302  case 5 : value=(COS(c*x)*c*c*c*c*c*x*x*x*x*x-K(5.0)*SIN(c*x)*c*c*c*c*x*x*x*x-K(20.0)*COS(c*x)*c*c*c*x*x*x+K(60.0)*SIN(c*x)*c*c*x*x+K(120.0)*COS(c*x)*c*x-K(120.0)*SIN(c*x))/(x*x*x*x*x*x); break;
303  case 6 : value=-(SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(6.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x-K(30.0)*SIN(c*x)*c*c*c*c*x*x*x*x-K(120.0)*COS(c*x)*c*c*c*x*x*x+K(360.0)*SIN(c*x)*c*c*x*x+K(720.0)*COS(c*x)*c*x-K(720.0)*SIN(c*x))/(x*x*x*x*x*x*x); break;
304  case 7 : value=-(COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(7.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(42.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x+K(210.0)*SIN(c*x)*c*c*c*c*x*x*x*x+K(840.0)*COS(c*x)*c*c*c*x*x*x-K(2520.0)*SIN(c*x)*c*c*x*x-K(5040.0)*COS(c*x)*c*x+K(5040.0)*SIN(c*x))/(x*x*x*x*x*x*x*x); break;
305  case 8 : value=(SIN(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+K(8.0)*COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(56.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(336.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x+K(1680.0)*SIN(c*x)*c*c*c*c*x*x*x*x+K(6720.0)*COS(c*x)*c*c*c*x*x*x-K(20160.0)*SIN(c*x)*c*c*x*x-K(40320.0)*COS(c*x)*c*x+K(40320.0)*SIN(c*x))/(x*x*x*x*x*x*x*x*x); break;
306  case 9 : value=(COS(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-K(9.0)*SIN(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-K(72.0)*COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+K(504.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(3024.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x-K(15120.0)*SIN(c*x)*c*c*c*c*x*x*x*x-K(60480.0)*COS(c*x)*c*c*c*x*x*x+K(181440.0)*SIN(c*x)*c*c*x*x+K(362880.0)*COS(c*x)*c*x-K(362880.0)*SIN(c*x))/POW(x,K(10.0)); break;
307  case 10 : value=-(SIN(c*x)*POW(c,K(10.0))*POW(x,K(10.0))+K(10.0)*COS(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-K(90.0)*SIN(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-K(720.0)*COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+K(5040.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(30240.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x-K(151200.0)*SIN(c*x)*c*c*c*c*x*x*x*x-K(604800.0)*COS(c*x)*c*c*c*x*x*x+K(1814400.0)*SIN(c*x)*c*c*x*x+K(3628800.0)*COS(c*x)*c*x-K(3628800.0)*SIN(c*x))/POW(x,K(11.0)); break;
308  case 11 : value=-(COS(c*x)*POW(c,K(11.0))*POW(x,K(11.0))-K(11.0)*SIN(c*x)*POW(c,K(10.0))*POW(x,K(10.0))-K(110.0)*COS(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+K(990.0)*SIN(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+K(7920.0)*COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(55440.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(332640.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x+K(1663200.0)*SIN(c*x)*c*c*c*c*x*x*x*x+K(6652800.0)*COS(c*x)*c*c*c*x*x*x-K(19958400.0)*SIN(c*x)*c*c*x*x-K(39916800.0)*COS(c*x)*c*x+K(39916800.0)*SIN(c*x))/POW(x,K(12.0)); break;
309  case 12 : value=(SIN(c*x)*POW(c,K(12.0))*POW(x,K(12.0))+K(12.0)*COS(c*x)*POW(c,K(11.0))*POW(x,K(11.0))-K(132.0)*SIN(c*x)*POW(c,K(10.0))*POW(x,K(10.0))-K(1320.0)*COS(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+K(11880.0)*SIN(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+K(95040.0)*COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(665280.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(3991680.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x+K(19958400.0)*SIN(c*x)*c*c*c*c*x*x*x*x+K(79833600.0)*COS(c*x)*c*c*c*x*x*x-K(239500800.0)*SIN(c*x)*c*c*x*x-K(479001600.0)*COS(c*x)*c*x+K(479001600.0)*SIN(c*x))/POW(x,K(13.0)); break;
310  default : value=K(0.0);
311  }
312 
313  return value;
314 }
315 
316 C cosc(R x, int der, const R *param) /* K(x) = COS(cx)/x */
317 {
318  R c=param[0];
319  R value=K(0.0);
320  R sign;
321 
322  if (x<0) sign=-K(1.0); else sign=K(1.0);
323  x=FABS(x);
324 
325  if (FABS(x)<DBL_EPSILON) value=K(0.0);
326  else switch (der)
327  {
328  case 0 : value=COS(c*x)/x; break;
329  case 1 : value=-(SIN(c*x)*c*x+COS(c*x))/(x*x); break;
330  case 2 : value=(-COS(c*x)*c*c*x*x+K(2.0)*SIN(c*x)*c*x+K(2.0)*COS(c*x))/(x*x*x); break;
331  case 3 : value=(SIN(c*x)*c*c*c*x*x*x+K(3.0)*COS(c*x)*c*c*x*x-K(6.0)*SIN(c*x)*c*x-K(6.0)*COS(c*x))/(x*x*x*x); break;
332  case 4 : value=(COS(c*x)*c*c*c*c*x*x*x*x-K(4.0)*SIN(c*x)*c*c*c*x*x*x-K(12.0)*COS(c*x)*c*c*x*x+K(24.0)*SIN(c*x)*c*x+K(24.0)*COS(c*x))/(x*x*x*x*x); break;
333  case 5 : value=-(SIN(c*x)*c*c*c*c*c*x*x*x*x*x+K(5.0)*COS(c*x)*c*c*c*c*x*x*x*x-K(20.0)*SIN(c*x)*c*c*c*x*x*x-K(60.0)*COS(c*x)*c*c*x*x+K(120.0)*SIN(c*x)*c*x+K(120.0)*COS(c*x))/(x*x*x*x*x*x); break;
334  case 6 : value=-(COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(6.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x-K(30.0)*COS(c*x)*c*c*c*c*x*x*x*x+K(120.0)*SIN(c*x)*c*c*c*x*x*x+K(360.0)*COS(c*x)*c*c*x*x-K(720.0)*SIN(c*x)*c*x-K(720.0)*COS(c*x))/(x*x*x*x*x*x*x); break;
335  case 7 : value=(SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+K(7.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(42.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x-K(210.0)*COS(c*x)*c*c*c*c*x*x*x*x+K(840.0)*SIN(c*x)*c*c*c*x*x*x+K(2520.0)*COS(c*x)*c*c*x*x-K(5040.0)*SIN(c*x)*c*x-K(5040.0)*COS(c*x))/(x*x*x*x*x*x*x*x); break;
336  case 8 : value=(COS(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-K(8.0)*SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(56.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(336.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x+K(1680.0)*COS(c*x)*c*c*c*c*x*x*x*x-K(6720.0)*SIN(c*x)*c*c*c*x*x*x-K(20160.0)*COS(c*x)*c*c*x*x+K(40320.0)*SIN(c*x)*c*x+K(40320.0)*COS(c*x))/(x*x*x*x*x*x*x*x*x); break;
337  case 9 : value=-(SIN(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+K(9.0)*COS(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-K(72.0)*SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(504.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(3024.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x+K(15120.0)*COS(c*x)*c*c*c*c*x*x*x*x-K(60480.0)*SIN(c*x)*c*c*c*x*x*x-K(181440.0)*COS(c*x)*c*c*x*x+K(362880.0)*SIN(c*x)*c*x+K(362880.0)*COS(c*x))/POW(x,K(10.0)); break;
338  case 10 : value=-(COS(c*x)*POW(c,K(10.0))*POW(x,K(10.0))-K(10.0)*SIN(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-K(90.0)*COS(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+K(720.0)*SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+K(5040.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(30240.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x-K(151200.0)*COS(c*x)*c*c*c*c*x*x*x*x+K(604800.0)*SIN(c*x)*c*c*c*x*x*x+K(1814400.0)*COS(c*x)*c*c*x*x-K(3628800.0)*SIN(c*x)*c*x-K(3628800.0)*COS(c*x))/POW(x,K(11.0)); break;
339  case 11 : value=(SIN(c*x)*POW(c,K(11.0))*POW(x,K(11.0))+K(11.0)*COS(c*x)*POW(c,K(10.0))*POW(x,K(10.0))-K(110.0)*SIN(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-K(990.0)*COS(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+K(7920.0)*SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+K(55440.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(332640.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x-K(1663200.0)*COS(c*x)*c*c*c*c*x*x*x*x+K(6652800.0)*SIN(c*x)*c*c*c*x*x*x+K(19958400.0)*COS(c*x)*c*c*x*x-K(39916800.0)*SIN(c*x)*c*x-K(39916800.0)*COS(c*x))/POW(x,K(12.0)); break;
340  case 12 : value=(COS(c*x)*POW(c,K(12.0))*POW(x,K(12.0))-K(12.0)*SIN(c*x)*POW(c,K(11.0))*POW(x,K(11.0))-K(132.0)*COS(c*x)*POW(c,K(10.0))*POW(x,K(10.0))+K(1320.0)*SIN(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+K(11880.0)*COS(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-K(95040.0)*SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(665280.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(3991680.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x+K(19958400.0)*COS(c*x)*c*c*c*c*x*x*x*x-K(79833600.0)*SIN(c*x)*c*c*c*x*x*x-K(239500800.0)*COS(c*x)*c*c*x*x+K(479001600.0)*SIN(c*x)*c*x+K(479001600.0)*COS(c*x))/POW(x,K(13.0)); break;
341  default : value=K(0.0);
342  }
343  value *= POW(sign, (R)(der));
344 
345  return value;
346 }
347 
348 C kcot(R x, int der, const R *param) /* K(x) = cot(cx) */
349 {
350  R c=param[0];
351  R value=K(0.0);
352 
353  if (FABS(x)<DBL_EPSILON) value=K(0.0);
354  else switch (der)
355  {
356  case 0 : value = K(1.0)/TAN(c * x); break;
357  case 1 : value = -(K(1.0) + POW(K(1.0)/TAN(c * x), K(2.0))) * c; break;
358  case 2 : value = K(2.0) * K(1.0)/TAN(c * x) * (K(1.0) + POW(K(1.0)/TAN(c * x), K(2.0))) * c * c; break;
359  case 3 : value = -K(2.0) * (K(1.0) + POW(K(1.0)/TAN(c * x), K(2.0))) * POW(c, K(3.0)) * (K(1.0) + K(3.0) * POW(K(1.0)/TAN(c * x), K(2.0))); break;
360  case 4 : value = K(8.0) * (K(1.0) + POW(K(1.0)/TAN(c * x), K(2.0))) * POW(c, K(4.0)) * K(1.0)/TAN(c * x) * (K(2.0) + K(3.0) * POW(K(1.0)/TAN(c * x), K(2.0))); break;
361  case 5 : value = -K(0.8e1) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.5e1)) * (K(0.15e2) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.15e2) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.2e1)); break;
362  case 6 : value = K(0.16e2) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.6e1)) * K(1.0)/TAN(c * x) * (K(0.60e2) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.45e2) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.17e2)); break;
363  case 7 : value = -K(0.16e2) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.7e1)) * (K(0.525e3) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.315e3) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.231e3) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.17e2)); break;
364  case 8 : value = K(0.128e3) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.8e1)) * K(1.0)/TAN(c * x) * (K(0.630e3) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.315e3) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.378e3) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.62e2)); break;
365  case 9 : value = -K(0.128e3) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.9e1)) * (K(0.6615e4) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.2835e4) * POW(K(1.0)/TAN(c * x), K(0.8e1)) + K(0.5040e4) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.1320e4) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.62e2)); break;
366  case 10 : value = K(0.256e3) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.10e2)) * K(1.0)/TAN(c * x) * (K(0.37800e5) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.14175e5) * POW(K(1.0)/TAN(c * x), K(0.8e1)) + K(0.34965e5) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.12720e5) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.1382e4)); break;
367  case 11 : value = -K(0.256e3) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.11e2)) * (K(0.467775e6) * POW(K(1.0)/TAN(c * x), K(0.8e1)) + K(0.155925e6) * POW(K(1.0)/TAN(c * x), K(0.10e2)) + K(0.509355e6) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.238425e6) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.42306e5) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.1382e4)); break;
368  case 12 : value = K(0.1024e4) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.12e2)) * K(1.0)/TAN(c * x) * (K(0.1559250e7) * POW(K(1.0)/TAN(c * x), K(0.8e1)) + K(0.467775e6) * POW(K(1.0)/TAN(c * x), K(0.10e2)) + K(0.1954260e7) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.1121670e7) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.280731e6) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.21844e5)); break;
369  default : value=K(0.0);
370  }
371 
372  return value;
373 }
374 
375 
376 C one_over_cube(R x, int der, const R *param)
377 {
378  R value=K(0.0);
379  UNUSED(param);
380 
381  if (FABS(x)<DBL_EPSILON) value=K(0.0);
382  else switch (der)
383  {
384  case 0 : value = K(1.0)/(x*x*x); break;
385  case 1 : value = -K(3.0)/(x*x*x*x); break;
386  case 2 : value = K(12.0)/(x*x*x*x*x); break;
387  case 3 : value = -K(60.0)/(x*x*x*x*x*x); break;
388  case 4 : value = K(360.0)/(x*x*x*x*x*x*x); break;
389  case 5 : value = -K(2520.0)/(x*x*x*x*x*x*x*x); break;
390  case 6 : value = K(20160.0)/POW(x, K(9.0)); break;
391  case 7 : value = -K(181440.0)/POW(x, K(10.0)); break;
392  case 8 : value = K(1814400.0)/POW(x, K(11.0)); break;
393  case 9 : value = -K(19958400.0)/POW(x, K(12.0)); break;
394  case 10 : value = K(239500800.0)/POW(x, K(13.0)); break;
395  case 11 : value = -K(3113510400.0)/POW(x, K(14.0)); break;
396  case 12 : value = K(43589145600.0)/POW(x, K(15.0)); break;
397  default : value=K(0.0);
398  }
399 
400  return value;
401 }
402 
403 
404 /* \} */
405 
406 /* kernels.c */
Header file with predefined kernels for the fast summation algorithm.
#define UNUSED(x)
Dummy use of unused parameters to silence compiler warnings.
Definition: infft.h:1383