From f6ac805ba459601765673b8346e9fd727c7a0f67 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 14 May 2025 09:49:23 -0300
Subject: Pre-compute some exponentials and Γ function
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 log-fourier.cpp | 35 +++++++++++++++++++++--------------
 log-fourier.hpp |  3 +++
 2 files changed, 24 insertions(+), 14 deletions(-)

diff --git a/log-fourier.cpp b/log-fourier.cpp
index a18b61a..e6a84fe 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -4,7 +4,16 @@
 #include <fstream>
 #include <types.hpp>
 
-LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad, Real shift) : N(N), pad(pad), k(k), Δτ(Δτ), shift(shift) {
+Complex Γ(Complex z) {
+  gsl_sf_result logΓ;
+  gsl_sf_result argΓ;
+
+  gsl_sf_lngamma_complex_e((double)z.real(), (double)z.imag(), &logΓ, &argΓ);
+
+  return std::exp((Real)logΓ.val + II * (Real)argΓ.val);
+}
+
+LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad, Real shift) : N(N), pad(pad), k(k), Δτ(Δτ), ts(N), νs(N), Γs(pad * N), shift(shift) {
   τₛ = -0.5 * N;
   ωₛ = -0.5 * N;
   sₛ = -0.5 * pad * N;
@@ -14,6 +23,13 @@ LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Rea
   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");
+  for (unsigned n = 0; n < N; n++) {
+    ts[n] = std::exp(τ(n)) / shift;
+    νs[n] = std::exp(ω(n)) * shift;
+  }
+  for (unsigned n = 0; n < pad * N; n++) {
+    Γs[n] = Γ(k - II * s(n));
+  }
 }
 
 LogarithmicFourierTransform::~LogarithmicFourierTransform() {
@@ -37,20 +53,11 @@ Real LogarithmicFourierTransform::s(unsigned n) const {
 }
 
 Real LogarithmicFourierTransform::t(unsigned n) const {
-  return std::exp(τ(n)) / shift;
+  return ts[n];
 }
 
 Real LogarithmicFourierTransform::ν(unsigned n) const {
-  return std::exp(ω(n)) * shift;
-}
-
-Complex Γ(Complex z) {
-  gsl_sf_result logΓ;
-  gsl_sf_result argΓ;
-
-  gsl_sf_lngamma_complex_e((double)z.real(), (double)z.imag(), &logΓ, &argΓ);
-
-  return std::exp((Real)logΓ.val + II * (Real)argΓ.val);
+  return νs[n];
 }
 
 std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric) {
@@ -72,7 +79,7 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real
     }
     FFTW_EXECUTE(a_to_â);
     for (unsigned n = 0; n < pad*N; n++) {
-      â[(pad*N / 2 + n) % (pad*N)] *= std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n));
+      â[(pad*N / 2 + n) % (pad*N)] *= std::pow(II * σ, II * s(n) - k) * Γs[n];
     }
     FFTW_EXECUTE(â_to_a);
     for (unsigned n = 0; n < N; n++) {
@@ -104,7 +111,7 @@ std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex
     }
     FFTW_EXECUTE(a_to_â);
     for (unsigned n = 0; n < pad*N; n++) {
-      â[(pad*N / 2 + n) % (pad*N)] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n));
+      â[(pad*N / 2 + n) % (pad*N)] *= std::pow(-II * σ, II * s(n) - k) * Γs[n];
     }
     FFTW_EXECUTE(â_to_a);
     for (unsigned n = 0; n < N; n++) {
diff --git a/log-fourier.hpp b/log-fourier.hpp
index ee85775..730639c 100644
--- a/log-fourier.hpp
+++ b/log-fourier.hpp
@@ -20,6 +20,9 @@ private:
   Real τₛ;
   Real ωₛ;
   Real sₛ;
+  std::vector<Real> ts;
+  std::vector<Real> νs;
+  std::vector<Complex> Γs;
 public:
   Real shift;
   LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad = 4, Real shift = 0.5);
-- 
cgit v1.2.3-70-g09d2