Support long double in the integrator

6 files changed, 84 insertions, 36 deletions

diff --git a/.gitignore b/.gitignore
index fe5f9c3..35be3ad 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,6 +3,7 @@
 correlations
 fourier_integrator
 log-fourier_integrator
+log-fourier_integrator_long
 vgcore.*
 walk
 get_energy
diff --git a/Makefile b/Makefile
index 6cf8a52..c5cac94 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log_get_energy
+all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log-fourier_integrator_long log_get_energy
 
 CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1
 
@@ -20,12 +20,21 @@ p-spin.o: p-spin.cpp p-spin.hpp types.hpp
 log-fourier.o: log-fourier.cpp log-fourier.hpp types.hpp
 	$(CC) log-fourier.cpp -c -o log-fourier.o
 
+p-spin_long.o: p-spin.cpp p-spin.hpp types.hpp
+	$(CC) p-spin.cpp -D QUAD_PRECISION -c -o p-spin_long.o
+
+log-fourier_long.o: log-fourier.cpp log-fourier.hpp types.hpp
+	$(CC) log-fourier.cpp -D QUAD_PRECISION -c -o log-fourier_long.o
+
 fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o p-spin.o
 	$(CC) fourier_integrator.cpp fourier.o p-spin.o -lfftw3 -lgsl -o fourier_integrator
 
 log-fourier_integrator: log-fourier_integrator.cpp log-fourier.hpp log-fourier.o p-spin.o
 	$(CC) log-fourier_integrator.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log-fourier_integrator
 
+log-fourier_integrator_long: log-fourier_integrator.cpp log-fourier.hpp log-fourier_long.o p-spin_long.o
+	$(CC) log-fourier_integrator.cpp log-fourier_long.o p-spin_long.o -D QUAD_PRECISION -lfftw3l -lgsl -o log-fourier_integrator_long
+
 get_energy: get_energy.cpp fourier.hpp fourier.o p-spin.o
 	$(CC) get_energy.cpp fourier.o p-spin.o -lfftw3 -lgsl -o get_energy
 
diff --git a/log-fourier.cpp b/log-fourier.cpp
index 47d5c5c..9ae70be 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -2,25 +2,26 @@
 #include "p-spin.hpp"
 #include <complex>
 #include <fstream>
+#include <types.hpp>
 
 LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad) : N(N), pad(pad), k(k), Δτ(Δτ) {
   τₛ = -0.5 * N;
   ωₛ = -0.5 * N;
   sₛ = -0.5 * pad * N;
-  a = reinterpret_cast<Complex*>(fftw_alloc_complex(pad*N));
-  â = reinterpret_cast<Complex*>(fftw_alloc_complex(pad*N));
-  fftw_import_wisdom_from_filename("fftw.wisdom");
-  a_to_â = fftw_plan_dft_1d(pad*N, reinterpret_cast<fftw_complex*>(a), reinterpret_cast<fftw_complex*>(â), FFTW_BACKWARD, 0);
-  â_to_a = fftw_plan_dft_1d(pad*N, reinterpret_cast<fftw_complex*>(â), reinterpret_cast<fftw_complex*>(a), FFTW_BACKWARD, 0);
-  fftw_export_wisdom_to_filename("fftw.wisdom");
+  a = reinterpret_cast<Complex*>(FFTW_ALLOC_COMPLEX(pad*N));
+  â = reinterpret_cast<Complex*>(FFTW_ALLOC_COMPLEX(pad*N));
+  FFTW_IMPORT_WISDOM("fftw.wisdom");
+  a_to_â = FFTW_PLAN_DFT_1D(pad*N, reinterpret_cast<FFTW_COMPLEX*>(a), reinterpret_cast<FFTW_COMPLEX*>(â), FFTW_BACKWARD, 0);
+  â_to_a = FFTW_PLAN_DFT_1D(pad*N, reinterpret_cast<FFTW_COMPLEX*>(â), reinterpret_cast<FFTW_COMPLEX*>(a), FFTW_BACKWARD, 0);
+  FFTW_EXPORT_WISDOM("fftw.wisdom");
 }
 
 LogarithmicFourierTransform::~LogarithmicFourierTransform() {
-  fftw_destroy_plan(a_to_â);
-  fftw_destroy_plan(â_to_a);
-  fftw_free(a);
-  fftw_free(â);
-  fftw_cleanup();
+  FFTW_DESTROY_PLAN(a_to_â);
+  FFTW_DESTROY_PLAN(â_to_a);
+  FFTW_FREE(a);
+  FFTW_FREE(â);
+  FFTW_CLEANUP();
 }
 
 Real LogarithmicFourierTransform::τ(unsigned n) const {
@@ -47,9 +48,9 @@ Complex Γ(Complex z) {
   gsl_sf_result logΓ;
   gsl_sf_result argΓ;
 
-  gsl_sf_lngamma_complex_e(z.real(), z.imag(), &logΓ, &argΓ);
+  gsl_sf_lngamma_complex_e((double)z.real(), (double)z.imag(), &logΓ, &argΓ);
 
-  return exp(logΓ.val + 1i * argΓ.val);
+  return exp((Real)logΓ.val + II * (Real)argΓ.val);
 }
 
 std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric) {
@@ -67,13 +68,13 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real
         a[n] = 0;
       }
     }
