From 238b7fd595ebc25be58905beaaf602838c770fb1 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 14 May 2025 10:23:47 -0300
Subject: Precompute some more exponentials

---
 log-fourier.cpp | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

(limited to 'log-fourier.cpp')

diff --git a/log-fourier.cpp b/log-fourier.cpp
index e6a84fe..cb57229 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -13,7 +13,7 @@ Complex Γ(Complex z) {
   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) {
+LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad, Real shift) : N(N), pad(pad), k(k), Δτ(Δτ), ts(N), νs(N), Γs(pad * N), exp1kω(N), exp1kτ(N), expkω(N), expkτ(N), shift(shift) {
   τₛ = -0.5 * N;
   ωₛ = -0.5 * N;
   sₛ = -0.5 * pad * N;
@@ -26,6 +26,10 @@ LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Rea
   for (unsigned n = 0; n < N; n++) {
     ts[n] = std::exp(τ(n)) / shift;
     νs[n] = std::exp(ω(n)) * shift;
+    exp1kτ[n] = std::exp((1 - k) * τ(n));
+    exp1kω[n] = std::exp((1 - k) * ω(n));
+    expkτ[n] = std::exp(-k * τ(n));
+    expkω[n] = std::exp(-k * ω(n));
   }
   for (unsigned n = 0; n < pad * N; n++) {
     Γs[n] = Γ(k - II * s(n));
@@ -70,9 +74,9 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real
   for (Real σ : σs) {
     for (unsigned n = 0; n < pad*N; n++) {
       if (n < N) {
-        a[n] = c[n] * std::exp((1 - k) * τ(n));
+        a[n] = c[n] * exp1kτ[n];
       } else if (n >= (pad - 1) * N) {
-        a[n] = c[pad*N-n-1] * std::exp((1 - k) * τ(pad*N-n-1));
+        a[n] = c[pad*N-n-1] * exp1kτ[pad*N-n-1];
       } else {
         a[n] = 0;
       }
@@ -83,7 +87,7 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real
     }
     FFTW_EXECUTE(â_to_a);
     for (unsigned n = 0; n < N; n++) {
-      ĉ[n] += std::exp(-k * ω(n)) * a[(pad - 1)*N+n] / (Real)(pad*N);
+      ĉ[n] += expkω[n] * a[(pad - 1)*N+n] / (Real)(pad*N);
     }
   }
 
@@ -102,9 +106,9 @@ std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex
   for (Real σ : σs) {
     for (unsigned n = 0; n < pad * N; n++) {
       if (n < N) {
-        a[n] = (ĉ[n].real() + II * σ * ĉ[n].imag()) * std::exp((1 - k) * ω(n));
+        a[n] = (ĉ[n].real() + II * σ * ĉ[n].imag()) * exp1kω[n];
       } else if (n >= (pad - 1) * N) {
-        a[n] = (ĉ[pad*N-n-1].real() - II * σ * ĉ[pad*N-n-1].imag()) * std::exp((1 - k) * ω(pad*N-n-1));
+        a[n] = (ĉ[pad*N-n-1].real() - II * σ * ĉ[pad*N-n-1].imag()) * exp1kω[pad*N-n-1];
       } else {
         a[n] = 0;
       }
@@ -115,7 +119,7 @@ std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex
     }
     FFTW_EXECUTE(â_to_a);
     for (unsigned n = 0; n < N; n++) {
-      c[n] += std::exp(-k * τ(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI);
+      c[n] += expkτ[n] * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI);
     }
   }
 
-- 
cgit v1.2.3-70-g09d2