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 +++++++++++++++++++++--------------
1 file changed, 21 insertions(+), 14 deletions(-)
(limited to 'log-fourier.cpp')
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++) {
--
cgit v1.2.3-70-g09d2