From 947daeb85321ed804bc3142623844b2617cb1b3e Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 30 Apr 2025 11:23:38 -0300
Subject: Support long double in the integrator

---
 log-fourier.cpp | 55 ++++++++++++++++++++++++++++---------------------------
 1 file changed, 28 insertions(+), 27 deletions(-)

(limited to 'log-fourier.cpp')

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₂ₙ₊₂
         );
   }
-- 
cgit v1.2.3-70-g09d2