-    fftw_execute(a_to_â);
+    FFTW_EXECUTE(a_to_â);
     for (unsigned n = 0; n < pad*N; n++) {
-      â[(pad*N / 2 + n) % (pad*N)] *= pow(1i * σ, 1i * s(n) - k) * Γ(k - 1i * s(n));
+      â[(pad*N / 2 + n) % (pad*N)] *= std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n));
     }
-    fftw_execute(â_to_a);
+    FFTW_EXECUTE(â_to_a);
     for (unsigned n = 0; n < N; n++) {
-      ĉ[n] += exp(-k * ω(n)) * a[(pad - 1)*N+n] / (Real)(pad*N);
+      ĉ[n] += std::exp(-k * ω(n)) * a[(pad - 1)*N+n] / (Real)(pad*N);
     }
   }
 
@@ -87,21 +88,21 @@ std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex
     for (unsigned n = 0; n < pad * N; n++) {
       if (n < N) {
         if (σ < 0) {
-          a[n] = std::conj(ĉ[n]) * exp((1 - k) * ω(n));
+          a[n] = std::conj(ĉ[n]) * std::exp((1 - k) * ω(n));
         } else {
-          a[n] = ĉ[n] * exp((1 - k) * ω(n));
+          a[n] = ĉ[n] * std::exp((1 - k) * ω(n));
         }
       } else {
         a[n] = 0;
       }
     }
-    fftw_execute(a_to_â);
+    FFTW_EXECUTE(a_to_â);
     for (unsigned n = 0; n < pad*N; n++) {
-      â[(pad*N / 2 + n) % (pad*N)] *= pow(-1i * σ, 1i * s(n) - k) * Γ(k - 1i * s(n));
+      â[(pad*N / 2 + n) % (pad*N)] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n));
     }
