Pre-compute some exponentials and Γ function

1 files changed, 21 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++) {