Make N be the smallest power of 2 larger than 3*K+1

4 files changed, 86 insertions, 5 deletions

diff --git a/docs/nstrophy_doc/nstrophy_doc.tex b/docs/nstrophy_doc/nstrophy_doc.tex
index 418520f..35c42ec 100644
--- a/docs/nstrophy_doc/nstrophy_doc.tex
+++ b/docs/nstrophy_doc/nstrophy_doc.tex
@@ -108,9 +108,9 @@ in which $\mathcal F^*$ is defined like $\mathcal F$ but with the opposite phase
 \end{equation}
 provided
 \begin{equation}
-  N_i>4K_i.
+  N_i>3K_i.
 \end{equation}
-Indeed, $\sum_{n_i=0}^{N_i}e^{-\frac{2i\pi}{N_i}n_im_i}$ vanishes unless $m_i=0\%N_i$ (in which $\%N_i$ means `modulo $N_i$'), and, if $p,q\in\mathcal K$, then $|p_i+q_i|\leqslant2K_i$.
+Indeed, $\sum_{n_i=0}^{N_i}e^{-\frac{2i\pi}{N_i}n_im_i}$ vanishes unless $m_i=0\%N_i$ (in which $\%N_i$ means `modulo $N_i$'), and, if $p,q,k\in\mathcal K$, then $|p_i+q_i-k_i|\leqslant3K_i$, so, as long as $N_i>3K_i$, then $(p_i+q_i-k_i)=0\%N_i$ implies $p_i+q_i=k_i$.
 Therefore,
 \begin{equation}
   T(\hat\varphi,k)
diff --git a/src/int_tools.c b/src/int_tools.c
new file mode 100644
index 0000000..13e2cf9
--- /dev/null
+++ b/src/int_tools.c
@@ -0,0 +1,26 @@
+#include <math.h>
+#include "int_tools.h"
+
+// return smallest power of 2 that is > x
+int smallest_pow2(
+  int x
+){
+  return ipow(2,((int)log2(x)+1));
+}
+
+// integer power
+int ipow(
+  int x,
+  int n
+){
+  int out=1;
+  while (n>0)
+  {
+    if (n%2==1){
+      out*=x; 
+    }
+    n/=2;
+    x*=x;
+  }
+  return out;
+}
diff --git a/src/int_tools.h b/src/int_tools.h
new file mode 100644
index 0000000..ab008f6
--- /dev/null
+++ b/src/int_tools.h
@@ -0,0 +1,10 @@
+#ifndef INTTOOLS_H
+#define INTTOOLS_H
+
+// return smallest power of 2 that is larger than x
+int smallest_pow2(int x);
+
+// integer power
+int ipow(int x, int n);
+
+#endif
diff --git a/src/main.c b/src/main.c
index 2fc2fcd..9999de5 100644
--- a/src/main.c
+++ b/src/main.c
@@ -9,6 +9,7 @@
 #include "navier-stokes.h"
 #include "driving.h"
 #include "init.h"
+#include "int_tools.h"
 
 // structure to store parameters, to make it easier to pass parameters to CLI functions
 typedef struct nstrophy_parameters {
@@ -26,6 +27,8 @@ typedef struct nstrophy_parameters {
 
 // usage message
 int print_usage();
+// print parameters
+int print_params(nstrophy_parameters parameters, unsigned int driving, unsigned int init, FILE* file);
 
 // read command line arguments
 int read_args(int argc, const char* argv[], char** params, unsigned int* driving_force, unsigned int* command, unsigned int* init, unsigned int* nthreads);
@@ -82,6 +85,10 @@ int main (
   // set initial condition
   u0=set_init(init, parameters);
 
+  // print parameters
+  print_params(parameters, driving, init, stderr);
+  print_params(parameters, driving, init, stdout);
+
   // run command
   if (command==COMMAND_UK){
     uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads);
@@ -112,6 +119,44 @@ int print_usage(){
   return(0);
 }
 
+// print parameters
+int print_params(
+  nstrophy_parameters parameters,
+  unsigned int driving,
+  unsigned int init,
+  FILE* file
+){
+  fprintf(file,"# K1=%d, K2=%d, N1=%d, N2=%d, nu=%.15e, delta=%.15e, L=%.15e", parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nu, parameters.delta, parameters.L);
+
+  switch(driving){
+  case DRIVING_TEST:
+    fprintf(file,", driving=test");
+    break;
+  case DRIVING_ZERO:
+    fprintf(file,", driving=zero");
+    break;
+  default:
+    fprintf(file,", driving=test");
+    break;
+  }
+
+  switch(init){
+  case INIT_RANDOM:
+    fprintf(file,", init=random");
+    break;
+  case INIT_GAUSSIAN:
+    fprintf(file,", init=gaussian");
+    break;
+  default:
+    fprintf(file,", init=gaussian");
+    break;
+  }
+
+  fprintf(file,"\n");
+
+  return 0;
+}
+
 // read command line arguments
 #define CP_FLAG_PARAMS 1
 #define CP_FLAG_DRIVING 2
@@ -312,12 +357,12 @@ int read_params(
     free(buffer_rhs);
   }
 
-  // if N not set explicitly, set it
+  // if N not set explicitly, set it to the smallest power of 2 that is >3*K+1 (the fft is faster on vectors whose length is a power of 2)
   if (!setN1){
-    parameters->N1=4*(parameters->K1)+1;
+    parameters->N1=smallest_pow2(3*(parameters->K1));
   }
   if (!setN2){
-    parameters->N2=4*(parameters->K2)+1;
+    parameters->N2=smallest_pow2(3*(parameters->K2));
   }
 
   return(0);