NFFT 3.6.0
infft.h
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#ifndef __INFFT_H__
23#define __INFFT_H__
24
25#include "config.h"
26
27#include <math.h>
28#include <float.h>
29#ifdef HAVE_COMPLEX_H
30#include <complex.h>
31#endif
32#include <stdio.h>
33#include <string.h>
34
35#include <stdlib.h> /* size_t */
36#include <stdarg.h> /* va_list */
37#include <stddef.h> /* ptrdiff_t */
38
39#if HAVE_SYS_TYPES_H
40#include <sys/types.h>
41#endif
42
43#if HAVE_STDINT_H
44#include <stdint.h> /* uintptr_t, maybe */
45#endif
46
47#if HAVE_INTTYPES_H
48#include <inttypes.h> /* uintptr_t, maybe */
49#endif
50
51#include <fftw3.h>
52
53#include "ticks.h"
54
66/* Determine precision and name-mangling scheme. */
67#define CONCAT(prefix, name) prefix ## name
68#if defined(NFFT_SINGLE)
69typedef float R;
70typedef float _Complex C;
71#define Y(name) CONCAT(nfftf_,name)
72#define FFTW(name) CONCAT(fftwf_,name)
73#define NFFT(name) CONCAT(nfftf_,name)
74#define NFCT(name) CONCAT(nfctf_,name)
75#define NFST(name) CONCAT(nfstf_,name)
76#define NFSFT(name) CONCAT(nfsftf_,name)
77#define SOLVER(name) CONCAT(solverf_,name)
78#elif defined(NFFT_LDOUBLE)
79typedef long double R;
80typedef long double _Complex C;
81#define Y(name) CONCAT(nfftl_,name)
82#define FFTW(name) CONCAT(fftwl_,name)
83#define NFFT(name) CONCAT(nfftl_,name)
84#define NFCT(name) CONCAT(nfctl_,name)
85#define NFST(name) CONCAT(nfstl_,name)
86#define NFSFT(name) CONCAT(nfsftl_,name)
87#define SOLVER(name) CONCAT(solverl_,name)
88#else
89typedef double R;
90typedef double _Complex C;
91#define Y(name) CONCAT(nfft_,name)
92#define FFTW(name) CONCAT(fftw_,name)
93#define NFFT(name) CONCAT(nfft_,name)
94#define NFCT(name) CONCAT(nfct_,name)
95#define NFST(name) CONCAT(nfst_,name)
96#define NFSFT(name) CONCAT(nfsft_,name)
97#define SOLVER(name) CONCAT(solver_,name)
98#endif
99#define X(name) Y(name)
100
101#define STRINGIZEx(x) #x
102#define STRINGIZE(x) STRINGIZEx(x)
103
104#ifdef NFFT_LDOUBLE
105# define K(x) ((R) x##L)
106#else
107# define K(x) ((R) x)
108#endif
109#define DK(name, value) const R name = K(value)
110
111#if defined __CYGWIN32__ && !defined __CYGWIN__
112 /* For backwards compatibility with Cygwin b19 and
113 earlier, we define __CYGWIN__ here, so that
114 we can rely on checking just for that macro. */
115# define __CYGWIN__ __CYGWIN32__
116#endif
117
118/* Integral type large enough to contain a stride (what ``int'' should have been
119 * in the first place) */
120typedef ptrdiff_t INT;
121
122#if defined(NFFT_LDOUBLE)
123 #define MANT_DIG LDBL_MANT_DIG
124 #define MIN_EXP LDBL_MIN_EXP
125 #define MAX_EXP LDBL_MAX_EXP
126 #define EPSILON LDBL_EPSILON
127#elif defined(NFFT_SINGLE)
128 #define MANT_DIG FLT_MANT_DIG
129 #define MIN_EXP FLT_MIN_EXP
130 #define MAX_EXP FLT_MAX_EXP
131 #define EPSILON FLT_EPSILON
132#else
133 #define MANT_DIG DBL_MANT_DIG
134 #define MIN_EXP DBL_MIN_EXP
135 #define MAX_EXP DBL_MAX_EXP
136 #define EPSILON DBL_EPSILON
137#endif
138
139#define KPI K(3.1415926535897932384626433832795028841971693993751)
140#define K2PI K(6.2831853071795864769252867665590057683943387987502)
141#define K4PI K(12.5663706143591729538505735331180115367886775975004)
142#define KE K(2.7182818284590452353602874713526624977572470937000)
143
144#define IF(x,a,b) ((x)?(a):(b))
145#define MIN(a,b) (((a)<(b))?(a):(b))
146#define MAX(a,b) (((a)>(b))?(a):(b))
147#define ABS(x) (((x)>K(0.0))?(x):(-(x)))
148#define SIGN(a) (((a)>=0)?1:-1)
149#define SIGN(a) (((a)>=0)?1:-1)
150#define SIGNF(a) IF((a)<K(0.0),K(-1.0),K(1.0))
151
152/* Size of array. */
153#define SIZE(x) sizeof(x)/sizeof(x[0])
154
156#define CSWAP(x,y) {C* NFFT_SWAP_temp__; \
157 NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;}
158
160#define RSWAP(x,y) {R* NFFT_SWAP_temp__; NFFT_SWAP_temp__=(x); \
161 (x)=(y); (y)=NFFT_SWAP_temp__;}
162
163/* macros for window functions */
164
165#if defined(DIRAC_DELTA)
166 #define PHI_HUT(n,k,d) K(1.0)
167 #define PHI(n,x,d) IF(FABS((x)) < K(10E-8),K(1.0),K(0.0))
168 #define WINDOW_HELP_INIT(d)
169 #define WINDOW_HELP_FINALIZE
170 #define WINDOW_HELP_ESTIMATE_m 0
171#elif defined(GAUSSIAN)
172 #define PHI_HUT(n,k,d) ((R)EXP(-(POW(KPI*(k)/n,K(2.0))*ths->b[d])))
173 #define PHI(n,x,d) ((R)EXP(-POW((x)*((R)n),K(2.0)) / \
174 ths->b[d])/SQRT(KPI*ths->b[d]))
175 #define WINDOW_HELP_INIT \
176 { \
177 int WINDOW_idx; \
178 ths->b = (R*) Y(malloc)(ths->d*sizeof(R)); \
179 for (WINDOW_idx = 0; WINDOW_idx < ths->d; WINDOW_idx++) \
180 ths->b[WINDOW_idx]=(K(2.0)*ths->sigma[WINDOW_idx]) / \
181 (K(2.0)*ths->sigma[WINDOW_idx] - K(1.0)) * (((R)ths->m) / KPI); \
182 }
183 #define WINDOW_HELP_FINALIZE {Y(free)(ths->b);}
184 #if MANT_DIG == 113
185 // IEEE 754 quadruple precision, 128 bits.
186 // TODO: Set good value for quadruple precision.
187 #define WINDOW_HELP_ESTIMATE_m 17
188 #elif MANT_DIG == 64
189 // Intel double extended, 80 bits.
190 #define WINDOW_HELP_ESTIMATE_m 17
191 #elif MANT_DIG == 53
192 // IEEE 754 double precision, 64 bits.
193 #define WINDOW_HELP_ESTIMATE_m 13
194 #elif MANT_DIG == 24
195 // IEEE 754 single precision, 32 bits.
196 #define WINDOW_HELP_ESTIMATE_m 5
197 #else
198 // Unknown floating-point type.
199 // Assume IEEE 754 double precision, 64 bits.
200 #define WINDOW_HELP_ESTIMATE_m 13
201 #endif
202#elif defined(B_SPLINE)
203 #define PHI_HUT(n,k,d) ((R)(((k) == 0) ? K(1.0) / n : \
204 POW(SIN((k) * KPI / n) / ((k) * KPI / n), \
205 K(2.0) * ths->m)/n))
206 #define PHI(n,x,d) (Y(bsplines)(2*ths->m,((x)*n) + \
207 (R)ths->m) / n)
208 #define WINDOW_HELP_INIT
209 #define WINDOW_HELP_FINALIZE
210 #if MANT_DIG == 113
211 // IEEE 754 quadruple precision, 128 bits.
212 // TODO: Set good value for quadruple precision.
213 #define WINDOW_HELP_ESTIMATE_m 11
214 #elif MANT_DIG == 64
215 // Intel double extended, 80 bits.
216 #define WINDOW_HELP_ESTIMATE_m 11
217 #elif MANT_DIG == 53
218 // IEEE 754 double precision, 64 bits.
219 #define WINDOW_HELP_ESTIMATE_m 11
220 #elif MANT_DIG == 24
221 // IEEE 754 single precision, 32 bits.
222 #define WINDOW_HELP_ESTIMATE_m 11
223 #else
224 // Unknown floating-point type.
225 // Assume IEEE 754 double precision, 64 bits.
226 #define WINDOW_HELP_ESTIMATE_m 11
227 #endif
228#elif defined(SINC_POWER)
229 #define PHI_HUT(n,k,d) (Y(bsplines)(2 * ths->m, (K(2.0) * ths->m*(k)) / \
230 ((K(2.0) * ths->sigma[(d)] - 1) * n / \
231 ths->sigma[(d)]) + (R)ths->m))
232 #define PHI(n,x,d) ((R)(n / ths->sigma[(d)] * \
233 (K(2.0) * ths->sigma[(d)] - K(1.0))/ (K(2.0)*ths->m) * \
234 POW(Y(sinc)(KPI * n / ths->sigma[(d)] * (x) * \
235 (K(2.0) * ths->sigma[(d)] - K(1.0)) / (K(2.0)*ths->m)) , 2*ths->m) / \
236 n))
237 #define WINDOW_HELP_INIT
238 #define WINDOW_HELP_FINALIZE
239 #if MANT_DIG == 113
240 // IEEE 754 quadruple precision, 128 bits.
241 // TODO: Set good value for quadruple precision.
242 #define WINDOW_HELP_ESTIMATE_m 13
243 #elif MANT_DIG == 64
244 // Intel double extended, 80 bits.
245 #define WINDOW_HELP_ESTIMATE_m 13
246 #elif MANT_DIG == 53
247 // IEEE 754 double precision, 64 bits.
248 #define WINDOW_HELP_ESTIMATE_m 11
249 #elif MANT_DIG == 24
250 // IEEE 754 single precision, 32 bits.
251 #define WINDOW_HELP_ESTIMATE_m 11
252 #else
253 // Unknown floating-point type.
254 // Assume IEEE 754 double precision, 64 bits.
255 #define WINDOW_HELP_ESTIMATE_m 11
256 #endif
257#else /* Kaiser-Bessel is the default. */
258 #define PHI_HUT(n,k,d) (Y(bessel_i0)((R)(ths->m) * SQRT(ths->b[d] * ths->b[d] - (K(2.0) * KPI * (R)(k) / (R)(n)) * (K(2.0) * KPI * (R)(k) / (R)(n)))))
259 #define PHI(n,x,d) ( (((R)(ths->m) * (R)(ths->m) - (x) * (R)(n) * (x) * (R)(n)) > K(0.0)) \
260 ? SINH(ths->b[d] * SQRT((R)(ths->m) * (R)(ths->m) - (x) * (R)(n) * (x) * (R)(n))) \
261 / (KPI * SQRT((R)(ths->m) * (R)(ths->m) - (x) * (R)(n) * (x) * (R)(n))) \
262 : ((((R)(ths->m) * (R)(ths->m) - (x) * (R)(n) * (x) * (R)(n)) < K(0.0)) \
263 ? SIN(ths->b[d] * SQRT((x) * (R)(n) * (x) * (R)(n) - (R)(ths->m) * (R)(ths->m))) \
264 / (KPI * SQRT((x) * (R)(n) * (x) * (R)(n) - (R)(ths->m) * (R)(ths->m))) \
265 : ths->b[d] / KPI))
266 #define WINDOW_HELP_INIT \
267 { \
268 int WINDOW_idx; \
269 ths->b = (R*) Y(malloc)((size_t)(ths->d) * sizeof(R)); \
270 for (WINDOW_idx = 0; WINDOW_idx < ths->d; WINDOW_idx++) \
271 ths->b[WINDOW_idx] = (KPI * (K(2.0) - K(1.0) / ths->sigma[WINDOW_idx])); \
272 }
273 #define WINDOW_HELP_FINALIZE {Y(free)(ths->b);}
274 #if MANT_DIG == 113
275 // IEEE 754 quadruple precision, 128 bits.
276 // TODO: Set good value for quadruple precision.
277 #define WINDOW_HELP_ESTIMATE_m 10
278 #elif MANT_DIG == 64
279 // Intel double extended, 80 bits.
280 #define WINDOW_HELP_ESTIMATE_m 9
281 #elif MANT_DIG == 53
282 // IEEE 754 double precision, 64 bits.
283 #define WINDOW_HELP_ESTIMATE_m 8
284 #elif MANT_DIG == 24
285 #define WINDOW_HELP_ESTIMATE_m 4
286 #else
287 // Unknown floating-point type.
288 // Assume IEEE 754 double precision, 64 bits.
289 #define WINDOW_HELP_ESTIMATE_m 8
290 #endif
291#endif
292
293/* window.c */
294INT Y(m2K)(const INT m);
295
296#if defined(NFFT_LDOUBLE)
297#if HAVE_DECL_COPYSIGNL == 0
298extern long double copysignl(long double, long double);
299#endif
300#if HAVE_DECL_NEXTAFTERL == 0
301extern long double nextafterl(long double, long double);
302#endif
303#if HAVE_DECL_NANL == 0
304extern long double nanl(const char *tag);
305#endif
306#if HAVE_DECL_CEILL == 0
307extern long double ceill(long double);
308#endif
309#if HAVE_DECL_FLOORL == 0
310extern long double floorl(long double);
311#endif
312#if HAVE_DECL_NEARBYINTL == 0
313extern long double nearbyintl(long double);
314#endif
315#if HAVE_DECL_RINTL == 0
316extern long double rintl(long double);
317#endif
318#if HAVE_DECL_ROUNDL == 0
319extern long double roundl(long double);
320#endif
321#if HAVE_DECL_LRINTL == 0
322extern long int lrintl(long double);
323#endif
324#if HAVE_DECL_LROUNDL == 0
325extern long int lroundl(long double);
326#endif
327#if HAVE_DECL_LLRINTL == 0
328extern long long int llrintl(long double);
329#endif
330#if HAVE_DECL_LLROUNDL == 0
331extern long long int llroundl(long double);
332#endif
333#if HAVE_DECL_TRUNCL == 0
334extern long double truncl(long double);
335#endif
336#if HAVE_DECL_FMODL == 0
337extern long double fmodl(long double, long double);
338#endif
339#if HAVE_DECL_REMAINDERL == 0
340extern long double remainderl(long double, long double);
341#endif
342#if HAVE_DECL_REMQUOL == 0
343extern long double remquol(long double x, long double y, int *);
344#endif
345#if HAVE_DECL_FDIML == 0
346extern long double fdiml(long double, long double);
347#endif
348#if HAVE_DECL_FMAXL == 0
349extern long double fmaxl(long double, long double);
350#endif
351#if HAVE_DECL_FMINL == 0
352extern long double fminl(long double, long double);
353#endif
354#if HAVE_DECL_FMAL == 0
355extern long double fmal(long double x, long double y, long double z);
356#endif
357#if HAVE_DECL_FABSL == 0
358extern long double fabsl(long double);
359#endif
360#if HAVE_DECL_SQRTL == 0
361extern long double sqrtl(long double);
362#endif
363#if HAVE_DECL_CBRTL == 0
364extern long double cbrtl(long double);
365#endif
366#if HAVE_DECL_HYPOTL == 0
367extern long double hypotl(long double, long double);
368#endif
369#if HAVE_DECL_EXPL == 0
370extern long double expl(long double);
371#endif
372#if HAVE_DECL_EXP2L == 0
373extern long double exp2l(long double);
374#endif
375#if HAVE_DECL_EXPM1L == 0
376extern long double expm1l(long double);
377#endif
378#if HAVE_DECL_LOGL == 0
379extern long double logl(long double);
380#endif
381#if HAVE_DECL_LOG2L == 0
382extern long double log2l(long double);
383#endif
384#if HAVE_DECL_LOG10L == 0
385extern long double log10l(long double);
386#endif
387#if HAVE_DECL_LOG1PL == 0
388extern long double log1pl(long double);
389#endif
390#if HAVE_DECL_LOGBL == 0
391extern long double logbl(long double);
392#endif
393#if HAVE_DECL_ILOGBL == 0
394extern int ilogbl(long double);
395#endif
396#if HAVE_DECL_MODFL == 0
397extern long double modfl(long double, long double *);
398#endif
399#if HAVE_DECL_FREXPL == 0
400extern long double frexpl(long double, int *);
401#endif
402#if HAVE_DECL_LDEXPL == 0
403extern long double ldexpl(long double, int);
404#endif
405#if HAVE_DECL_SCALBNL == 0
406extern long double scalbnl(long double, int);
407#endif
408#if HAVE_DECL_SCALBLNL == 0
409extern long double scalblnl(long double, long int);
410#endif
411#if HAVE_DECL_POWL == 0
412extern long double powl(long double, long double);
413#endif
414#if HAVE_DECL_COSL == 0
415extern long double cosl(long double);
416#endif
417#if HAVE_DECL_SINL == 0
418extern long double sinl(long double);
419#endif
420#if HAVE_DECL_TANL == 0
421extern long double tanl(long double);
422#endif
423#if HAVE_DECL_COSHL == 0
424extern long double coshl(long double);
425#endif
426#if HAVE_DECL_SINHL == 0
427extern long double sinhl(long double);
428#endif
429#if HAVE_DECL_TANHL == 0
430extern long double tanhl(long double);
431#endif
432#if HAVE_DECL_ACOSL == 0
433extern long double acosl(long double);
434#endif
435#if HAVE_DECL_ASINL == 0
436extern long double asinl(long double);
437#endif
438#if HAVE_DECL_ATANL == 0
439extern long double atanl(long double);
440#endif
441#if HAVE_DECL_ATAN2L == 0
442extern long double atan2l(long double, long double);
443#endif
444#if HAVE_DECL_ACOSHL == 0
445extern long double acoshl(long double);
446#endif
447#if HAVE_DECL_ASINHL == 0
448extern long double asinhl(long double);
449#endif
450#if HAVE_DECL_ATANHL == 0
451extern long double atanhl(long double);
452#endif
453#if HAVE_DECL_TGAMMAL == 0
454extern long double tgammal(long double);
455#endif
456#if HAVE_DECL_LGAMMAL == 0
457extern long double lgammal(long double);
458#endif
459#if HAVE_DECL_J0L == 0
460extern long double j0l(long double);
461#endif
462#if HAVE_DECL_J1L == 0
463extern long double j1l(long double);
464#endif
465#if HAVE_DECL_JNL == 0
466extern long double jnl(int, long double);
467#endif
468#if HAVE_DECL_Y0L == 0
469extern long double y0l(long double);
470#endif
471#if HAVE_DECL_Y1L == 0
472extern long double y1l(long double);
473#endif
474#if HAVE_DECL_YNL == 0
475extern long double ynl(int, long double);
476#endif
477#if HAVE_DECL_ERFL == 0
478extern long double erfl(long double);
479#endif
480#if HAVE_DECL_ERFCL == 0
481extern long double erfcl(long double);
482#endif
483#if HAVE_DECL_CREALL == 0
484extern long double creall(long double _Complex z);
485#endif
486#if HAVE_DECL_CIMAGL == 0
487extern long double cimagl(long double _Complex z);
488#endif
489#if HAVE_DECL_CABSL == 0
490extern long double cabsl(long double _Complex z);
491#endif
492#if HAVE_DECL_CARGL == 0
493extern long double cargl(long double _Complex z);
494#endif
495#if HAVE_DECL_CONJL == 0
496extern long double _Complex conjl(long double _Complex z);
497#endif
498#if HAVE_DECL_CPROJL == 0
499extern long double _Complex cprojl(long double _Complex z);
500#endif
501#if HAVE_DECL_CSQRTL == 0
502extern long double _Complex csqrtl(long double _Complex z);
503#endif
504#if HAVE_DECL_CEXPL == 0
505extern long double _Complex cexpl(long double _Complex z);
506#endif
507#if HAVE_DECL_CLOGL == 0
508extern long double _Complex clogl(long double _Complex z);
509#endif
510#if HAVE_DECL_CPOWL == 0
511extern long double _Complex cpowl(long double _Complex z, long double _Complex w);
512#endif
513#if HAVE_DECL_CSINL == 0
514extern long double _Complex csinl(long double _Complex z);
515#endif
516#if HAVE_DECL_CCOSL == 0
517extern long double _Complex ccosl(long double _Complex z);
518#endif
519#if HAVE_DECL_CTANL == 0
520extern long double _Complex ctanl(long double _Complex z);
521#endif
522#if HAVE_DECL_CASINL == 0
523extern long double _Complex casinl(long double _Complex z);
524#endif
525#if HAVE_DECL_CACOSL == 0
526extern long double _Complex cacosl(long double _Complex z);
527#endif
528#if HAVE_DECL_CATANL == 0
529extern long double _Complex catanl(long double _Complex z);
530#endif
531#if HAVE_DECL_CSINHL == 0
532extern long double _Complex csinhl(long double _Complex z);
533#endif
534#if HAVE_DECL_CCOSHL == 0
535extern long double _Complex ccoshl(long double _Complex z);
536#endif
537#if HAVE_DECL_CTANHL == 0
538extern long double _Complex ctanhl(long double _Complex z);
539#endif
540#if HAVE_DECL_CASINHL == 0
541extern long double _Complex casinhl(long double _Complex z);
542#endif
543#if HAVE_DECL_CACOSHL == 0
544extern long double _Complex cacoshl(long double _Complex z);
545#endif
546#if HAVE_DECL_CATANHL == 0
547extern long double _Complex catanhl(long double _Complex z);
548#endif
549#define COPYSIGN copysignl
550#define NEXTAFTER nextafterl
551#define MKNAN nanl
552#define CEIL ceill
553#define FLOOR floorl
554#define NEARBYINT nearbyintl
555#define RINT rintl
556#define ROUND roundl
557#define LRINT lrintl
558#define LROUND lroundl
559#define LLRINT llrintl
560#define LLROUND llroundl
561#define TRUNC truncl
562#define FMOD fmodl
563#define REMAINDER remainderl
564#define REMQUO remquol
565#define FDIM fdiml
566#define FMAX fmaxl
567#define FMIN fminl
568#define FFMA fmal
569#define FABS fabsl
570#define SQRT sqrtl
571#define CBRT cbrtl
572#define HYPOT hypotl
573#define EXP expl
574#define EXP2 exp2l
575#define EXPM1 expm1l
576#define LOG logl
577#define LOG2 log2l
578#define LOG10 log10l
579#define LOG1P log1pl
580#define LOGB logbl
581#define ILOGB ilogbl
582#define MODF modfl
583#define FREXP frexpl
584#define LDEXP ldexpl
585#define SCALBN scalbnl
586#define SCALBLN scalblnl
587#define POW powl
588#define COS cosl
589#define SIN sinl
590#define TAN tanl
591#define COSH coshl
592#define SINH sinhl
593#define TANH tanhl
594#define ACOS acosl
595#define ASIN asinl
596#define ATAN atanl
597#define ATAN2 atan2l
598#define ACOSH acoshl
599#define ASINH asinhl
600#define ATANH atanhl
601#define TGAMMA tgammal
602#define LGAMMA lgammal
603#define J0 j0l
604#define J1 j1l
605#define JN jnl
606#define Y0 y0l
607#define Y1 y1l
608#define YN ynl
609#define ERF erfl
610#define ERFC erfcl
611#define CREAL creall
612#define CIMAG cimagl
613#define CABS cabsl
614#define CARG cargl
615#define CONJ conjl
616#define CPROJ cprojl
617#define CSQRT csqrtl
618#define CEXP cexpl
619#define CLOG clogl
620#define CPOW cpowl
621#define CSIN csinl
622#define CCOS ccosl
623#define CTAN ctanl
624#define CASIN casinl
625#define CACOS cacosl
626#define CATAN catanl
627#define CSINH csinhl
628#define CCOSH ccoshl
629#define CTANH ctanhl
630#define CASINH casinhl
631#define CACOSH cacoshl
632#define CATANH catanhl
633#elif defined(NFFT_SINGLE)
634#if HAVE_DECL_COPYSIGNF == 0
635extern float copysignf(float, float);
636#endif
637#if HAVE_DECL_NEXTAFTERF == 0
638extern float nextafterf(float, float);
639#endif
640#if HAVE_DECL_NANF == 0
641extern float nanf(const char *tag);
642#endif
643#if HAVE_DECL_CEILF == 0
644extern float ceilf(float);
645#endif
646#if HAVE_DECL_FLOORF == 0
647extern float floorf(float);
648#endif
649#if HAVE_DECL_NEARBYINTF == 0
650extern float nearbyintf(float);
651#endif
652#if HAVE_DECL_RINTF == 0
653extern float rintf(float);
654#endif
655#if HAVE_DECL_ROUNDF == 0
656extern float roundf(float);
657#endif
658#if HAVE_DECL_LRINTF == 0
659extern long int lrintf(float);
660#endif
661#if HAVE_DECL_LROUNDF == 0
662extern long int lroundf(float);
663#endif
664#if HAVE_DECL_LLRINTF == 0
665extern long long int llrintf(float);
666#endif
667#if HAVE_DECL_LLROUNDF == 0
668extern long long int llroundf(float);
669#endif
670#if HAVE_DECL_TRUNCF == 0
671extern float truncf(float);
672#endif
673#if HAVE_DECL_FMODF == 0
674extern float fmodf(float, float);
675#endif
676#if HAVE_DECL_REMAINDERF == 0
677extern float remainderf(float, float);
678#endif
679#if HAVE_DECL_REMQUOF == 0
680extern float remquof(float x, float y, int *);
681#endif
682#if HAVE_DECL_FDIMF == 0
683extern float fdimf(float, float);
684#endif
685#if HAVE_DECL_FMAXF == 0
686extern float fmaxf(float, float);
687#endif
688#if HAVE_DECL_FMINF == 0
689extern float fminf(float, float);
690#endif
691#if HAVE_DECL_FMAF == 0
692extern float fmaf(float x, float y, float z);
693#endif
694#if HAVE_DECL_FABSF == 0
695extern float fabsf(float);
696#endif
697#if HAVE_DECL_SQRTF == 0
698extern float sqrtf(float);
699#endif
700#if HAVE_DECL_CBRTF == 0
701extern float cbrtf(float);
702#endif
703#if HAVE_DECL_HYPOTF == 0
704extern float hypotf(float, float);
705#endif
706#if HAVE_DECL_EXPF == 0
707extern float expf(float);
708#endif
709#if HAVE_DECL_EXP2F == 0
710extern float exp2f(float);
711#endif
712#if HAVE_DECL_EXPM1F == 0
713extern float expm1f(float);
714#endif
715#if HAVE_DECL_LOGF == 0
716extern float logf(float);
717#endif
718#if HAVE_DECL_LOG2F == 0
719extern float log2f(float);
720#endif
721#if HAVE_DECL_LOG10F == 0
722extern float log10f(float);
723#endif
724#if HAVE_DECL_LOG1PF == 0
725extern float log1pf(float);
726#endif
727#if HAVE_DECL_LOGBF == 0
728extern float logbf(float);
729#endif
730#if HAVE_DECL_ILOGBF == 0
731extern int ilogbf(float);
732#endif
733#if HAVE_DECL_MODFF == 0
734extern float modff(float, float *);
735#endif
736#if HAVE_DECL_FREXPF == 0
737extern float frexpf(float, int *);
738#endif
739#if HAVE_DECL_LDEXPF == 0
740extern float ldexpf(float, int);
741#endif
742#if HAVE_DECL_SCALBNF == 0
743extern float scalbnf(float, int);
744#endif
745#if HAVE_DECL_SCALBLNF == 0
746extern float scalblnf(float, long int);
747#endif
748#if HAVE_DECL_POWF == 0
749extern float powf(float, float);
750#endif
751#if HAVE_DECL_COSF == 0
752extern float cosf(float);
753#endif
754#if HAVE_DECL_SINF == 0
755extern float sinf(float);
756#endif
757#if HAVE_DECL_TANF == 0
758extern float tanf(float);
759#endif
760#if HAVE_DECL_COSHF == 0
761extern float coshf(float);
762#endif
763#if HAVE_DECL_SINHF == 0
764extern float sinhf(float);
765#endif
766#if HAVE_DECL_TANHF == 0
767extern float tanhf(float);
768#endif
769#if HAVE_DECL_ACOSF == 0
770extern float acosf(float);
771#endif
772#if HAVE_DECL_ASINF == 0
773extern float asinf(float);
774#endif
775#if HAVE_DECL_ATANF == 0
776extern float atanf(float);
777#endif
778#if HAVE_DECL_ATAN2F == 0
779extern float atan2f(float, float);
780#endif
781#if HAVE_DECL_ACOSHF == 0
782extern float acoshf(float);
783#endif
784#if HAVE_DECL_ASINHF == 0
785extern float asinhf(float);
786#endif
787#if HAVE_DECL_ATANHF == 0
788extern float atanhf(float);
789#endif
790#if HAVE_DECL_TGAMMAF == 0
791extern float tgammaf(float);
792#endif
793#if HAVE_DECL_LGAMMAF == 0
794extern float lgammaf(float);
795#endif
796#if HAVE_DECL_J0F == 0
797extern float j0f(float);
798#endif
799#if HAVE_DECL_J1F == 0
800extern float j1f(float);
801#endif
802#if HAVE_DECL_JNF == 0
803extern float jnf(int, float);
804#endif
805#if HAVE_DECL_Y0F == 0
806extern float y0f(float);
807#endif
808#if HAVE_DECL_Y1F == 0
809extern float y1f(float);
810#endif
811#if HAVE_DECL_YNF == 0
812extern float ynf(int, float);
813#endif
814#if HAVE_DECL_ERFF == 0
815extern float erff(float);
816#endif
817#if HAVE_DECL_ERFCF == 0
818extern float erfcf(float);
819#endif
820#if HAVE_DECL_CREALF == 0
821extern float crealf(float _Complex z);
822#endif
823#if HAVE_DECL_CIMAGF == 0
824extern float cimagf(float _Complex z);
825#endif
826#if HAVE_DECL_CABSF == 0
827extern float cabsf(float _Complex z);
828#endif
829#if HAVE_DECL_CARGF == 0
830extern float cargf(float _Complex z);
831#endif
832#if HAVE_DECL_CONJF == 0
833extern float _Complex conjf(float _Complex z);
834#endif
835#if HAVE_DECL_CPROJF == 0
836extern float _Complex cprojf(float _Complex z);
837#endif
838#if HAVE_DECL_CSQRTF == 0
839extern float _Complex csqrtf(float _Complex z);
840#endif
841#if HAVE_DECL_CEXPF == 0
842extern float _Complex cexpf(float _Complex z);
843#endif
844#if HAVE_DECL_CLOGF == 0
845extern float _Complex clogf(float _Complex z);
846#endif
847#if HAVE_DECL_CPOWF == 0
848extern float _Complex cpowf(float _Complex z, float _Complex w);
849#endif
850#if HAVE_DECL_CSINF == 0
851extern float _Complex csinf(float _Complex z);
852#endif
853#if HAVE_DECL_CCOSF == 0
854extern float _Complex ccosf(float _Complex z);
855#endif
856#if HAVE_DECL_CTANF == 0
857extern float _Complex ctanf(float _Complex z);
858#endif
859#if HAVE_DECL_CASINF == 0
860extern float _Complex casinf(float _Complex z);
861#endif
862#if HAVE_DECL_CACOSF == 0
863extern float _Complex cacosf(float _Complex z);
864#endif
865#if HAVE_DECL_CATANF == 0
866extern float _Complex catanf(float _Complex z);
867#endif
868#if HAVE_DECL_CSINHF == 0
869extern float _Complex csinhf(float _Complex z);
870#endif
871#if HAVE_DECL_CCOSHF == 0
872extern float _Complex ccoshf(float _Complex z);
873#endif
874#if HAVE_DECL_CTANHF == 0
875extern float _Complex ctanhf(float _Complex z);
876#endif
877#if HAVE_DECL_CASINHF == 0
878extern float _Complex casinhf(float _Complex z);
879#endif
880#if HAVE_DECL_CACOSHF == 0
881extern float _Complex cacoshf(float _Complex z);
882#endif
883#if HAVE_DECL_CATANHF == 0
884extern float _Complex catanhf(float _Complex z);
885#endif
886#define COPYSIGN copysignf
887#define NEXTAFTER nextafterf
888#define MKNAN nanf
889#define CEIL ceilf
890#define FLOOR floorf
891#define NEARBYINT nearbyintf
892#define RINT rintf
893#define ROUND roundf
894#define LRINT lrintf
895#define LROUND lroundf
896#define LLRINT llrintf
897#define LLROUND llroundf
898#define TRUNC truncf
899#define FMOD fmodf
900#define REMAINDER remainderf
901#define REMQUO remquof
902#define FDIM fdimf
903#define FMAX fmaxf
904#define FMIN fminf
905#define FFMA fmaf
906#define FABS fabsf
907#define SQRT sqrtf
908#define CBRT cbrtf
909#define HYPOT hypotf
910#define EXP expf
911#define EXP2 exp2f
912#define EXPM1 expm1f
913#define LOG logf
914#define LOG2 log2f
915#define LOG10 log10f
916#define LOG1P log1pf
917#define LOGB logbf
918#define ILOGB ilogbf
919#define MODF modff
920#define FREXP frexpf
921#define LDEXP ldexpf
922#define SCALBN scalbnf
923#define SCALBLN scalblnf
924#define POW powf
925#define COS cosf
926#define SIN sinf
927#define TAN tanf
928#define COSH coshf
929#define SINH sinhf
930#define TANH tanhf
931#define ACOS acosf
932#define ASIN asinf
933#define ATAN atanf
934#define ATAN2 atan2f
935#define ACOSH acoshf
936#define ASINH asinhf
937#define ATANH atanhf
938#define TGAMMA tgammaf
939#define LGAMMA lgammaf
940#define J0 j0f
941#define J1 j1f
942#define JN jnf
943#define Y0 y0f
944#define Y1 y1f
945#define YN ynf
946#define ERF erff
947#define ERFC erfcf
948#define CREAL crealf
949#define CIMAG cimagf
950#define CABS cabsf
951#define CARG cargf
952#define CONJ conjf
953#define CPROJ cprojf
954#define CSQRT csqrtf
955#define CEXP cexpf
956#define CLOG clogf
957#define CPOW cpowf
958#define CSIN csinf
959#define CCOS ccosf
960#define CTAN ctanf
961#define CASIN casinf
962#define CACOS cacosf
963#define CATAN catanf
964#define CSINH csinhf
965#define CCOSH ccoshf
966#define CTANH ctanhf
967#define CASINH casinhf
968#define CACOSH cacoshf
969#define CATANH catanhf
970#else
971#if HAVE_DECL_COPYSIGN == 0
972extern double copysign(double, double);
973#endif
974#if HAVE_DECL_NEXTAFTER == 0
975extern double nextafter(double, double);
976#endif
977#if HAVE_DECL_NAN == 0
978extern double nan(const char *tag);
979#endif
980#if HAVE_DECL_CEIL == 0
981extern double ceil(double);
982#endif
983#if HAVE_DECL_FLOOR == 0
984extern double floor(double);
985#endif
986#if HAVE_DECL_NEARBYINT == 0
987extern double nearbyint(double);
988#endif
989#if HAVE_DECL_RINT == 0
990extern double rint(double);
991#endif
992#if HAVE_DECL_ROUND == 0
993extern double round(double);
994#endif
995#if HAVE_DECL_LRINT == 0
996extern long int lrint(double);
997#endif
998#if HAVE_DECL_LROUND == 0
999extern long int lround(double);
1000#endif
1001#if HAVE_DECL_LLRINT == 0
1002extern long long int llrint(double);
1003#endif
1004#if HAVE_DECL_LLROUND == 0
1005extern long long int llround(double);
1006#endif
1007#if HAVE_DECL_TRUNC == 0
1008extern double trunc(double);
1009#endif
1010#if HAVE_DECL_FMOD == 0
1011extern double fmod(double, double);
1012#endif
1013#if HAVE_DECL_REMAINDER == 0
1014extern double remainder(double, double);
1015#endif
1016#if HAVE_DECL_REMQUO == 0
1017extern double remquo(double x, double y, int *);
1018#endif
1019#if HAVE_DECL_FDIM == 0
1020extern double fdim(double, double);
1021#endif
1022#if HAVE_DECL_FMAX == 0
1023extern double fmax(double, double);
1024#endif
1025#if HAVE_DECL_FMIN == 0
1026extern double fmin(double, double);
1027#endif
1028#if HAVE_DECL_FMA == 0
1029extern double fma(double x, double y, double z);
1030#endif
1031#if HAVE_DECL_FABS == 0
1032extern double fabs(double);
1033#endif
1034#if HAVE_DECL_SQRT == 0
1035extern double sqrt(double);
1036#endif
1037#if HAVE_DECL_CBRT == 0
1038extern double cbrt(double);
1039#endif
1040#if HAVE_DECL_HYPOT == 0
1041extern double hypot(double, double);
1042#endif
1043#if HAVE_DECL_EXP == 0
1044extern double exp(double);
1045#endif
1046#if HAVE_DECL_EXP2 == 0
1047extern double exp2(double);
1048#endif
1049#if HAVE_DECL_EXPM1 == 0
1050extern double expm1(double);
1051#endif
1052#if HAVE_DECL_LOG == 0
1053extern double log(double);
1054#endif
1055#if HAVE_DECL_LOG2 == 0
1056extern double log2(double);
1057#endif
1058#if HAVE_DECL_LOG10 == 0
1059extern double log10(double);
1060#endif
1061#if HAVE_DECL_LOG1P == 0
1062extern double log1p(double);
1063#endif
1064#if HAVE_DECL_LOGB == 0
1065extern double logb(double);
1066#endif
1067#if HAVE_DECL_ILOGB == 0
1068extern int ilogb(double);
1069#endif
1070#if HAVE_DECL_MODF == 0
1071extern double modf(double, double *);
1072#endif
1073#if HAVE_DECL_FREXP == 0
1074extern double frexp(double, int *);
1075#endif
1076#if HAVE_DECL_LDEXP == 0
1077extern double ldexp(double, int);
1078#endif
1079#if HAVE_DECL_SCALBN == 0
1080extern double scalbn(double, int);
1081#endif
1082#if HAVE_DECL_SCALBLN == 0
1083extern double scalbln(double, long int);
1084#endif
1085#if HAVE_DECL_POW == 0
1086extern double pow(double, double);
1087#endif
1088#if HAVE_DECL_COS == 0
1089extern double cos(double);
1090#endif
1091#if HAVE_DECL_SIN == 0
1092extern double sin(double);
1093#endif
1094#if HAVE_DECL_TAN == 0
1095extern double tan(double);
1096#endif
1097#if HAVE_DECL_COSH == 0
1098extern double cosh(double);
1099#endif
1100#if HAVE_DECL_SINH == 0
1101extern double sinh(double);
1102#endif
1103#if HAVE_DECL_TANH == 0
1104extern double tanh(double);
1105#endif
1106#if HAVE_DECL_ACOS == 0
1107extern double acos(double);
1108#endif
1109#if HAVE_DECL_ASIN == 0
1110extern double asin(double);
1111#endif
1112#if HAVE_DECL_ATAN == 0
1113extern double atan(double);
1114#endif
1115#if HAVE_DECL_ATAN2 == 0
1116extern double atan2(double, double);
1117#endif
1118#if HAVE_DECL_ACOSH == 0
1119extern double acosh(double);
1120#endif
1121#if HAVE_DECL_ASINH == 0
1122extern double asinh(double);
1123#endif
1124#if HAVE_DECL_ATANH == 0
1125extern double atanh(double);
1126#endif
1127#if HAVE_DECL_TGAMMA == 0
1128extern double tgamma(double);
1129#endif
1130#if HAVE_DECL_LGAMMA == 0
1131extern double lgamma(double);
1132#endif
1133#if HAVE_DECL_J0 == 0
1134extern double j0(double);
1135#endif
1136#if HAVE_DECL_J1 == 0
1137extern double j1(double);
1138#endif
1139#if HAVE_DECL_JN == 0
1140extern double jn(int, double);
1141#endif
1142#if HAVE_DECL_Y0 == 0
1143extern double y0(double);
1144#endif
1145#if HAVE_DECL_Y1 == 0
1146extern double y1(double);
1147#endif
1148#if HAVE_DECL_YN == 0
1149extern double yn(int, double);
1150#endif
1151#if HAVE_DECL_ERF == 0
1152extern double erf(double);
1153#endif
1154#if HAVE_DECL_ERFC == 0
1155extern double erfc(double);
1156#endif
1157#if HAVE_DECL_CREAL == 0
1158extern double creal(double _Complex z);
1159#endif
1160#if HAVE_DECL_CIMAG == 0
1161extern double cimag(double _Complex z);
1162#endif
1163#if HAVE_DECL_CABS == 0
1164extern double cabs(double _Complex z);
1165#endif
1166#if HAVE_DECL_CARG == 0
1167extern double carg(double _Complex z);
1168#endif
1169#if HAVE_DECL_CONJ == 0
1170extern double _Complex conj(double _Complex z);
1171#endif
1172#if HAVE_DECL_CPROJ == 0
1173extern double _Complex cproj(double _Complex z);
1174#endif
1175#if HAVE_DECL_CSQRT == 0
1176extern double _Complex csqrt(double _Complex z);
1177#endif
1178#if HAVE_DECL_CEXP == 0
1179extern double _Complex cexp(double _Complex z);
1180#endif
1181#if HAVE_DECL_CLOG == 0
1182extern double _Complex clog(double _Complex z);
1183#endif
1184#if HAVE_DECL_CPOW == 0
1185extern double _Complex cpow(double _Complex z, double _Complex w);
1186#endif
1187#if HAVE_DECL_CSIN == 0
1188extern double _Complex csin(double _Complex z);
1189#endif
1190#if HAVE_DECL_CCOS == 0
1191extern double _Complex ccos(double _Complex z);
1192#endif
1193#if HAVE_DECL_CTAN == 0
1194extern double _Complex ctan(double _Complex z);
1195#endif
1196#if HAVE_DECL_CASIN == 0
1197extern double _Complex casin(double _Complex z);
1198#endif
1199#if HAVE_DECL_CACOS == 0
1200extern double _Complex cacos(double _Complex z);
1201#endif
1202#if HAVE_DECL_CATAN == 0
1203extern double _Complex catan(double _Complex z);
1204#endif
1205#if HAVE_DECL_CSINH == 0
1206extern double _Complex csinh(double _Complex z);
1207#endif
1208#if HAVE_DECL_CCOSH == 0
1209extern double _Complex ccosh(double _Complex z);
1210#endif
1211#if HAVE_DECL_CTANH == 0
1212extern double _Complex ctanh(double _Complex z);
1213#endif
1214#if HAVE_DECL_CASINH == 0
1215extern double _Complex casinh(double _Complex z);
1216#endif
1217#if HAVE_DECL_CACOSH == 0
1218extern double _Complex cacosh(double _Complex z);
1219#endif
1220#if HAVE_DECL_CATANH == 0
1221extern double _Complex catanh(double _Complex z);
1222#endif
1223#define COPYSIGN copysign
1224#define NEXTAFTER nextafter
1225#define MKNAN nan
1226#define CEIL ceil
1227#define FLOOR floor
1228#define NEARBYINT nearbyint
1229#define RINT rint
1230#define ROUND round
1231#define LRINT lrint
1232#define LROUND lround
1233#define LLRINT llrint
1234#define LLROUND llround
1235#define TRUNC trunc
1236#define FMOD fmod
1237#define REMAINDER remainder
1238#define REMQUO remquo
1239#define FDIM fdim
1240#define FMAX fmax
1241#define FMIN fmin
1242#define FFMA fma
1243#define FABS fabs
1244#define SQRT sqrt
1245#define CBRT cbrt
1246#define HYPOT hypot
1247#define EXP exp
1248#define EXP2 exp2
1249#define EXPM1 expm1
1250#define LOG log
1251#define LOG2 log2
1252#define LOG10 log10
1253#define LOG1P log1p
1254#define LOGB logb
1255#define ILOGB ilogb
1256#define MODF modf
1257#define FREXP frexp
1258#define LDEXP ldexp
1259#define SCALBN scalbn
1260#define SCALBLN scalbln
1261#define POW pow
1262#define COS cos
1263#define SIN sin
1264#define TAN tan
1265#define COSH cosh
1266#define SINH sinh
1267#define TANH tanh
1268#define ACOS acos
1269#define ASIN asin
1270#define ATAN atan
1271#define ATAN2 atan2
1272#define ACOSH acosh
1273#define ASINH asinh
1274#define ATANH atanh
1275#define TGAMMA tgamma
1276#define LGAMMA lgamma
1277#define J0 j0
1278#define J1 j1
1279#define JN jn
1280#define Y0 y0
1281#define Y1 y1
1282#define YN yn
1283#define ERF erf
1284#define ERFC erfc
1285#define CREAL creal
1286#define CIMAG cimag
1287#define CABS cabs
1288#define CARG carg
1289#define CONJ conj
1290#define CPROJ cproj
1291#define CSQRT csqrt
1292#define CEXP cexp
1293#define CLOG clog
1294#define CPOW cpow
1295#define CSIN csin
1296#define CCOS ccos
1297#define CTAN ctan
1298#define CASIN casin
1299#define CACOS cacos
1300#define CATAN catan
1301#define CSINH csinh
1302#define CCOSH ccosh
1303#define CTANH ctanh
1304#define CASINH casinh
1305#define CACOSH cacosh
1306#define CATANH catanh
1307#endif
1308
1309#if defined(FLT_ROUND)
1310 #if FLT_ROUND != -1
1311 #define FLTROUND 1.0
1312 #else
1313 #define FLTROUND 0.0
1314 #endif
1315#else
1316 #define FLTROUND 0.0
1317#endif
1318
1319#if HAVE_DECL_DRAND48 == 0
1320 extern double drand48(void);
1321#endif
1322#if HAVE_DECL_SRAND48 == 0
1323 extern void srand48(long int);
1324#endif
1325#define R_RADIX FLT_RADIX
1326#define II _Complex_I
1327
1328/* format strings */
1329#if defined(NFFT_LDOUBLE)
1330# define __FGS__ "Lg"
1331# define __FES__ "LE"
1332# define __FI__ "%Lf"
1333# define __FIS__ "Lf"
1334# define __FR__ "%Le"
1335#elif defined(NFFT_SINGLE)
1336# define __FGS__ "g"
1337# define __FES__ "E"
1338# define __FI__ "%f"
1339# define __FIS__ "f"
1340# define __FR__ "%e"
1341#else // double precision
1342# define __FGS__ "lg"
1343# define __FES__ "lE"
1344# define __FI__ "%lf"
1345# define __FIS__ "lf"
1346# define __FR__ "%le"
1347#endif
1348
1349#if MANT_DIG == 113
1350 // IEEE 754 quadruple precision, 128 bits.
1351 #define __FE__ "% 36.32LE"
1352#elif MANT_DIG == 64
1353 // Intel double extended, 80 bits.
1354 #define __FE__ "% 24.20LE"
1355#elif MANT_DIG == 53
1356 // IEEE 754 double precision, 64 bits.
1357 #define __FE__ "% 20.16lE"
1358#elif MANT_DIG == 24
1359 #define __FE__ "% 12.8E"
1360#else
1361 // Unknown floating-point type.
1362 // Assume IEEE 754 double precision, 64 bits.
1363 #define __FE__ "% 20.16LE"
1364#endif
1365
1366#define TRUE 1
1367#define FALSE 0
1368
1369#if defined(_WIN32) || defined(_WIN64)
1370# define __D__ "%Id"
1371#else
1372# define __D__ "%td"
1373#endif
1374
1376#define UNUSED(x) (void)x
1377
1378#ifdef HAVE_ALLOCA
1379 /* Use alloca if available. */
1380 #ifndef alloca
1381 #ifdef __GNUC__
1382 /* No alloca defined but can use GCC's builtin version. */
1383 #define alloca __builtin_alloca
1384 #else
1385 /* No alloca defined and not using GCC. */
1386 #ifdef _MSC_VER
1387 /* Using Microsoft's C compiler. Include header file and use _alloca
1388 * defined therein. */
1389 #include <malloc.h>
1390 #define alloca _alloca
1391 #else
1392 /* Also not using Microsoft's C compiler. */
1393 #if HAVE_ALLOCA_H
1394 /* Alloca header is available. */
1395 #include <alloca.h>
1396 #else
1397 /* No alloca header available. */
1398 #ifdef _AIX
1399 /* We're using the AIX C compiler. Use pragma. */
1400 #pragma alloca
1401 #else
1402 /* Not using AIX compiler. */
1403 #ifndef alloca /* HP's cc +Olibcalls predefines alloca. */
1404 void *alloca(size_t);
1405 #endif
1406 #endif
1407 #endif
1408 #endif
1409 #endif
1410 #endif
1411 /* So we have alloca. */
1412 #define STACK_MALLOC(T, p, x) p = (T)alloca(x)
1413 #define STACK_FREE(x) /* Nothing. Cleanup done automatically. */
1414#else /* ! HAVE_ALLOCA */
1415 /* Use malloc instead of alloca. So we allocate memory on the heap instead of
1416 * on the stack which is slower. */
1417 #define STACK_MALLOC(T, p, x) p = (T)Y(malloc)(x)
1418 #define STACK_FREE(x) Y(free)(x)
1419#endif /* ! HAVE_ALLOCA */
1420
1422R Y(elapsed_seconds)(ticks t1, ticks t0);
1423
1425#define UNUSED(x) (void)x
1426
1433#ifdef MEASURE_TIME
1434 int MEASURE_TIME_r;
1435 double MEASURE_TIME_tt;
1436 ticks MEASURE_TIME_t0, MEASURE_TIME_t1;
1437
1438#define TIC(a) \
1439 ths->MEASURE_TIME_t[(a)]=0; \
1440 MEASURE_TIME_r=0; \
1441 /* DISABLED LOOP due to code blocks causing segfault when repeatedly run */ \
1442 /*while(ths->MEASURE_TIME_t[(a)]<0.01)*/ \
1443 { \
1444 MEASURE_TIME_r++; \
1445 MEASURE_TIME_t0 = getticks(); \
1446
1447/* THE MEASURED FUNCTION IS CALLED REPEATEDLY */
1448
1449#define TOC(a) \
1450 MEASURE_TIME_t1 = getticks(); \
1451 MEASURE_TIME_tt = Y(elapsed_seconds)(MEASURE_TIME_t1,MEASURE_TIME_t0);\
1452 ths->MEASURE_TIME_t[(a)]+=MEASURE_TIME_tt; \
1453 } \
1454 ths->MEASURE_TIME_t[(a)]/=MEASURE_TIME_r; \
1455
1456#else
1457#define TIC(a)
1458#define TOC(a)
1459#endif
1460
1461#ifdef MEASURE_TIME_FFTW
1462#define TIC_FFTW(a) TIC(a)
1463#define TOC_FFTW(a) TOC(a)
1464#else
1465#define TIC_FFTW(a)
1466#define TOC_FFTW(a)
1467#endif
1468
1469/* sinc.c: */
1470
1471/* Sinus cardinalis. */
1472R Y(sinc)(R x);
1473
1474/* lambda.c: */
1475
1476/* lambda(z, eps) = gamma(z + eps) / gamma(z + 1) */
1477R Y(lambda)(R z, R eps);
1478
1479/* lambda2(mu, nu) = sqrt(gamma(mu + nu + 1) / (gamma(mu + 1) * gamma(nu + 1))) */
1480R Y(lambda2)(R mu, R nu);
1481
1482/* bessel_i0.c: */
1483R Y(bessel_i0)(R x);
1484
1485/* bspline.c: */
1486R Y(bsplines)(const INT, const R x);
1487
1488/* float.c: */
1489typedef enum {NFFT_EPSILON = 0, NFFT_SAFE__MIN = 1, NFFT_BASE = 2,
1490 NFFT_PRECISION = 3, NFFT_MANT_DIG = 4, NFFT_FLTROUND = 5, NFFT_E_MIN = 6,
1491 NFFT_R_MIN = 7, NFFT_E_MAX = 8, NFFT_R_MAX = 9 } float_property;
1492
1493R Y(float_property)(float_property);
1494R Y(prod_real)(R *vec, INT d);
1495
1496/* int.c: */
1497INT Y(log2i)(const INT m);
1498void Y(next_power_of_2_exp)(const INT N, INT *N2, INT *t);
1499void Y(next_power_of_2_exp_int)(const int N, int *N2, int *t);
1500
1501/* error.c: */
1502/* not used */ R Y(error_l_infty_double)(const R *x, const R *y, const INT n);
1503/* not used */ R Y(error_l_infty_1_double)(const R *x, const R *y, const INT n, const R *z,
1504 const INT m);
1505R Y(error_l_2_complex)(const C *x, const C *y, const INT n);
1506/* not used */ R Y(error_l_2_double)(const R *x, const R *y, const INT n);
1507
1508/* sort.c: */
1509void Y(sort_node_indices_radix_msdf)(INT n, INT *keys0, INT *keys1, INT rhigh);
1510void Y(sort_node_indices_radix_lsdf)(INT n, INT *keys0, INT *keys1, INT rhigh);
1511
1512/* assert.c */
1513void Y(assertion_failed)(const char *s, int line, const char *file);
1514
1515/* vector1.c */
1517R Y(dot_double)(R *x, INT n);
1519R Y(dot_w_complex)(C *x, R *w, INT n);
1521R Y(dot_w_double)(R *x, R *w, INT n);
1523R Y(dot_w_w2_complex)(C *x, R *w, R *w2, INT n);
1525R Y(dot_w2_complex)(C *x, R *w2, INT n);
1526
1527/* vector2.c */
1529void Y(cp_complex)(C *x, C *y, INT n);
1531void Y(cp_double)(R *x, R *y, INT n);
1533void Y(cp_a_complex)(C *x, R a, C *y, INT n);
1535void Y(cp_a_double)(R *x, R a, R *y, INT n);
1537void Y(cp_w_complex)(C *x, R *w, C *y, INT n);
1539void Y(cp_w_double)(R *x, R *w, R *y, INT n);
1540
1541/* vector3.c */
1543void Y(upd_axpy_double)(R *x, R a, R *y, INT n);
1545void Y(upd_xpay_complex)(C *x, R a, C *y, INT n);
1547void Y(upd_xpay_double)(R *x, R a, R *y, INT n);
1549void Y(upd_axpby_complex)(C *x, R a, C *y, R b, INT n);
1551void Y(upd_axpby_double)(R *x, R a, R *y, R b, INT n);
1553void Y(upd_xpawy_complex)(C *x, R a, R *w, C *y, INT n);
1555void Y(upd_xpawy_double)(R *x, R a, R *w, R *y, INT n);
1557void Y(upd_axpwy_complex)(C *x, R a, R *w, C *y, INT n);
1559void Y(upd_axpwy_double)(R *x, R a, R *w, R *y, INT n);
1560
1561/* voronoi.c */
1562void Y(voronoi_weights_1d)(R *w, R *x, const INT M);
1563
1564/* damp.c */
1569R Y(modified_fejer)(const INT N, const INT kk);
1571R Y(modified_jackson2)(const INT N, const INT kk);
1573R Y(modified_jackson4)(const INT N, const INT kk);
1575R Y(modified_sobolev)(const R mu, const INT kk);
1577R Y(modified_multiquadric)(const R mu, const R c, const INT kk);
1578
1579/* always check */
1580#define CK(ex) \
1581 (void)((ex) || (Y(assertion_failed)(#ex, __LINE__, __FILE__), 0))
1582
1583#ifdef NFFT_DEBUG
1584 /* check only if debug enabled */
1585 #define A(ex) \
1586 (void)((ex) || (Y(assertion_failed)(#ex, __LINE__, __FILE__), 0))
1587#else
1588 #define A(ex) /* nothing */
1589#endif
1590
1594#endif