De-macroize step 1/n: Probably uncontroversial switches from macros to constexpr

include/cufinufft/common.h

@@ -41,8 +41,8 @@ template<typename T>
 void cufinufft_setup_binsize(int type, int ns, int dim, cufinufft_opts *opts);
 
 template<typename T, typename V>
-auto cufinufft_set_shared_memory(V *kernel, const int dim,
-                                 const cufinufft_plan_t<T> &d_plan) {
+int cufinufft_set_shared_memory(V *kernel, const int dim,
+                                const cufinufft_plan_t<T> &d_plan) {
   /**
    * WARNING: this function does not handle cuda errors. The caller should check them.
    */

include/finufft.h

@@ -35,10 +35,6 @@
 #include <stdint.h>
 #define FINUFFT_BIGINT int64_t
 
-#ifndef __cplusplus
-#include <stdbool.h> // for bool type in C (needed for item in plan struct)
-#endif
-
 // this macro name has to be safe since exposed to user
 #define FINUFFT_SINGLE
 #include <finufft_eitherprec.h>

include/finufft/defs.h

@@ -24,18 +24,16 @@
 
 // All indexing in library that potentially can exceed 2^31 uses 64-bit signed.
 // This includes all calling arguments (eg M,N) that could be huge someday.
-#define BIGINT  int64_t
-#define UBIGINT uint64_t
+using BIGINT  = int64_t;
+using UBIGINT = uint64_t;
 // Precision-independent real and complex types, for private lib/test compile
 #ifdef SINGLE
-#define FLT float
+using FLT = float;
 #else
-#define FLT double
+using FLT = double;
 #endif
-// next line possibly obsolete...
-#define _USE_MATH_DEFINES
 #include <complex> // we define C++ complex type only
-#define CPX std::complex<FLT>
+using CPX = std::complex<FLT>;
 
 // inline macro, to force inlining of small functions
 // this avoids the use of macros to implement functions
@@ -65,44 +63,49 @@
 // ------------- Library-wide algorithm parameter settings ----------------
 
 // Library version (is a string)
-#define FINUFFT_VER          "2.3.0"
+#define FINUFFT_VER "2.3.0"
 
 // Smallest possible kernel spread width per dimension, in fine grid points
 // (used only in spreadinterp.cpp)
-#define MIN_NSPREAD          2
+inline constexpr int MIN_NSPREAD = 2;
 
 // Largest possible kernel spread width per dimension, in fine grid points
 // (used only in spreadinterp.cpp)
-#define MAX_NSPREAD          16
+inline constexpr int MAX_NSPREAD = 16;
 
 // Fraction growth cut-off in utils:arraywidcen, sets when translate in type-3
-#define ARRAYWIDCEN_GROWFRAC 0.1
+inline constexpr double ARRAYWIDCEN_GROWFRAC = 0.1;
 
 // Max number of positive quadr nodes for kernel FT (used only in common.cpp)
-#define MAX_NQUAD            100
+inline constexpr int MAX_NQUAD = 100;
 
 // Internal (nf1 etc) array allocation size that immediately raises error.
 // (Note: next235 takes 1s for 1e11, so it is also to prevent hang here.)
 // Increase this if you need >10TB (!) RAM...
-#define MAX_NF               (BIGINT)1e12
+inline constexpr BIGINT MAX_NF = BIGINT(1e12);
 
 // Maximum allowed number M of NU points; useful to catch incorrectly cast int32
 // values for M = nj (also nk in type 3)...
-#define MAX_NU_PTS           (BIGINT)1e14
+inline constexpr BIGINT MAX_NU_PTS = BIGINT(1e14);
 
 // -------------- Math consts (not in math.h) and useful math macros ----------
 #include <math.h>
 
 // either-precision unit imaginary number...
 #define IMA (CPX(0.0, 1.0))
-// using namespace std::complex_literals;  // needs C++14, provides 1i, 1if
+
+// MR: In the longer term I suggest to move
+// away from M_PI, which was never part of the standard.
+// Perhaps a constexpr pi in the namespace finufft, or a constexpr finufft_pi
+// if no namespaces are used?
+// In C++20 these constants will be part of the language, and the problem will go away.
 #ifndef M_PI // Windows apparently doesn't have this const
 #define M_PI 3.14159265358979329
 #endif
 #define M_1_2PI 0.159154943091895336
 #define M_2PI   6.28318530717958648
 // to avoid mixed precision operators in eg i*pi, an either-prec PI...
-#define PI      (FLT) M_PI
+#define PI      FLT(M_PI)
 
 // machine epsilon for decisions of achievable tolerance...
 #ifdef SINGLE
@@ -115,36 +118,45 @@
 // These macros should probably be replaced by modern C++ std lib or random123.
 // (RAND_MAX is in stdlib.h)
 #include <stdlib.h>
-// #define rand01() (((FLT)(rand()%RAND_MAX))/RAND_MAX)
-#define rand01()     ((FLT)rand() / (FLT)RAND_MAX)
+static inline FLT rand01() { return FLT(rand()) / FLT(RAND_MAX); }
 // unif[-1,1]:
-#define randm11()    (2 * rand01() - (FLT)1.0)
+static inline FLT randm11() { return 2 * rand01() - FLT(1); }
 // complex unif[-1,1] for Re and Im:
-#define crandm11()   (randm11() + IMA * randm11())
+static inline CPX crandm11() { return randm11() + IMA * randm11(); }
 
 // Thread-safe seed-carrying versions of above (x is ptr to seed)...
+// MR: we have to leave those as macros for now, as "rand_r" is deprecated
+// and apparently no longer available on Windows.
+#if 1
 #define rand01r(x)   ((FLT)rand_r(x) / (FLT)RAND_MAX)
 // unif[-1,1]:
 #define randm11r(x)  (2 * rand01r(x) - (FLT)1.0)
 // complex unif[-1,1] for Re and Im:
 #define crandm11r(x) (randm11r(x) + IMA * randm11r(x))
+#else
+static inline FLT rand01r(unsigned int *x) { return FLT(rand_r(x)) / FLT(RAND_MAX); }
+// unif[-1,1]:
+static inline FLT randm11r(unsigned int *x) { return 2 * rand01r(x) - FLT(1); }
+// complex unif[-1,1] for Re and Im:
+static inline CPX crandm11r(unsigned int *x) { return randm11r(x) + IMA * randm11r(x); }
+#endif
 
 // ----- OpenMP macros which also work when omp not present -----
 // Allows compile-time switch off of openmp, so compilation without any openmp
 // is done (Note: _OPENMP is automatically set by -fopenmp compile flag)
 #ifdef _OPENMP
 #include <omp.h>
 // point to actual omp utils
-#define MY_OMP_GET_NUM_THREADS()  omp_get_num_threads()
-#define MY_OMP_GET_MAX_THREADS()  omp_get_max_threads()
-#define MY_OMP_GET_THREAD_NUM()   omp_get_thread_num()
-#define MY_OMP_SET_NUM_THREADS(x) omp_set_num_threads(x)
+static inline int MY_OMP_GET_NUM_THREADS() { return omp_get_num_threads(); }
+static inline int MY_OMP_GET_MAX_THREADS() { return omp_get_max_threads(); }
+static inline int MY_OMP_GET_THREAD_NUM() { return omp_get_thread_num(); }
+static inline void MY_OMP_SET_NUM_THREADS(int x) { omp_set_num_threads(x); }
 #else
 // non-omp safe dummy versions of omp utils...
-#define MY_OMP_GET_NUM_THREADS() 1
-#define MY_OMP_GET_MAX_THREADS() 1
-#define MY_OMP_GET_THREAD_NUM()  0
-#define MY_OMP_SET_NUM_THREADS(x)
+static inline int MY_OMP_GET_NUM_THREADS() { return 1; }
+static inline int MY_OMP_GET_MAX_THREADS() { return 1; }
+static inline int MY_OMP_GET_THREAD_NUM() { return 0; }
+static inline void MY_OMP_SET_NUM_THREADS(int) {}
 #endif
 
 // Prec-switching name macros (respond to SINGLE), used in lib & test sources
@@ -194,23 +206,22 @@
 #include <finufft/fft.h> // (must come after complex.h)
 
 // group together a bunch of type 3 rescaling/centering/phasing parameters:
-#define TYPE3PARAMS FINUFFTIFY(_type3Params)
-typedef struct {
-  FLT X1, C1, D1, h1, gam1; // x dim: X=halfwid C=center D=freqcen h,gam=rescale
-  FLT X2, C2, D2, h2, gam2; // y
-  FLT X3, C3, D3, h3, gam3; // z
-} TYPE3PARAMS;
+template<typename T> struct type3params {
+  T X1, C1, D1, h1, gam1; // x dim: X=halfwid C=center D=freqcen h,gam=rescale
+  T X2, C2, D2, h2, gam2; // y
+  T X3, C3, D3, h3, gam3; // z
+};
 
 typedef struct FINUFFT_PLAN_S { // the main plan object, fully C++
 
   int type;                     // transform type (Rokhlin naming): 1,2 or 3
   int dim;                      // overall dimension: 1,2 or 3
   int ntrans;          // how many transforms to do at once (vector or "many" mode)
   BIGINT nj;           // num of NU pts in type 1,2 (for type 3, num input x pts)
   BIGINT nk;           // number of NU freq pts (type 3 only)
   FLT tol;             // relative user tolerance
   int batchSize;       // # strength vectors to group together for FFTW, etc
   int nbatch;          // how many batches done to cover all ntrans vectors
 
   BIGINT ms;           // number of modes in x (1) dir (historical CMCL name) = N1
   BIGINT mt;           // number of modes in y (2) direction = N2
@@ -222,7 +233,7 @@
   BIGINT nf3;          // " z (3)
   BIGINT nf;           // total # fine grid points (product of the above three)
 
   int fftSign;         // sign in exponential for NUFFT defn, guaranteed to be +-1
 
   FLT *phiHat1;        // FT of kernel in t1,2, on x-axis mode grid
   FLT *phiHat2;        // " y-axis.
@@ -232,7 +243,7 @@
                        // & act on. Usually the largest working array
 
   BIGINT *sortIndices; // precomputed NU pt permutation, speeds spread/interp
   bool didSort;        // whether binsorting used (false: identity perm used)
 
   FLT *X, *Y, *Z;      // for t1,2: ptr to user-supplied NU pts (no new allocs).
                        // for t3: allocated as "primed" (scaled) src pts x'_j, etc
@@ -243,7 +254,7 @@
   CPX *deconv;              // reciprocal of kernel FT, phase, all output NU pts
   CPX *CpBatch;             // working array of prephased strengths
   FLT *Sp, *Tp, *Up;        // internal primed targs (s'_k, etc), allocated
-  TYPE3PARAMS t3P;          // groups together type 3 shift, scale, phase, parameters
+  type3params<FLT> t3P;     // groups together type 3 shift, scale, phase, parameters
   FINUFFT_PLAN innerT2plan; // ptr used for type 2 in step 2 of type 3
 
   // other internal structs; each is C-compatible of course
@@ -255,6 +266,4 @@
 
 } FINUFFT_PLAN_S;
 
-#undef TYPE3PARAMS
-
 #endif // DEFS_H

include/finufft/fft.h

@@ -3,10 +3,10 @@
 
 #ifdef FINUFFT_USE_DUCC0
 #include "ducc0/fft/fftnd_impl.h"
-#define FFTW_FORGET_WISDOM()   // temporary hack since some tests call this unconditionally
-#define FFTW_CLEANUP()         // temporary hack since some tests call this unconditionally
-#define FFTW_CLEANUP_THREADS() // temporary hack since some tests call this
-                               // unconditionally
+// temporary hacks to allow compilation of tests that assume FFTW is used
+static inline void FFTW_FORGET_WISDOM() {}
+static inline void FFTW_CLEANUP() {}
+static inline void FFTW_CLEANUP_THREADS() {}
 #else
 #include "fftw_defs.h"
 #endif

include/finufft/fftw_defs.h

@@ -22,11 +22,7 @@
 // now use this tool (note we replaced typedefs v<=2.0.4, in favor of macros):
 #define FFTW_CPX           FFTWIFY(complex)
 #define FFTW_PLAN          FFTWIFY(plan)
-#define FFTW_ALLOC_RE      FFTWIFY(alloc_real)
 #define FFTW_ALLOC_CPX     FFTWIFY(alloc_complex)
-#define FFTW_PLAN_1D       FFTWIFY(plan_dft_1d)
-#define FFTW_PLAN_2D       FFTWIFY(plan_dft_2d)
-#define FFTW_PLAN_3D       FFTWIFY(plan_dft_3d)
 #define FFTW_PLAN_MANY_DFT FFTWIFY(plan_many_dft)
 #define FFTW_EX            FFTWIFY(execute)
 #define FFTW_DE            FFTWIFY(destroy_plan)

include/finufft/spreadinterp.h

@@ -20,10 +20,12 @@
     NOTE: non-zero values are for experts only, since
     NUMERICAL OUTPUT MAY BE INCORRECT UNLESS finufft_spread_opts.flags=0 !
 */
-#define TF_OMIT_WRITE_TO_GRID        1 // don't add subgrids to out grid (dir=1)
-#define TF_OMIT_EVALUATE_KERNEL      2 // don't evaluate the kernel at all
-#define TF_OMIT_EVALUATE_EXPONENTIAL 4 // omit exp() in kernel (kereval=0 only)
-#define TF_OMIT_SPREADING            8 // don't interp/spread (dir=1: to subgrids)
+enum {
+  TF_OMIT_WRITE_TO_GRID        = 1, // don't add subgrids to out grid (dir=1)
+  TF_OMIT_EVALUATE_KERNEL      = 2, // don't evaluate the kernel at all
+  TF_OMIT_EVALUATE_EXPONENTIAL = 4, // omit exp() in kernel (kereval=0 only)
+  TF_OMIT_SPREADING            = 8  // don't interp/spread (dir=1: to subgrids)
+};
 
 namespace finufft {
 namespace spreadinterp {

include/finufft_eitherprec.h

@@ -24,7 +24,6 @@
 
 // decide which kind of complex numbers FINUFFT_CPX is (four options)
 #ifdef __cplusplus
-#define _USE_MATH_DEFINES
 #include <complex> // C++ type
 #define FINUFFT_COMPLEXIFY(X) std::complex<X>
 #else
@@ -183,4 +182,3 @@ FINUFFT_EXPORT int FINUFFT_CDECL FINUFFTIFY(3d3many)(
 #undef FINUFFT_CPX
 #undef FINUFFT_PLAN
 #undef FINUFFT_PLAN_S
-#undef FINUFFT_TYPE3PARAMS

include/finufft_errors.h

@@ -3,24 +3,26 @@
 
 // ---------- Global error/warning output codes for the library ---------------
 // All documentation is at ../docs/errors.rst (not here):
-#define FINUFFT_WARN_EPS_TOO_SMALL         1
-#define FINUFFT_ERR_MAXNALLOC              2
-#define FINUFFT_ERR_SPREAD_BOX_SMALL       3
-#define FINUFFT_ERR_SPREAD_PTS_OUT_RANGE   4 // DEPRECATED
-#define FINUFFT_ERR_SPREAD_ALLOC           5
-#define FINUFFT_ERR_SPREAD_DIR             6
-#define FINUFFT_ERR_UPSAMPFAC_TOO_SMALL    7
-#define FINUFFT_ERR_HORNER_WRONG_BETA      8
-#define FINUFFT_ERR_NTRANS_NOTVALID        9
-#define FINUFFT_ERR_TYPE_NOTVALID          10
-#define FINUFFT_ERR_ALLOC                  11
-#define FINUFFT_ERR_DIM_NOTVALID           12
-#define FINUFFT_ERR_SPREAD_THREAD_NOTVALID 13
-#define FINUFFT_ERR_NDATA_NOTVALID         14
-#define FINUFFT_ERR_CUDA_FAILURE           15
-#define FINUFFT_ERR_PLAN_NOTVALID          16
-#define FINUFFT_ERR_METHOD_NOTVALID        17
-#define FINUFFT_ERR_BINSIZE_NOTVALID       18
-#define FINUFFT_ERR_INSUFFICIENT_SHMEM     19
-#define FINUFFT_ERR_NUM_NU_PTS_INVALID     20
+enum {
+  FINUFFT_WARN_EPS_TOO_SMALL         = 1,
+  FINUFFT_ERR_MAXNALLOC              = 2,
+  FINUFFT_ERR_SPREAD_BOX_SMALL       = 3,
+  FINUFFT_ERR_SPREAD_PTS_OUT_RANGE   = 4, // DEPRECATED
+  FINUFFT_ERR_SPREAD_ALLOC           = 5,
+  FINUFFT_ERR_SPREAD_DIR             = 6,
+  FINUFFT_ERR_UPSAMPFAC_TOO_SMALL    = 7,
+  FINUFFT_ERR_HORNER_WRONG_BETA      = 8,
+  FINUFFT_ERR_NTRANS_NOTVALID        = 9,
+  FINUFFT_ERR_TYPE_NOTVALID          = 10,
+  FINUFFT_ERR_ALLOC                  = 11,
+  FINUFFT_ERR_DIM_NOTVALID           = 12,
+  FINUFFT_ERR_SPREAD_THREAD_NOTVALID = 13,
+  FINUFFT_ERR_NDATA_NOTVALID         = 14,
+  FINUFFT_ERR_CUDA_FAILURE           = 15,
+  FINUFFT_ERR_PLAN_NOTVALID          = 16,
+  FINUFFT_ERR_METHOD_NOTVALID        = 17,
+  FINUFFT_ERR_BINSIZE_NOTVALID       = 18,
+  FINUFFT_ERR_INSUFFICIENT_SHMEM     = 19,
+  FINUFFT_ERR_NUM_NU_PTS_INVALID     = 20
+};
 #endif

src/finufft.cpp

@@ -91,13 +91,8 @@ namespace common {
 static std::mutex fftw_lock;
 #endif
 
-// We macro because it has no FLT args but gets compiled for both prec's...
-#ifdef SINGLE
-#define SET_NF_TYPE12 set_nf_type12f
-#else
-#define SET_NF_TYPE12 set_nf_type12
-#endif
-int SET_NF_TYPE12(BIGINT ms, finufft_opts opts, finufft_spread_opts spopts, BIGINT *nf)
+static int set_nf_type12(BIGINT ms, finufft_opts opts, finufft_spread_opts spopts,
+                         BIGINT *nf)
 // Type 1 & 2 recipe for how to set 1d size of upsampled array, nf, given opts
 // and requested number of Fourier modes ms. Returns 0 if success, else an
 // error code if nf was unreasonably big (& tell the world).
@@ -684,16 +679,16 @@ int FINUFFT_MAKEPLAN(int type, int dim, BIGINT *n_modes, int iflag, int ntrans,
     }
 
     // determine fine grid sizes, sanity check..
-    int nfier = SET_NF_TYPE12(p->ms, p->opts, p->spopts, &(p->nf1));
+    int nfier = set_nf_type12(p->ms, p->opts, p->spopts, &(p->nf1));
     if (nfier) return nfier; // nf too big; we're done
     p->phiHat1 = (FLT *)malloc(sizeof(FLT) * (p->nf1 / 2 + 1));
     if (dim > 1) {
-      nfier = SET_NF_TYPE12(p->mt, p->opts, p->spopts, &(p->nf2));
+      nfier = set_nf_type12(p->mt, p->opts, p->spopts, &(p->nf2));
       if (nfier) return nfier;
       p->phiHat2 = (FLT *)malloc(sizeof(FLT) * (p->nf2 / 2 + 1));
     }
     if (dim > 2) {
-      nfier = SET_NF_TYPE12(p->mu, p->opts, p->spopts, &(p->nf3));
+      nfier = set_nf_type12(p->mu, p->opts, p->spopts, &(p->nf3));
       if (nfier) return nfier;
       p->phiHat3 = (FLT *)malloc(sizeof(FLT) * (p->nf3 / 2 + 1));
     }

test/CMakeLists.txt

@@ -30,6 +30,7 @@ add_executable(testutils testutils.cpp)
 if(FINUFFT_USE_DUCC0)
   target_compile_definitions(testutils PRIVATE -DFINUFFT_USE_DUCC0)
 endif()
+target_compile_features(testutils PRIVATE cxx_std_17)
 finufft_link_test(testutils)
 add_test(
   NAME run_testutils