-    fftw_execute(â_to_a);
+    FFTW_EXECUTE(â_to_a);
     for (unsigned n = 0; n < N; n++) {
-      c[n] += exp(-k * τ(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI);
+      c[n] += std::exp(-k * τ(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI);
     }
   }
 
@@ -113,7 +114,7 @@ std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ,
 }
 
 void logFourierSave(const std::vector<Real>& C, const std::vector<Real>& R, const std::vector<Complex>& Ct, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k) {
-    unsigned N = pow(2, log2n);
+    unsigned N = std::pow(2, log2n);
     std::ofstream outfile(logFourierFile("C", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::out | std::ios::binary);
     outfile.write((const char*)(C.data()), N * sizeof(Real));
     outfile.close();
@@ -141,7 +142,7 @@ bool logFourierLoad(std::vector<Real>& C, std::vector<Real>& R, std::vector<Comp
     return false;
   }
 
-  unsigned N = pow(2, log2n);
+  unsigned N = std::pow(2, log2n);
 
   cfile.read((char*)(C.data()), N * sizeof(Real));
   cfile.close();
@@ -175,7 +176,7 @@ Real estimateZ(LogarithmicFourierTransform& fft, const std::vector<Real>& C, con
   auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
   Real Γ₀ = 1.0;
 
-  return ((2 * Γ₀ * std::conj(Rt[0]) + pow(β, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
+  return ((2 * Γ₀ * std::conj(Rt[0]) + std::pow(β, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
 }
 
 Real energy(const LogarithmicFourierTransform& fft, std::vector<Real>&  C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real β) {
@@ -188,7 +189,7 @@ Real energy(const LogarithmicFourierTransform& fft, std::vector<Real>&  C, const
     Real f₂ₙ₊₂ = R[2*n+2] * df(λ, p, s, C[2*n+2]);
     E += (h₂ₙ + h₂ₙ₊₁) / 6 * (
           (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ
-          + pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
+          + std::pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
           + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂
         );
   }
diff --git a/log-fourier.hpp b/log-fourier.hpp
index e232bd6..97c0e15 100644
--- a/log-fourier.hpp
+++ b/log-fourier.hpp
@@ -11,8 +11,8 @@ class LogarithmicFourierTransform {
 private:
   Complex* a;
   Complex* â;
-  fftw_plan a_to_â;
-  fftw_plan â_to_a;
+  FFTW_PLAN a_to_â;
+  FFTW_PLAN â_to_a;
   unsigned N;
   unsigned pad;
   Real k;
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 9bbf554..7ac76b5 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -79,7 +79,7 @@ int main(int argc, char* argv[]) {
 
   LogarithmicFourierTransform fft(N, k, Δτ, 4);
 
-  Real Γ₀ = 1.0;
+  Real Γ₀ = 1;
   Real μₜ₋₁ = Γ₀;
   if (τ₀ > 0) {
     μₜ₋₁ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
@@ -105,7 +105,7 @@ int main(int argc, char* argv[]) {
       Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n));
 
       Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μₜ₋₁, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
-      Ȓₜ₋₁[n] = 1.0 / (μₜ₋₁ + 1i * fft.ν(n));
+      Ȓₜ₋₁[n] = (Real)1.0 / (μₜ₋₁ + II * fft.ν(n));
     }
   } else {
     logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, k);
@@ -130,8 +130,8 @@ int main(int argc, char* argv[]) {
       std::vector<Complex> Ĉₜ₊₁(N);
       std::vector<Complex> Ȓₜ₊₁(N);
       for (unsigned n = 0; n < N; n++) {
-        Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + 1i * fft.ν(n));
-        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + 1i * fft.ν(n));
+        Ȓₜ₊₁[n] = ((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + II * fft.ν(n));
+        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n));
       }
       std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
       std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
@@ -153,7 +153,7 @@ int main(int argc, char* argv[]) {
       }
 
       if (stepsUp > stepsToRespond) {
-        γ = std::max(γ/2, 1e-4);
+        γ = std::max(γ/2, (Real)1e-4);
         stepsUp = 0;
         ΔCmin = ΔCₜ;
       }
diff --git a/types.hpp b/types.hpp
index 05eceb4..71556fd 100644
--- a/types.hpp
+++ b/types.hpp
@@ -1,6 +1,43 @@
 #pragma once
 #include <complex>
+#include <fftw3.h>
+
+using namespace std::complex_literals;
+
+#ifndef QUAD_PRECISION
+
+#define FFTW_PLAN fftw_plan
+#define FFTW_ALLOC_COMPLEX fftw_alloc_complex
+#define FFTW_IMPORT_WISDOM fftw_import_wisdom_from_filename
+#define FFTW_EXPORT_WISDOM fftw_export_wisdom_to_filename
+#define FFTW_PLAN_DFT_1D fftw_plan_dft_1d
+#define FFTW_DESTROY_PLAN fftw_destroy_plan
+#define FFTW_FREE fftw_free
+#define FFTW_CLEANUP fftw_cleanup
+#define FFTW_EXECUTE fftw_execute
+#define FFTW_COMPLEX fftw_complex
 
 using Real = double;
+
+#define II 1i
+
+#else
+
+#define FFTW_PLAN fftwl_plan
+#define FFTW_ALLOC_COMPLEX fftwl_alloc_complex
+#define FFTW_IMPORT_WISDOM fftwl_import_wisdom_from_filename
+#define FFTW_EXPORT_WISDOM fftwl_export_wisdom_to_filename
+#define FFTW_PLAN_DFT_1D fftwl_plan_dft_1d
+#define FFTW_DESTROY_PLAN fftwl_destroy_plan
+#define FFTW_FREE fftwl_free
+#define FFTW_CLEANUP fftwl_cleanup
+#define FFTW_EXECUTE fftwl_execute
+#define FFTW_COMPLEX fftwl_complex
+
+using Real = long double;
+
+#define II 1il
+
+#endif
+
 using Complex = std::complex<Real>;
-using namespace std::complex_literals;