NFFT  3.4.1
kernels.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2017 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 
22 #include "config.h"
23 
24 #include <stdio.h>
25 #include <math.h>
26 #include <float.h>
27 #ifdef HAVE_COMPLEX_H
28 #include <complex.h>
29 #endif
30 
31 #include "kernels.h"
32 
38 C gaussian(R x, int der, const R *param) /* K(x)=EXP(-x^2/c^2) */
39 {
40  R c = param[0];
41  R value = K(0.0);
42 
43  switch (der)
44  {
45  case 0 : value = EXP(-x*x/(c*c)); break;
46  case 1 : value = -K(2.0)*x/(c*c)*EXP(-x*x/(c*c)); break;
47  case 2 : value = K(2.0)*EXP(-x*x/(c*c))*(-c*c+K(2.0)*x*x)/(c*c*c*c); break;
48  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;
49  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;
50  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;
51  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;
52  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;
53  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;
54  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;
55  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;
56  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;
57  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;
58  default : value = K(0.0);
59  }
60 
61  return value;
62 }
63 
64 C multiquadric(R x, int der, const R *param) /* K(x)=SQRT(x^2+c^2) */
65 {
66  R c=param[0];
67  R value=K(0.0);
68 
69  switch (der)
70  {
71  case 0 : value=SQRT(x*x+c*c); break;
72  case 1 : value=K(1.0)/(SQRT(x*x+c*c))*x; break;
73  case 2 : value=c*c/SQRT(POW(x*x+c*c,K(3.0))); break;
74  case 3 : value=-K(3.0)*x*c*c/SQRT(POW(x*x+c*c,K(5.0))); break;
75  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;
76  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;
77  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;
78  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;
79  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;
80  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;
81  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;
82  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;
83  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;
84  default : value=K(0.0);
85  }
86 
87  return value;
88 }
89 
90 C inverse_multiquadric(R x, int der, const R *param) /* K(x)=1/SQRT(x^2+c^2) */
91 {
92  R c=param[0];
93  R value=K(0.0);
94 
95  switch (der)
96  {
97  case 0 : value=K(1.0)/SQRT(x*x+c*c); break;
98  case 1 : value=-K(1.0)/(SQRT(POW(x*x+c*c,K(3.0))))*x; break;
99  case 2 : value=(K(2.0)*x*x-c*c)/SQRT(POW(x*x+c*c,K(5.0))); break;
100  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;
101  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;
102  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;
103  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;
104  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;
105  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;
106  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;
107  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;
108  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;
109  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;
110  default : value=K(0.0);
111  }
112 
113  return value;
114 }
115 
116 C logarithm(R x, int der, const R *param) /* K(x)=LOG |x| */
117 {
118  R value=K(0.0);
119 
120  (void)param;
121 
122  if (FABS(x)<DBL_EPSILON) value=K(0.0);
123  else switch (der)
124  {
125  case 0 : value=LOG(FABS(x)); break;
126  case 1 : value=(x<0 ? -1 : 1)/FABS(x); break;
127  case 2 : value=-1/(x*x); break;
128  case 3 : value=K(2.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(3.0)); break;
129  case 4 : value=-K(6.0)/(x*x*x*x); break;
130  case 5 : value=K(24.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(5.0)); break;
131  case 6 : value=-K(120.0)/(x*x*x*x*x*x); break;
132  case 7 : value=K(720.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(7.0)); break;
133  case 8 : value=-K(5040.0)/(x*x*x*x*x*x*x*x); break;
134  case 9 : value=K(40320.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(9.0)); break;
135  case 10 : value=-K(362880.0)/POW(x,K(10.0)); break;
136  case 11 : value=K(3628800.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(11.0)); break;
137  case 12 : value=-K(39916800.0)/POW(x,K(12.0)); break;
138  case 13 : value=K(479001600.0)/POW(x,K(13.0)); break;
139  case 14 : value=-K(6227020800.0)/POW(x,K(14.0)); break;
140  case 15 : value=K(87178291200.0)/POW(x,K(15.0)); break;
141  case 16 : value=-K(1307674368000.0)/POW(x,K(16.0)); break;
142  case 17 : value=K(20922789888000.0)/POW(x,K(17.0)); break;
143  default : value=K(0.0);
144  }
145 
146  return value;
147 }
148 
149 C thinplate_spline(R x, int der, const R *param) /* K(x) = x^2 LOG |x| */
150 {
151  R value=K(0.0);
152 
153  (void)param;
154 
155  if (FABS(x)<DBL_EPSILON) value=K(0.0);
156  else switch (der)
157  {
158  case 0 : value=x*x*LOG(FABS(x)); break;
159  case 1 : value=K(2.0)*x*LOG(FABS(x))+x; break;
160  case 2 : value=K(2.0)*LOG(FABS(x))+K(3.0); break;
161  case 3 : value=K(2.0)/x; break;
162  case 4 : value=-K(2.0)/(x*x); break;
163  case 5 : value=K(4.0)/(x*x*x); break;
164  case 6 : value=-K(12.0)/(x*x*x*x); break;
165  case 7 : value=K(48.0)/(x*x*x*x*x); break;
166  case 8 : value=-K(240.0)/(x*x*x*x*x*x); break;
167  case 9 : value=K(1440.0)/(x*x*x*x*x*x*x); break;
168  case 10 : value=-K(10080.0)/(x*x*x*x*x*x*x*x); break;
169  case 11 : value=K(80640.0)/(x*x*x*x*x*x*x*x*x); break;
170  case 12 : value=-K(725760.0)/POW(x,K(10.0)); break;
171  default : value=K(0.0);
172  }
173 
174  return value;
175 }
176 
177 C one_over_square(R x, int der, const R *param) /* K(x) = 1/x^2 */
178 {
179  R value=K(0.0);
180 
181  (void)param;
182 
183  if (FABS(x)<DBL_EPSILON) value=K(0.0);
184  else switch (der)
185  {
186  case 0 : value=K(1.0)/(x*x); break;
187  case 1 : value=-K(2.0)/(x*x*x); break;
188  case 2 : value=K(6.0)/(x*x*x*x); break;
189  case 3 : value=-K(24.0)/(x*x*x*x*x); break;
190  case 4 : value=K(120.0)/(x*x*x*x*x*x); break;
191  case 5 : value=-K(720.0)/(x*x*x*x*x*x*x); break;
192  case 6 : value=K(5040.0)/(x*x*x*x*x*x*x*x); break;
193  case 7 : value=-K(40320.0)/(x*x*x*x*x*x*x*x*x); break;
194  case 8 : value=K(362880.0)/POW(x,K(10.0)); break;
195  case 9 : value=-K(3628800.0)/POW(x,K(11.0)); break;
196  case 10 : value=K(39916800.0)/POW(x,K(12.0)); break;
197  case 11 : value=-K(479001600.0)/POW(x,K(13.0)); break;
198  case 12 : value=K(6227020800.0)/POW(x,K(14.0)); break;
199  default : value=K(0.0);
200  }
201 
202  return value;
203 }
204 
205 C one_over_modulus(R x, int der, const R *param) /* K(x) = 1/|x| */
206 {
207  R value=K(0.0);
208 
209  (void)param;
210 
211  if (FABS(x)<DBL_EPSILON) value=K(0.0);
212  else switch (der)
213  {
214  case 0 : value=K(1.0)/FABS(x); break;
215  case 1 : value=-1/x/FABS(x); break;
216  case 2 : value=K(2.0)/POW(FABS(x),K(3.0)); break;
217  case 3 : value=-K(6.0)/(x*x*x)/FABS(x); break;
218  case 4 : value=K(24.0)/POW(FABS(x),K(5.0)); break;
219  case 5 : value=-K(120.0)/(x*x*x*x*x)/FABS(x); break;
220  case 6 : value=K(720.0)/POW(FABS(x),K(7.0)); break;
221  case 7 : value=-K(5040.0)/(x*x*x*x*x*x*x)/FABS(x); break;
222  case 8 : value=K(40320.0)/POW(FABS(x),K(9.0)); break;
223  case 9 : value=-K(362880.0)/(x*x*x*x*x*x*x*x*x)/FABS(x); break;
224  case 10 : value=K(3628800.0)/POW(FABS(x),K(11.0)); break;
225  case 11 : value=-K(39916800.0)/POW(x,K(11.0))/FABS(x); break;
226  case 12 : value=K(479001600.0)/POW(FABS(x),K(13.0)); break;
227  default : value=K(0.0);
228  }
229 
230  return value;
231 }
232 
233 C one_over_x(R x, int der, const R *param) /* K(x) = 1/x */
234 {
235  R value=K(0.0);
236 
237  (void)param;
238 
239  if (FABS(x)<DBL_EPSILON) value=K(0.0);
240  else switch (der)
241  {
242  case 0 : value=K(1.0)/x; break;
243  case 1 : value=-K(1.0)/(x*x); break;
244  case 2 : value=K(2.0)/(x*x*x); break;
245  case 3 : value=-K(6.0)/(x*x*x*x); break;
246  case 4 : value=K(24.0)/(x*x*x*x*x); break;
247  case 5 : value=-K(120.0)/(x*x*x*x*x*x); break;
248  case 6 : value=K(720.0)/(x*x*x*x*x*x*x); break;
249  case 7 : value=-K(5040.0)/(x*x*x*x*x*x*x*x); break;
250  case 8 : value=K(40320.0)/(x*x*x*x*x*x*x*x*x); break;
251  case 9 : value=-K(362880.0)/POW(x,K(10.0)); break;
252  case 10 : value=K(3628800.0)/POW(x,K(11.0)); break;
253  case 11 : value=-K(39916800.0)/POW(x,K(12.0)); break;
254  case 12 : value=K(479001600.0)/POW(x,K(13.0)); break;
255  default : value=K(0.0);
256  }
257 
258  return value;
259 }
260 
261 C inverse_multiquadric3(R x, int der, const R *param) /* K(x) = 1/SQRT(x^2+c^2)^3 */
262 {
263  R c=param[0];
264  R value=K(0.0);
265 
266  switch (der)
267  {
268  case 0 : value=K(1.0)/(SQRT(POW(x*x+c*c,K(3.0)))); break;
269  case 1 : value=-K(3.0)/SQRT(POW(x*x+c*c,K(5.0)))*x; break;
270  case 2 : value=K(3.0)*(K(4.0)*x*x-c*c)/SQRT(POW(x*x+c*c,K(7.0))); break;
271  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;
272  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;
273  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;
274  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;
275  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;
276  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;
277  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;
278  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;
279  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;
280  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;
281  default : value=K(0.0);
282  }
283 
284  return value;
285 }
286 
287 C sinc_kernel(R x, int der, const R *param) /* K(x) = SIN(cx)/x */
288 {
289  R c=param[0];
290  R value=K(0.0);
291 
292  if (FABS(x)<DBL_EPSILON) value=K(0.0);
293  else switch (der)
294  {
295  case 0 : value=SIN(c*x)/x; break;
296  case 1 : value=(COS(c*x)*c*x-SIN(c*x))/(x*x); break;
297  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;
298  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;
299  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;
300  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;
301  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;
302  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;
303  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;
304  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;
305  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;
306  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;
307  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;
308  default : value=K(0.0);
309  }
310 
311  return value;
312 }
313 
314 C cosc(R x, int der, const R *param) /* K(x) = COS(cx)/x */
315 {
316  R c=param[0];
317  R value=K(0.0);
318  R sign;
319 
320  if (x<0) sign=-K(1.0); else sign=K(1.0);
321  x=FABS(x);
322 
323  if (FABS(x)<DBL_EPSILON) value=K(0.0);
324  else switch (der)
325  {
326  case 0 : value=COS(c*x)/x; break;
327  case 1 : value=-(SIN(c*x)*c*x+COS(c*x))/(x*x); break;
328  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;
329  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;
330  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;
331  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;
332  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;
333  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;
334  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;
335  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;
336  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;
337  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;
338  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;
339  default : value=K(0.0);
340  }
341  value *= POW(sign, (R)(der));
342 
343  return value;
344 }
345 
346 C kcot(R x, int der, const R *param) /* K(x) = cot(cx) */
347 {
348  R c=param[0];
349  R value=K(0.0);
350 
351  if (FABS(x)<DBL_EPSILON) value=K(0.0);
352  else switch (der)
353  {
354  case 0 : value = K(1.0)/TAN(c * x); break;
355  case 1 : value = -(K(1.0) + POW(K(1.0)/TAN(c * x), K(2.0))) * c; break;
356  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;
357  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;
358  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;
359  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;
360  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;
361  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;
362  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;
363  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;
364  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;
365  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;
366  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;
367  default : value=K(0.0);
368  }
369 
370  return value;
371 }
372 
373 
374 C one_over_cube(R x, int der, const R *param)
375 {
376  R value=K(0.0);
377  UNUSED(param);
378 
379  if (FABS(x)<DBL_EPSILON) value=K(0.0);
380  else switch (der)
381  {
382  case 0 : value = K(1.0)/(x*x*x); break;
383  case 1 : value = -K(3.0)/(x*x*x*x); break;
384  case 2 : value = K(12.0)/(x*x*x*x*x); break;
385  case 3 : value = -K(60.0)/(x*x*x*x*x*x); break;
386  case 4 : value = K(360.0)/(x*x*x*x*x*x*x); break;
387  case 5 : value = -K(2520.0)/(x*x*x*x*x*x*x*x); break;
388  case 6 : value = K(20160.0)/POW(x, K(9.0)); break;
389  case 7 : value = -K(181440.0)/POW(x, K(10.0)); break;
390  case 8 : value = K(1814400.0)/POW(x, K(11.0)); break;
391  case 9 : value = -K(19958400.0)/POW(x, K(12.0)); break;
392  case 10 : value = K(239500800.0)/POW(x, K(13.0)); break;
393  case 11 : value = -K(3113510400.0)/POW(x, K(14.0)); break;
394  case 12 : value = K(43589145600.0)/POW(x, K(15.0)); break;
395  default : value=K(0.0);
396  }
397 
398  return value;
399 }
400 
401 
402 C log_sin(R x, int der, const R *param) /* K(x) = log(|sin(cx)|) */
403 {
404  R c=param[0];
405  R value=K(0.0);
406 
407  if (FABS(x)<DBL_EPSILON) value=K(0.0);
408  else
409  {
410  if (der == 0) value = LOG(FABS(SIN(c * x)));
411  else value = c * kcot(x, der-1, param);
412  }
413 
414  return value;
415 }
416 
417 
418 /* \} */
419 
420 /* kernels.c */
inverse_multiquadric3
C inverse_multiquadric3(R x, int der, const R *param)
K(x) = 1/sqrt(x^2+c^2)^3.
Definition: kernels.c:261
one_over_cube
C one_over_cube(R x, int der, const R *param)
K(x) = 1/x^3.
Definition: kernels.c:374
gaussian
C gaussian(R x, int der, const R *param)
K(x)=exp(-x^2/c^2)
Definition: kernels.c:38
kernels.h
Header file with predefined kernels for the fast summation algorithm.
cosc
C cosc(R x, int der, const R *param)
K(x) = cos(cx)/x.
Definition: kernels.c:314
one_over_x
C one_over_x(R x, int der, const R *param)
K(x) = 1/x.
Definition: kernels.c:233
thinplate_spline
C thinplate_spline(R x, int der, const R *param)
K(x) = x^2 log |x|.
Definition: kernels.c:149
sinc_kernel
C sinc_kernel(R x, int der, const R *param)
K(x) = sin(cx)/x.
Definition: kernels.c:287
one_over_square
C one_over_square(R x, int der, const R *param)
K(x) = 1/x^2.
Definition: kernels.c:177
UNUSED
#define UNUSED(x)
Dummy use of unused parameters to silence compiler warnings.
Definition: infft.h:1365
one_over_modulus
C one_over_modulus(R x, int der, const R *param)
K(x) = 1/|x|.
Definition: kernels.c:205
logarithm
C logarithm(R x, int der, const R *param)
K(x)=log |x|.
Definition: kernels.c:116
inverse_multiquadric
C inverse_multiquadric(R x, int der, const R *param)
K(x)=1/sqrt(x^2+c^2)
Definition: kernels.c:90
kcot
C kcot(R x, int der, const R *param)
K(x) = cot(cx)
Definition: kernels.c:346
multiquadric
C multiquadric(R x, int der, const R *param)
K(x)=sqrt(x^2+c^2)
Definition: kernels.c:64
log_sin
C log_sin(R x, int der, const R *param)
K(x) = log(|sin(cx)|)
Definition: kernels.c:402
Generated on Thu Apr 12 2018 21:00:49 for NFFT by   doxygen 1.8.6