NFFT Logo 3.2.3
kernels.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2012 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: kernels.c 3775 2012-06-02 16:39:48Z keiner $ */
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 double _Complex gaussian(double x, int der, const double *param) /* K(x)=exp(-x^2/c^2) */
41 {
42  double c=param[0];
43  double value=0.0;
44 
45  switch (der)
46  {
47  case 0 : value=exp(-x*x/(c*c)); break;
48  case 1 : value=-2.0*x/(c*c)*exp(-x*x/(c*c)); break;
49  case 2 : value=2.0*exp(-x*x/(c*c))*(-c*c+2.0*x*x)/(c*c*c*c); break;
50  case 3 : value=-4.0*x*exp(-x*x/(c*c))*(-3.0*c*c+2.0*x*x)/(c*c*c*c*c*c); break;
51  case 4 : value=4.0*exp(-x*x/(c*c))*(3.0*c*c*c*c-12.0*c*c*x*x+4.0*x*x*x*x)/(c*c*c*c*c*c*c*c); break;
52  case 5 : value=-8.0*x*exp(-x*x/(c*c))*(15.0*c*c*c*c-20.0*c*c*x*x+4.0*x*x*x*x)/pow(c,10.0); break;
53  case 6 : value=8.0*exp(-x*x/(c*c))*(-15.0*c*c*c*c*c*c+90.0*x*x*c*c*c*c-60.0*x*x*x*x*c*c+8.0*x*x*x*x*x*x)/pow(c,12.0); break;
54  case 7 : value=-16.0*x*exp(-x*x/(c*c))*(-105.0*c*c*c*c*c*c+210.0*x*x*c*c*c*c-84.0*x*x*x*x*c*c+8.0*x*x*x*x*x*x)/pow(c,14.0); break;
55  case 8 : value=16.0*exp(-x*x/(c*c))*(105.0*c*c*c*c*c*c*c*c-840.0*x*x*c*c*c*c*c*c+840.0*x*x*x*x*c*c*c*c-224.0*x*x*x*x*x*x*c*c+16.0*x*x*x*x*x*x*x*x)/pow(c,16.0); break;
56  case 9 : value=-32.0*x*exp(-x*x/(c*c))*(945.0*c*c*c*c*c*c*c*c-2520.0*x*x*c*c*c*c*c*c+1512.0*x*x*x*x*c*c*c*c-288.0*x*x*x*x*x*x*c*c+16.0*x*x*x*x*x*x*x*x)/pow(c,18.0); break;
57  case 10 : value=32.0*exp(-x*x/(c*c))*(-945.0*pow(c,10.0)+9450.0*x*x*c*c*c*c*c*c*c*c-12600.0*x*x*x*x*c*c*c*c*c*c+5040.0*x*x*x*x*x*x*c*c*c*c-720.0*x*x*x*x*x*x*x*x*c*c+32.0*pow(x,10.0))/pow(c,20.0); break;
58  case 11 : value=-64.0*x*exp(-x*x/(c*c))*(-10395.0*pow(c,10.0)+34650.0*x*x*c*c*c*c*c*c*c*c-27720.0*x*x*x*x*c*c*c*c*c*c+7920.0*x*x*x*x*x*x*c*c*c*c-880.0*x*x*x*x*x*x*x*x*c*c+32.0*pow(x,10.0))/pow(c,22.0); break;
59  case 12 : value=64.0*exp(-x*x/(c*c))*(10395.0*pow(c,12.0)-124740.0*x*x*pow(c,10.0)+207900.0*x*x*x*x*c*c*c*c*c*c*c*c-110880.0*x*x*x*x*x*x*c*c*c*c*c*c+23760.0*x*x*x*x*x*x*x*x*c*c*c*c-2112.0*pow(x,10.0)*c*c+64.0*pow(x,12.0))/pow(c,24.0); break;
60  default : value=0.0;
61  }
62 
63  return value;
64 }
65 
66 double _Complex multiquadric(double x, int der, const double *param) /* K(x)=sqrt(x^2+c^2) */
67 {
68  double c=param[0];
69  double value=0.0;
70 
71  switch (der)
72  {
73  case 0 : value=sqrt(x*x+c*c); break;
74  case 1 : value=1.0/(sqrt(x*x+c*c))*x; break;
75  case 2 : value=c*c/sqrt(pow(x*x+c*c,3.0)); break;
76  case 3 : value=-3.0*x*c*c/sqrt(pow(x*x+c*c,5.0)); break;
77  case 4 : value=3.0*c*c*(4.0*x*x-c*c)/sqrt(pow(x*x+c*c,7.0)); break;
78  case 5 : value=-15.0*x*c*c*(4.0*x*x-3.0*c*c)/sqrt(pow(x*x+c*c,9.0)); break;
79  case 6 : value=45.0*c*c*(8.0*x*x*x*x-12.0*x*x*c*c+c*c*c*c)/sqrt(pow(x*x+c*c,11.0)); break;
80  case 7 : value=-315.0*x*c*c*(8.0*x*x*x*x-20.0*x*x*c*c+5.0*c*c*c*c)/sqrt(pow(x*x+c*c,13.0)); break;
81  case 8 : value=315.0*c*c*(64.0*x*x*x*x*x*x-240.0*x*x*x*x*c*c+120.0*x*x*c*c*c*c-5.0*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,15.0)); break;
82  case 9 : value=-2835.0*x*c*c*(64.0*x*x*x*x*x*x-336.0*x*x*x*x*c*c+280.0*x*x*c*c*c*c-35.0*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,17.0)); break;
83  case 10 : value=14175.0*c*c*(128.0*x*x*x*x*x*x*x*x-896.0*x*x*x*x*x*x*c*c+1120.0*x*x*x*x*c*c*c*c-280.0*x*x*c*c*c*c*c*c+7.0*c*c*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,19.0)); break;
84  case 11 : value=-155925.0*x*c*c*(128.0*x*x*x*x*x*x*x*x-1152.0*x*x*x*x*x*x*c*c+2016.0*x*x*x*x*c*c*c*c-840.0*x*x*c*c*c*c*c*c+63.0*c*c*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,21.0)); break;
85  case 12 : value=467775.0*c*c*(1260.0*x*x*c*c*c*c*c*c*c*c-21.0*pow(c,10.0)+512.0*pow(x,10.0)-5760.0*x*x*x*x*x*x*x*x*c*c+13440.0*x*x*x*x*x*x*c*c*c*c-8400.0*x*x*x*x*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,23.0)); break;
86  default : value=0.0;
87  }
88 
89  return value;
90 }
91 
92 double _Complex inverse_multiquadric(double x, int der, const double *param) /* K(x)=1/sqrt(x^2+c^2) */
93 {
94  double c=param[0];
95  double value=0.0;
96 
97  switch (der)
98  {
99  case 0 : value=1.0/sqrt(x*x+c*c); break;
100  case 1 : value=-1.0/(sqrt(pow(x*x+c*c,3.0)))*x; break;
101  case 2 : value=(2.0*x*x-c*c)/sqrt(pow(x*x+c*c,5.0)); break;
102  case 3 : value=-3.0*x*(2.0*x*x-3.0*c*c)/sqrt(pow(x*x+c*c,7.0)); break;
103  case 4 : value=3.0*(8.0*x*x*x*x-24.0*x*x*c*c+3.0*c*c*c*c)/sqrt(pow(x*x+c*c,9.0)); break;
104  case 5 : value=-15.0*x*(8.0*x*x*x*x-40.0*x*x*c*c+15.0*c*c*c*c)/sqrt(pow(x*x+c*c,11.0)); break;
105  case 6 : value=45.0*(16.0*x*x*x*x*x*x-120.0*x*x*x*x*c*c+90.0*x*x*c*c*c*c-5.0*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,13.0)); break;
106  case 7 : value=-315.0*x*(16.0*x*x*x*x*x*x-168.0*x*x*x*x*c*c+210.0*x*x*c*c*c*c-35.0*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,15.0)); break;
107  case 8 : value=315.0*(128.0*x*x*x*x*x*x*x*x-1792.0*x*x*x*x*x*x*c*c+3360.0*x*x*x*x*c*c*c*c-1120.0*x*x*c*c*c*c*c*c+35.0*c*c*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,17.0)); break;
108  case 9 : value=-2835.0*x*(128.0*x*x*x*x*x*x*x*x-2304.0*x*x*x*x*x*x*c*c+6048.0*x*x*x*x*c*c*c*c-3360.0*x*x*c*c*c*c*c*c+315.0*c*c*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,19.0)); break;
109  case 10 : value=14175.0*(256.0*pow(x,10.0)-5760.0*x*x*x*x*x*x*x*x*c*c+20160.0*x*x*x*x*x*x*c*c*c*c-16800.0*x*x*x*x*c*c*c*c*c*c+3150.0*x*x*c*c*c*c*c*c*c*c-63.0*pow(c,10.0))/sqrt(pow(x*x+c*c,21.0)); break;
110  case 11 : value=-155925.0*x*(256.0*pow(x,10.0)-7040.0*x*x*x*x*x*x*x*x*c*c+31680.0*x*x*x*x*x*x*c*c*c*c-36960.0*x*x*x*x*c*c*c*c*c*c+11550.0*x*x*c*c*c*c*c*c*c*c-693.0*pow(c,10.0))/sqrt(pow(x*x+c*c,23.0)); break;
111  case 12 : value=467775.0*(231.0*pow(c,12.0)+190080.0*x*x*x*x*x*x*x*x*c*c*c*c-16632.0*x*x*pow(c,10.0)-295680.0*x*x*x*x*x*x*c*c*c*c*c*c+138600.0*x*x*x*x*c*c*c*c*c*c*c*c+1024.0*pow(x,12.0)-33792.0*pow(x,10.0)*c*c)/sqrt(pow(x*x+c*c,25.0)); break;
112  default : value=0.0;
113  }
114 
115  return value;
116 }
117 
118 double _Complex logarithm(double x, int der, const double *param) /* K(x)=log |x| */
119 {
120  double value=0.0;
121 
122  (void)param;
123 
124  if (fabs(x)<DBL_EPSILON) value=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=2.0*(x<0 ? -1 : 1)/pow(fabs(x),3.0); break;
131  case 4 : value=-6.0/(x*x*x*x); break;
132  case 5 : value=24.0*(x<0 ? -1 : 1)/pow(fabs(x),5.0); break;
133  case 6 : value=-120.0/(x*x*x*x*x*x); break;
134  case 7 : value=720.0*(x<0 ? -1 : 1)/pow(fabs(x),7.0); break;
135  case 8 : value=-5040.0/(x*x*x*x*x*x*x*x); break;
136  case 9 : value=40320.0*(x<0 ? -1 : 1)/pow(fabs(x),9.0); break;
137  case 10 : value=-362880.0/pow(x,10.0); break;
138  case 11 : value=3628800.0*(x<0 ? -1 : 1)/pow(fabs(x),11.0); break;
139  case 12 : value=-39916800.0/pow(x,12.0); break;
140  case 13 : value=479001600.0/pow(x,13.0); break;
141  case 14 : value=-6227020800.0/pow(x,14.0); break;
142  case 15 : value=87178291200.0/pow(x,15.0); break;
143  case 16 : value=-1307674368000.0/pow(x,16.0); break;
144  case 17 : value=20922789888000.0/pow(x,17.0); break;
145  default : value=0.0;
146  }
147 
148  return value;
149 }
150 
151 double _Complex thinplate_spline(double x, int der, const double *param) /* K(x) = x^2 log |x| */
152 {
153  double value=0.0;
154 
155  (void)param;
156 
157  if (fabs(x)<DBL_EPSILON) value=0.0;
158  else switch (der)
159  {
160  case 0 : value=x*x*log(fabs(x)); break;
161  case 1 : value=2.0*x*log(fabs(x))+x; break;
162  case 2 : value=2.0*log(fabs(x))+3.0; break;
163  case 3 : value=2.0/x; break;
164  case 4 : value=-2.0/(x*x); break;
165  case 5 : value=4.0/(x*x*x); break;
166  case 6 : value=-12.0/(x*x*x*x); break;
167  case 7 : value=48.0/(x*x*x*x*x); break;
168  case 8 : value=-240.0/(x*x*x*x*x*x); break;
169  case 9 : value=1440.0/(x*x*x*x*x*x*x); break;
170  case 10 : value=-10080.0/(x*x*x*x*x*x*x*x); break;
171  case 11 : value=80640.0/(x*x*x*x*x*x*x*x*x); break;
172  case 12 : value=-725760.0/pow(x,10.0); break;
173  default : value=0.0;
174  }
175 
176  return value;
177 }
178 
179 double _Complex one_over_square(double x, int der, const double *param) /* K(x) = 1/x^2 */
180 {
181  double value=0.0;
182 
183  (void)param;
184 
185  if (fabs(x)<DBL_EPSILON) value=0.0;
186  else switch (der)
187  {
188  case 0 : value=1.0/(x*x); break;
189  case 1 : value=-2.0/(x*x*x); break;
190  case 2 : value=6.0/(x*x*x*x); break;
191  case 3 : value=-24.0/(x*x*x*x*x); break;
192  case 4 : value=120.0/(x*x*x*x*x*x); break;
193  case 5 : value=-720.0/(x*x*x*x*x*x*x); break;
194  case 6 : value=5040.0/(x*x*x*x*x*x*x*x); break;
195  case 7 : value=-40320.0/(x*x*x*x*x*x*x*x*x); break;
196  case 8 : value=362880.0/pow(x,10.0); break;
197  case 9 : value=-3628800.0/pow(x,11.0); break;
198  case 10 : value=39916800.0/pow(x,12.0); break;
199  case 11 : value=-479001600.0/pow(x,13.0); break;
200  case 12 : value=6227020800.0/pow(x,14.0); break;
201  default : value=0.0;
202  }
203 
204  return value;
205 }
206 
207 double _Complex one_over_modulus(double x, int der, const double *param) /* K(x) = 1/|x| */
208 {
209  double value=0.0;
210 
211  (void)param;
212 
213  if (fabs(x)<DBL_EPSILON) value=0.0;
214  else switch (der)
215  {
216  case 0 : value=1.0/fabs(x); break;
217  case 1 : value=-1/x/fabs(x); break;
218  case 2 : value=2.0/pow(fabs(x),3.0); break;
219  case 3 : value=-6.0/(x*x*x)/fabs(x); break;
220  case 4 : value=24.0/pow(fabs(x),5.0); break;
221  case 5 : value=-120.0/(x*x*x*x*x)/fabs(x); break;
222  case 6 : value=720.0/pow(fabs(x),7.0); break;
223  case 7 : value=-5040.0/(x*x*x*x*x*x*x)/fabs(x); break;
224  case 8 : value=40320.0/pow(fabs(x),9.0); break;
225  case 9 : value=-362880.0/(x*x*x*x*x*x*x*x*x)/fabs(x); break;
226  case 10 : value=3628800.0/pow(fabs(x),11.0); break;
227  case 11 : value=-39916800.0/pow(x,11.0)/fabs(x); break;
228  case 12 : value=479001600.0/pow(fabs(x),13.0); break;
229  default : value=0.0;
230  }
231 
232  return value;
233 }
234 
235 double _Complex one_over_x(double x, int der, const double *param) /* K(x) = 1/x */
236 {
237  double value=0.0;
238 
239  (void)param;
240 
241  if (fabs(x)<DBL_EPSILON) value=0.0;
242  else switch (der)
243  {
244  case 0 : value=1.0/x; break;
245  case 1 : value=-1.0/(x*x); break;
246  case 2 : value=2.0/(x*x*x); break;
247  case 3 : value=-6.0/(x*x*x*x); break;
248  case 4 : value=24.0/(x*x*x*x*x); break;
249  case 5 : value=-120.0/(x*x*x*x*x*x); break;
250  case 6 : value=720.0/(x*x*x*x*x*x*x); break;
251  case 7 : value=-5040.0/(x*x*x*x*x*x*x*x); break;
252  case 8 : value=40320.0/(x*x*x*x*x*x*x*x*x); break;
253  case 9 : value=-362880.0/pow(x,10.0); break;
254  case 10 : value=3628800.0/pow(x,11.0); break;
255  case 11 : value=-39916800.0/pow(x,12.0); break;
256  case 12 : value=479001600.0/pow(x,13.0); break;
257  default : value=0.0;
258  }
259 
260  return value;
261 }
262 
263 double _Complex inverse_multiquadric3(double x, int der, const double *param) /* K(x) = 1/sqrt(x^2+c^2)^3 */
264 {
265  double c=param[0];
266  double value=0.0;
267 
268  switch (der)
269  {
270  case 0 : value=1.0/(sqrt(pow(x*x+c*c,3.0))); break;
271  case 1 : value=-3.0/sqrt(pow(x*x+c*c,5.0))*x; break;
272  case 2 : value=3.0*(4.0*x*x-c*c)/sqrt(pow(x*x+c*c,7.0)); break;
273  case 3 : value=-15.0*x*(4.0*x*x-3.0*c*c)/sqrt(pow(x*x+c*c,9.0)); break;
274  case 4 : value=45.0*(8.0*x*x*x*x-12.0*x*x*c*c+c*c*c*c)/sqrt(pow(x*x+c*c,11.0)); break;
275  case 5 : value=-315.0*x*(8.0*x*x*x*x-20.0*x*x*c*c+5.0*c*c*c*c)/sqrt(pow(x*x+c*c,13.0)); break;
276  case 6 : value=315.0*(64.0*x*x*x*x*x*x-240.0*x*x*x*x*c*c+120.0*x*x*c*c*c*c-5.0*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,15.0)); break;
277  case 7 : value=-2835.0*x*(64.0*x*x*x*x*x*x-336.0*x*x*x*x*c*c+280.0*x*x*c*c*c*c-35.0*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,17.0)); break;
278  case 8 : value=14175.0*(128.0*x*x*x*x*x*x*x*x-896.0*x*x*x*x*x*x*c*c+1120.0*x*x*x*x*c*c*c*c-280.0*x*x*c*c*c*c*c*c+7.0*c*c*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,19.0)); break;
279  case 9 : value=-155925.0*x*(128.0*x*x*x*x*x*x*x*x-1152.0*x*x*x*x*x*x*c*c+2016.0*x*x*x*x*c*c*c*c-840.0*x*x*c*c*c*c*c*c+63.0*c*c*c*c*c*c*c*c)/sqrt(pow(x*x+c*c,21.0)); break;
280  case 10 : value=467775.0*(512.0*pow(x,10.0)-5760.0*x*x*x*x*x*x*x*x*c*c+13440.0*x*x*x*x*x*x*c*c*c*c-8400.0*x*x*x*x*c*c*c*c*c*c+1260.0*x*x*c*c*c*c*c*c*c*c-21.0*pow(c,10.0))/sqrt(pow(x*x+c*c,23.0)); break;
281  case 11 : value=-6081075.0*x*(512.0*pow(x,10.0)-7040.0*x*x*x*x*x*x*x*x*c*c+21120.0*x*x*x*x*x*x*c*c*c*c-18480.0*x*x*x*x*c*c*c*c*c*c+4620.0*x*x*c*c*c*c*c*c*c*c-231.0*pow(c,10.0))/sqrt(pow(x*x+c*c,25.0)); break;
282  case 12 : value=42567525.0*(1024.0*pow(x,12.0)+27720.0*x*x*x*x*c*c*c*c*c*c*c*c+33.0*pow(c,12.0)-2772.0*x*x*pow(c,10.0)-73920.0*x*x*x*x*x*x*c*c*c*c*c*c+63360.0*x*x*x*x*x*x*x*x*c*c*c*c-16896.0*pow(x,10.0)*c*c)/sqrt(pow(x*x+c*c,27.0)); break;
283  default : value=0.0;
284  }
285 
286  return value;
287 }
288 
289 double _Complex sinc_kernel(double x, int der, const double *param) /* K(x) = sin(cx)/x */
290 {
291  double c=param[0];
292  double value=0.0;
293 
294  if (fabs(x)<DBL_EPSILON) value=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+2.0*cos(c*x)*c*x-2.0*sin(c*x))/(x*x*x); break;
300  case 3 : value=-(cos(c*x)*c*c*c*x*x*x-3.0*sin(c*x)*c*c*x*x-6.0*cos(c*x)*c*x+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+4.0*cos(c*x)*c*c*c*x*x*x-12.0*sin(c*x)*c*c*x*x-24.0*cos(c*x)*c*x+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-5.0*sin(c*x)*c*c*c*c*x*x*x*x-20.0*cos(c*x)*c*c*c*x*x*x+60.0*sin(c*x)*c*c*x*x+120.0*cos(c*x)*c*x-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+6.0*cos(c*x)*c*c*c*c*c*x*x*x*x*x-30.0*sin(c*x)*c*c*c*c*x*x*x*x-120.0*cos(c*x)*c*c*c*x*x*x+360.0*sin(c*x)*c*c*x*x+720.0*cos(c*x)*c*x-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-7.0*sin(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-42.0*cos(c*x)*c*c*c*c*c*x*x*x*x*x+210.0*sin(c*x)*c*c*c*c*x*x*x*x+840.0*cos(c*x)*c*c*c*x*x*x-2520.0*sin(c*x)*c*c*x*x-5040.0*cos(c*x)*c*x+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+8.0*cos(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-56.0*sin(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-336.0*cos(c*x)*c*c*c*c*c*x*x*x*x*x+1680.0*sin(c*x)*c*c*c*c*x*x*x*x+6720.0*cos(c*x)*c*c*c*x*x*x-20160.0*sin(c*x)*c*c*x*x-40320.0*cos(c*x)*c*x+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-9.0*sin(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-72.0*cos(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+504.0*sin(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+3024.0*cos(c*x)*c*c*c*c*c*x*x*x*x*x-15120.0*sin(c*x)*c*c*c*c*x*x*x*x-60480.0*cos(c*x)*c*c*c*x*x*x+181440.0*sin(c*x)*c*c*x*x+362880.0*cos(c*x)*c*x-362880.0*sin(c*x))/pow(x,10.0); break;
307  case 10 : value=-(sin(c*x)*pow(c,10.0)*pow(x,10.0)+10.0*cos(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-90.0*sin(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-720.0*cos(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+5040.0*sin(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+30240.0*cos(c*x)*c*c*c*c*c*x*x*x*x*x-151200.0*sin(c*x)*c*c*c*c*x*x*x*x-604800.0*cos(c*x)*c*c*c*x*x*x+1814400.0*sin(c*x)*c*c*x*x+3628800.0*cos(c*x)*c*x-3628800.0*sin(c*x))/pow(x,11.0); break;
308  case 11 : value=-(cos(c*x)*pow(c,11.0)*pow(x,11.0)-11.0*sin(c*x)*pow(c,10.0)*pow(x,10.0)-110.0*cos(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+990.0*sin(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+7920.0*cos(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-55440.0*sin(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-332640.0*cos(c*x)*c*c*c*c*c*x*x*x*x*x+1663200.0*sin(c*x)*c*c*c*c*x*x*x*x+6652800.0*cos(c*x)*c*c*c*x*x*x-19958400.0*sin(c*x)*c*c*x*x-39916800.0*cos(c*x)*c*x+39916800.0*sin(c*x))/pow(x,12.0); break;
309  case 12 : value=(sin(c*x)*pow(c,12.0)*pow(x,12.0)+12.0*cos(c*x)*pow(c,11.0)*pow(x,11.0)-132.0*sin(c*x)*pow(c,10.0)*pow(x,10.0)-1320.0*cos(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+11880.0*sin(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+95040.0*cos(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-665280.0*sin(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-3991680.0*cos(c*x)*c*c*c*c*c*x*x*x*x*x+19958400.0*sin(c*x)*c*c*c*c*x*x*x*x+79833600.0*cos(c*x)*c*c*c*x*x*x-239500800.0*sin(c*x)*c*c*x*x-479001600.0*cos(c*x)*c*x+479001600.0*sin(c*x))/pow(x,13.0); break;
310  default : value=0.0;
311  }
312 
313  return value;
314 }
315 
316 double _Complex cosc(double x, int der, const double *param) /* K(x) = cos(cx)/x */
317 {
318  double c=param[0];
319  double value=0.0;
320  double sign;
321 
322  if (x<0) sign=-1.0; else sign=1.0;
323  x=fabs(x);
324 
325  if (fabs(x)<DBL_EPSILON) value=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+2.0*sin(c*x)*c*x+2.0*cos(c*x))/(x*x*x); break;
331  case 3 : value=(sin(c*x)*c*c*c*x*x*x+3.0*cos(c*x)*c*c*x*x-6.0*sin(c*x)*c*x-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-4.0*sin(c*x)*c*c*c*x*x*x-12.0*cos(c*x)*c*c*x*x+24.0*sin(c*x)*c*x+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+5.0*cos(c*x)*c*c*c*c*x*x*x*x-20.0*sin(c*x)*c*c*c*x*x*x-60.0*cos(c*x)*c*c*x*x+120.0*sin(c*x)*c*x+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-6.0*sin(c*x)*c*c*c*c*c*x*x*x*x*x-30.0*cos(c*x)*c*c*c*c*x*x*x*x+120.0*sin(c*x)*c*c*c*x*x*x+360.0*cos(c*x)*c*c*x*x-720.0*sin(c*x)*c*x-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+7.0*cos(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-42.0*sin(c*x)*c*c*c*c*c*x*x*x*x*x-210.0*cos(c*x)*c*c*c*c*x*x*x*x+840.0*sin(c*x)*c*c*c*x*x*x+2520.0*cos(c*x)*c*c*x*x-5040.0*sin(c*x)*c*x-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-8.0*sin(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-56.0*cos(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+336.0*sin(c*x)*c*c*c*c*c*x*x*x*x*x+1680.0*cos(c*x)*c*c*c*c*x*x*x*x-6720.0*sin(c*x)*c*c*c*x*x*x-20160.0*cos(c*x)*c*c*x*x+40320.0*sin(c*x)*c*x+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+9.0*cos(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-72.0*sin(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-504.0*cos(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+3024.0*sin(c*x)*c*c*c*c*c*x*x*x*x*x+15120.0*cos(c*x)*c*c*c*c*x*x*x*x-60480.0*sin(c*x)*c*c*c*x*x*x-181440.0*cos(c*x)*c*c*x*x+362880.0*sin(c*x)*c*x+362880.0*cos(c*x))/pow(x,10.0); break;
338  case 10 : value=-(cos(c*x)*pow(c,10.0)*pow(x,10.0)-10.0*sin(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-90.0*cos(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+720.0*sin(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+5040.0*cos(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-30240.0*sin(c*x)*c*c*c*c*c*x*x*x*x*x-151200.0*cos(c*x)*c*c*c*c*x*x*x*x+604800.0*sin(c*x)*c*c*c*x*x*x+1814400.0*cos(c*x)*c*c*x*x-3628800.0*sin(c*x)*c*x-3628800.0*cos(c*x))/pow(x,11.0); break;
339  case 11 : value=(sin(c*x)*pow(c,11.0)*pow(x,11.0)+11.0*cos(c*x)*pow(c,10.0)*pow(x,10.0)-110.0*sin(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-990.0*cos(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+7920.0*sin(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+55440.0*cos(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-332640.0*sin(c*x)*c*c*c*c*c*x*x*x*x*x-1663200.0*cos(c*x)*c*c*c*c*x*x*x*x+6652800.0*sin(c*x)*c*c*c*x*x*x+19958400.0*cos(c*x)*c*c*x*x-39916800.0*sin(c*x)*c*x-39916800.0*cos(c*x))/pow(x,12.0); break;
340  case 12 : value=(cos(c*x)*pow(c,12.0)*pow(x,12.0)-12.0*sin(c*x)*pow(c,11.0)*pow(x,11.0)-132.0*cos(c*x)*pow(c,10.0)*pow(x,10.0)+1320.0*sin(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+11880.0*cos(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-95040.0*sin(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-665280.0*cos(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+3991680.0*sin(c*x)*c*c*c*c*c*x*x*x*x*x+19958400.0*cos(c*x)*c*c*c*c*x*x*x*x-79833600.0*sin(c*x)*c*c*c*x*x*x-239500800.0*cos(c*x)*c*c*x*x+479001600.0*sin(c*x)*c*x+479001600.0*cos(c*x))/pow(x,13.0); break;
341  default : value=0.0;
342  }
343  value*=pow(sign,der);
344 
345  return value;
346 }
347 
348 double _Complex kcot(double x, int der, const double *param) /* K(x) = cot(cx) */
349 {
350  double c=param[0];
351  double value=0.0;
352 
353  if (fabs(x)<DBL_EPSILON) value=0.0;
354  else switch (der)
355  {
356  case 0 : value = 1.0/tan(c * x); break;
357  case 1 : value = -(1.0 + pow(1.0/tan(c * x), 2.0)) * c; break;
358  case 2 : value = 2.0 * 1.0/tan(c * x) * (1.0 + pow(1.0/tan(c * x), 2.0)) * c * c; break;
359  case 3 : value = -2.0 * (1.0 + pow(1.0/tan(c * x), 2.0)) * pow(c, 3.0) * (1.0 + 3.0 * pow(1.0/tan(c * x), 2.0)); break;
360  case 4 : value = 8.0 * (1.0 + pow(1.0/tan(c * x), 2.0)) * pow(c, 4.0) * 1.0/tan(c * x) * (2.0 + 3.0 * pow(1.0/tan(c * x), 2.0)); break;
361  case 5 : value = -0.8e1 * (0.1e1 + pow(1.0/tan(c * x), 0.2e1)) * pow(c, 0.5e1) * (0.15e2 * pow(1.0/tan(c * x), 0.2e1) + 0.15e2 * pow(1.0/tan(c * x), 0.4e1) + 0.2e1); break;
362  case 6 : value = 0.16e2 * (0.1e1 + pow(1.0/tan(c * x), 0.2e1)) * pow(c, 0.6e1) * 1.0/tan(c * x) * (0.60e2 * pow(1.0/tan(c * x), 0.2e1) + 0.45e2 * pow(1.0/tan(c * x), 0.4e1) + 0.17e2); break;
363  case 7 : value = -0.16e2 * (0.1e1 + pow(1.0/tan(c * x), 0.2e1)) * pow(c, 0.7e1) * (0.525e3 * pow(1.0/tan(c * x), 0.4e1) + 0.315e3 * pow(1.0/tan(c * x), 0.6e1) + 0.231e3 * pow(1.0/tan(c * x), 0.2e1) + 0.17e2); break;
364  case 8 : value = 0.128e3 * (0.1e1 + pow(1.0/tan(c * x), 0.2e1)) * pow(c, 0.8e1) * 1.0/tan(c * x) * (0.630e3 * pow(1.0/tan(c * x), 0.4e1) + 0.315e3 * pow(1.0/tan(c * x), 0.6e1) + 0.378e3 * pow(1.0/tan(c * x), 0.2e1) + 0.62e2); break;
365  case 9 : value = -0.128e3 * (0.1e1 + pow(1.0/tan(c * x), 0.2e1)) * pow(c, 0.9e1) * (0.6615e4 * pow(1.0/tan(c * x), 0.6e1) + 0.2835e4 * pow(1.0/tan(c * x), 0.8e1) + 0.5040e4 * pow(1.0/tan(c * x), 0.4e1) + 0.1320e4 * pow(1.0/tan(c * x), 0.2e1) + 0.62e2); break;
366  case 10 : value = 0.256e3 * (0.1e1 + pow(1.0/tan(c * x), 0.2e1)) * pow(c, 0.10e2) * 1.0/tan(c * x) * (0.37800e5 * pow(1.0/tan(c * x), 0.6e1) + 0.14175e5 * pow(1.0/tan(c * x), 0.8e1) + 0.34965e5 * pow(1.0/tan(c * x), 0.4e1) + 0.12720e5 * pow(1.0/tan(c * x), 0.2e1) + 0.1382e4); break;
367  case 11 : value = -0.256e3 * (0.1e1 + pow(1.0/tan(c * x), 0.2e1)) * pow(c, 0.11e2) * (0.467775e6 * pow(1.0/tan(c * x), 0.8e1) + 0.155925e6 * pow(1.0/tan(c * x), 0.10e2) + 0.509355e6 * pow(1.0/tan(c * x), 0.6e1) + 0.238425e6 * pow(1.0/tan(c * x), 0.4e1) + 0.42306e5 * pow(1.0/tan(c * x), 0.2e1) + 0.1382e4); break;
368  case 12 : value = 0.1024e4 * (0.1e1 + pow(1.0/tan(c * x), 0.2e1)) * pow(c, 0.12e2) * 1.0/tan(c * x) * (0.1559250e7 * pow(1.0/tan(c * x), 0.8e1) + 0.467775e6 * pow(1.0/tan(c * x), 0.10e2) + 0.1954260e7 * pow(1.0/tan(c * x), 0.6e1) + 0.1121670e7 * pow(1.0/tan(c * x), 0.4e1) + 0.280731e6 * pow(1.0/tan(c * x), 0.2e1) + 0.21844e5); break;
369  default : value=0.0;
370  }
371 
372  return value;
373 }
374 
375 
376 double _Complex one_over_cube(double x, int der, const double *param)
377 {
378  double value=0.0;
379 
380  if (fabs(x)<DBL_EPSILON) value=0.0;
381  else switch (der)
382  {
383  case 0 : value = 1.0/(x*x*x); break;
384  case 1 : value = -3.0/(x*x*x*x); break;
385  case 2 : value = 12.0/(x*x*x*x*x); break;
386  case 3 : value = -60.0/(x*x*x*x*x*x); break;
387  case 4 : value = 360.0/(x*x*x*x*x*x*x); break;
388  case 5 : value = -2520.0/(x*x*x*x*x*x*x*x); break;
389  case 6 : value = 20160.0/pow(x, 9.0); break;
390  case 7 : value = -181440.0/pow(x, 10.0); break;
391  case 8 : value = 1814400.0/pow(x, 11.0); break;
392  case 9 : value = -19958400.0/pow(x, 12.0); break;
393  case 10 : value = 239500800.0/pow(x, 13.0); break;
394  case 11 : value = -3113510400.0/pow(x, 14.0); break;
395  case 12 : value = 43589145600.0/pow(x, 15.0); break;
396  default : value=0.0;
397  }
398 
399  return value;
400 }
401 
402 
403 /* \} */
404 
405 /* kernels.c */
Generated on Tue Apr 30 2013 by Doxygen 1.8.1