summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitignore10
-rw-r--r--Makefile38
-rw-r--r--fourier.cpp86
-rw-r--r--fourier.hpp32
-rw-r--r--fourier_integrator.cpp172
-rw-r--r--get_energy.cpp87
-rw-r--r--integrator.cpp198
-rw-r--r--log-fourier.cpp278
-rw-r--r--log-fourier.hpp58
-rw-r--r--log-fourier_integrator.cpp238
-rw-r--r--log_get_energy.cpp100
-rw-r--r--p-spin.cpp25
-rw-r--r--p-spin.hpp7
-rw-r--r--print_correlations.cpp86
-rw-r--r--types.hpp43
15 files changed, 873 insertions, 585 deletions
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..56f8120
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,10 @@
+*.dat
+*.o
+vgcore.*
+fftw.wisdom
+walk
+correlations
+log-fourier_integrator
+log-fourier_integrator_long
+log_get_energy
+log_get_energy_long
diff --git a/Makefile b/Makefile
index d4c15bf..2bcdee1 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
-all: walk correlations integrator fourier_integrator get_energy
+all: walk correlations log-fourier_integrator log-fourier_integrator_long log_get_energy log_get_energy_long print_correlations print_correlations_long
-CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp
+CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1
walk: walk.cpp
$(CC) walk.cpp -o walk
@@ -8,14 +8,32 @@ walk: walk.cpp
correlations: correlations.cpp
$(CC) correlations.cpp -o correlations
-integrator: integrator.cpp
- $(CC) integrator.cpp -o integrator
+p-spin.o: p-spin.cpp p-spin.hpp types.hpp
+ $(CC) p-spin.cpp -c -o p-spin.o
-fourier.o: fourier.cpp fourier.hpp
- $(CC) fourier.cpp -c -o fourier.o
+p-spin_long.o: p-spin.cpp p-spin.hpp types.hpp
+ $(CC) p-spin.cpp -D QUAD_PRECISION -c -o p-spin_long.o
-fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o
- $(CC) fourier_integrator.cpp fourier.o -lfftw3 -o fourier_integrator
+log-fourier.o: log-fourier.cpp log-fourier.hpp types.hpp
+ $(CC) log-fourier.cpp -c -o log-fourier.o
-get_energy: get_energy.cpp fourier.hpp fourier.o
- $(CC) get_energy.cpp fourier.o -lfftw3 -o get_energy
+log-fourier_long.o: log-fourier.cpp log-fourier.hpp types.hpp
+ $(CC) log-fourier.cpp -D QUAD_PRECISION -c -o log-fourier_long.o
+
+log-fourier_integrator: log-fourier_integrator.cpp log-fourier.hpp log-fourier.o p-spin.o
+ $(CC) log-fourier_integrator.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log-fourier_integrator
+
+log-fourier_integrator_long: log-fourier_integrator.cpp log-fourier.hpp log-fourier_long.o p-spin_long.o
+ $(CC) log-fourier_integrator.cpp log-fourier_long.o p-spin_long.o -D QUAD_PRECISION -lfftw3l -lgsl -o log-fourier_integrator_long
+
+log_get_energy: log_get_energy.cpp log-fourier.hpp log-fourier.o p-spin.o
+ $(CC) log_get_energy.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log_get_energy
+
+log_get_energy_long: log_get_energy.cpp log-fourier.hpp log-fourier_long.o p-spin_long.o
+ $(CC) log_get_energy.cpp log-fourier_long.o p-spin_long.o -D QUAD_PRECISION -lfftw3l -lgsl -o log_get_energy_long
+
+print_correlations: print_correlations.cpp log-fourier.hpp log-fourier.o p-spin.o
+ $(CC) print_correlations.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o print_correlations
+
+print_correlations_long: print_correlations.cpp log-fourier.hpp log-fourier_long.o p-spin_long.o
+ $(CC) print_correlations.cpp log-fourier_long.o p-spin_long.o -D QUAD_PRECISION -lfftw3l -lgsl -o print_correlations_long
diff --git a/fourier.cpp b/fourier.cpp
deleted file mode 100644
index 7808989..0000000
--- a/fourier.cpp
+++ /dev/null
@@ -1,86 +0,0 @@
-#include "fourier.hpp"
-
-inline Real fP(unsigned p, Real q) {
- return 0.5 * pow(q, p);
-}
-
-inline Real dfP(unsigned p, Real q) {
- return 0.5 * p * pow(q, p - 1);
-}
-
-inline Real ddfP(unsigned p, Real q) {
- return 0.5 * p * (p - 1) * pow(q, p - 2);
-}
-
-Real f(Real λ, unsigned p, unsigned s, Real q) {
- return (1 - λ) * fP(p, q) + λ * fP(s, q);
-}
-
-Real df(Real λ, unsigned p, unsigned s, Real q) {
- return (1 - λ) * dfP(p, q) + λ * dfP(s, q);
-}
-
-Real ddf(Real λ, unsigned p, unsigned s, Real q) {
- return (1 - λ) * ddfP(p, q) + λ * ddfP(s, q);
-}
-
-FourierTransform::FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags) : a(2 * n), â(n + 1), Δω(Δω), Δτ(Δτ) {
- plan_r2c = fftw_plan_dft_r2c_1d(2 * n, a.data(), reinterpret_cast<fftw_complex*>(â.data()), flags);
- plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â.data()), a.data(), flags);
-}
-
-FourierTransform::~FourierTransform() {
- fftw_destroy_plan(plan_r2c);
- fftw_destroy_plan(plan_c2r);
- fftw_cleanup();
-}
-
-std::vector<Complex> FourierTransform::fourier(const std::vector<Real>& c) {
- a = c;
- fftw_execute(plan_r2c);
- std::vector<Complex> ĉ(â.size());
- for (unsigned i = 0; i < â.size(); i++) {
- ĉ[i] = â[i] * (Δτ * M_PI);
- }
- return ĉ;
-}
-
-std::vector<Real> FourierTransform::inverse(const std::vector<Complex>& ĉ) {
- â = ĉ;
- fftw_execute(plan_c2r);
- std::vector<Real> c(a.size());
- for (unsigned i = 0; i < a.size(); i++) {
- c[i] = a[i] * (Δω / (2 * M_PI));
- }
- return c;
-}
-
-std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ) {
- return prefix + "_" + std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(y) + "_" + std::to_string(log2n) + "_" + std::to_string(τₘₐₓ) + ".dat";
-}
-
-Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ) {
- Real e = 0;
-
- for (unsigned i = 0; i < C.size() / 2; i++) {
- e += y * R[i] * df(λ, p, s, C[i]) * M_PI * Δτ;
- }
-
- return e;
-}
-
-std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ) {
- std::vector<Real> RddfC(C.size());
- for (unsigned i = 0; i < C.size(); i++) {
- RddfC[i] = R[i] * ddf(λ, p, s, C[i]);
- }
- std::vector<Complex> RddfCt = fft.fourier(RddfC);
-
- std::vector<Real> dfC(C.size());
- for (unsigned i = 0; i < C.size(); i++) {
- dfC[i] = df(λ, p, s, C[i]);
- }
- std::vector<Complex> dfCt = fft.fourier(dfC);
-
- return {RddfCt, dfCt};
-}
diff --git a/fourier.hpp b/fourier.hpp
deleted file mode 100644
index 86335f7..0000000
--- a/fourier.hpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#pragma once
-#include <cmath>
-#include <fftw3.h>
-#include <complex>
-#include <vector>
-
-using Real = double;
-using Complex = std::complex<Real>;
-using namespace std::complex_literals;
-
-Real f(Real λ, unsigned p, unsigned s, Real q);
-Real df(Real λ, unsigned p, unsigned s, Real q);
-Real ddf(Real λ, unsigned p, unsigned s, Real q);
-
-class FourierTransform {
-private:
- std::vector<Real> a;
- std::vector<Complex> â;
- fftw_plan plan_r2c;
- fftw_plan plan_c2r;
- Real Δω;
- Real Δτ;
-public:
- FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags = 0);
- ~FourierTransform();
- std::vector<Complex> fourier(const std::vector<Real>& c);
- std::vector<Real> inverse(const std::vector<Complex>& ĉ);
-};
-
-std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ);
-Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ);
-std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ);
diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
deleted file mode 100644
index 86c3b24..0000000
--- a/fourier_integrator.cpp
+++ /dev/null
@@ -1,172 +0,0 @@
-#include "fourier.hpp"
-#include <getopt.h>
-#include <iostream>
-#include <fstream>
-
-int main(int argc, char* argv[]) {
- unsigned p = 2;
- unsigned s = 2;
- Real λ = 0.5;
- Real τ₀ = 0;
- Real y₀ = 0;
- Real yₘₐₓ = 0.5;
- Real Δy = 0.05;
-
- unsigned log2n = 8;
- Real τₘₐₓ = 20;
-
- unsigned maxIterations = 1000;
- Real ε = 1e-14;
- Real γ = 0;
-
- bool loadData = false;
-
- int opt;
-
- while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:I:g:l")) != -1) {
- switch (opt) {
- case 'p':
- p = atoi(optarg);
- break;
- case 's':
- s = atoi(optarg);
- break;
- case '2':
- log2n = atoi(optarg);
- break;
- case 'T':
- τₘₐₓ = atof(optarg);
- break;
- case 't':
- τ₀ = atof(optarg);
- break;
- case '0':
- y₀ = atof(optarg);
- break;
- case 'y':
- yₘₐₓ = atof(optarg);
- break;
- case 'd':
- Δy = atof(optarg);
- break;
- case 'I':
- maxIterations = (unsigned)atof(optarg);
- break;
- case 'g':
- γ = atof(optarg);
- break;
- case 'l':
- loadData = true;
- break;
- default:
- exit(1);
- }
- }
-
- unsigned n = pow(2, log2n);
-
- Real Δτ = τₘₐₓ / M_PI / n;
- Real Δω = M_PI / τₘₐₓ;
-
- Real z = (-1+sqrt(1+2*τ₀)) / (2 * τ₀);
- Real Γ₀ = 1;
-
- std::vector<Real> C(2 * n);
- std::vector<Real> R(2 * n);
-
- FourierTransform fft(n, Δω, Δτ);
- std::vector<Complex> Ct;
- std::vector<Complex> Rt;
-
- Real y = y₀;
-
- if (!loadData) {
- // start from the exact solution for τ₀ = 0
- for (unsigned i = 0; i < n; i++) {
- Real τ = i * Δτ * M_PI;
- C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
- if (i > 0) {
- C[2 * n - i] = C[i];
- }
- R[i] = exp(-z * τ);
- }
- Ct = fft.fourier(C);
- Rt = fft.fourier(R);
- } else {
- std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
- cfile.read((char*)(C.data()), C.size() * sizeof(Real));
- cfile.close();
- std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
- rfile.read((char*)(R.data()), R.size() * sizeof(Real));
- rfile.close();
-
- Ct = fft.fourier(C);
- Rt = fft.fourier(R);
-
- auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-
- z = ((Γ₀ * std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
- }
-
- while (y += Δy, y <= yₘₐₓ) {
- Real ΔC = 1;;
- unsigned it = 0;
- while (sqrt(ΔC / C.size()) > ε) {
- it++;
- auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-
- for (unsigned i = 0; i < Rt.size(); i++) {
- Real ω = i * Δω;
- Rt[i] = (1.0 + pow(y, 2) * RddfCt[i] * Rt[i]) / (z + 1i * ω);
- }
-
- for (unsigned i = 0; i < Ct.size(); i++) {
- Real ω = i * Δω;
- Ct[i] = (Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(y, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (z + 1i * ω);
- }
-
- std::vector<Real> Cnew = fft.inverse(Ct);
- std::vector<Real> Rnew = fft.inverse(Rt);
- for (unsigned i = n; i < 2 * n; i++) {
- Rnew[i] = 0;
- }
-
- ΔC = 0;
- for (unsigned i = 0; i < Cnew.size(); i++) {
- ΔC += pow(Cnew[i] - C[i], 2);
- ΔC += pow(Rnew[i] - R[i], 2);
- }
-
- for (unsigned i = 0; i < Cnew.size(); i++) {
- C[i] += γ * (Cnew[i] - C[i]);
- }
-
- for (unsigned i = 0; i < Rnew.size(); i++) {
- R[i] += γ * (Rnew[i] - R[i]);
- }
-
- z *= Cnew[0];
-
- std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << sqrt(ΔC / C.size()) << std::endl;
-
- if (it > maxIterations) {
- it = 0;
- γ /= 2;
- }
- }
-
- Real e = energy(C, R, p, s, λ, y, Δτ);
-
- std::cerr << "y " << y << " " << e << " " << z << std::endl;
-
- std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
- outfile.write((const char*)(C.data()), C.size() * sizeof(Real));
- outfile.close();
-
- std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
- outfileR.write((const char*)(R.data()), R.size() * sizeof(Real));
- outfileR.close();
- }
-
- return 0;
-}
diff --git a/get_energy.cpp b/get_energy.cpp
deleted file mode 100644
index a388e0c..0000000
--- a/get_energy.cpp
+++ /dev/null
@@ -1,87 +0,0 @@
-#include "fourier.hpp"
-#include <getopt.h>
-#include <iostream>
-#include <fstream>
-
-int main(int argc, char* argv[]) {
- unsigned p = 2;
- unsigned s = 2;
- Real λ = 0.5;
- Real τ₀ = 0;
- Real y₀ = 0;
- Real yₘₐₓ = 0.5;
- Real Δy = 0.05;
-
- unsigned log2n = 8;
- Real τₘₐₓ = 20;
-
- int opt;
-
- while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:")) != -1) {
- switch (opt) {
- case 'p':
- p = atoi(optarg);
- break;
- case 's':
- s = atoi(optarg);
- break;
- case '2':
- log2n = atoi(optarg);
- break;
- case 'T':
- τₘₐₓ = atof(optarg);
- break;
- case 't':
- τ₀ = atof(optarg);
- break;
- case '0':
- y₀ = atof(optarg);
- break;
- case 'y':
- yₘₐₓ = atof(optarg);
- break;
- case 'd':
- Δy = atof(optarg);
- break;
- default:
- exit(1);
- }
- }
-
- unsigned n = pow(2, log2n);
-
- Real Δτ = τₘₐₓ / M_PI / n;
- Real Δω = M_PI / τₘₐₓ;
-
- Real y = y₀;
-
- FourierTransform fft(n, Δω, Δτ, FFTW_ESTIMATE);
-
- while (y += Δy, y <= yₘₐₓ) {
- std::vector<Real> C(2 * n);
- std::vector<Real> R(2 * n);
-
- std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
- if (cfile.is_open()) {
- cfile.read((char*)(C.data()), C.size() * sizeof(Real));
- cfile.close();
-
- std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
- rfile.read((char*)(R.data()), R.size() * sizeof(Real));
- rfile.close();
-
- Real e = energy(C, R, p, s, λ, y, Δτ);
-
- std::vector<Complex> Ct = fft.fourier(C);
- std::vector<Complex> Rt = fft.fourier(R);
-
- auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-
- Real z = ((std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
-
- std::cout << y << " " << e << " " << Ct[0].real() << " " << z << std::endl;
- }
- }
-
- return 0;
-}
diff --git a/integrator.cpp b/integrator.cpp
deleted file mode 100644
index 8b13a9b..0000000
--- a/integrator.cpp
+++ /dev/null
@@ -1,198 +0,0 @@
-#include <getopt.h>
-#include <vector>
-#include <cmath>
-#include <iostream>
-
-using Real = double;
-
-unsigned p = 2;
-
-Real f(Real q) {
- return 0.5 * pow(q, p);
-}
-
-Real df(Real q) {
- return 0.5 * p * pow(q, p - 1);
-}
-
-Real ddf(Real q) {
- return 0.5 * p * (p - 1) * pow(q, p - 2);
-}
-
-Real integrate(const std::vector<Real>& C, signed τ = std::numeric_limits<unsigned>::max()) {
- Real I = 0;
- if (τ > C.size() - 1) {
- τ = C.size() - 1;
- }
-#pragma omp parallel for reduction(+:I)
- for (unsigned σ = 0; σ < τ; σ++) {
- unsigned τ_σ = τ - σ;
- Real Cτ_σ = (C[τ_σ] + C[τ_σ - 1]) / 2;
- Real dCσ = C[σ + 1] - C[σ];
-
- I += df(Cτ_σ) * dCσ;
- }
- return I;
-}
-
-Real integratePast(const std::vector<Real>& C, signed τ) {
- Real I = 0;
-#pragma omp parallel for reduction(+:I)
- for (signed σ = -C.size() + τ + 3; σ < τ - 2; σ++) {
- signed τ_σ = τ - σ;
-
- Real Cτ_σ = (C[abs(τ_σ)] + C[abs(τ_σ) - 1]) / 2;
- Real Cσ = (C[abs(σ) + 1] + C[abs(σ)]) / 2;
- Real dddC;
- if (τ_σ != 0) {
- dddC = (τ_σ / abs(τ_σ)) * (C[abs(τ_σ)+2] - 2 * C[abs(τ_σ)+1] + 2 * C[abs(τ_σ)-1] - C[abs(τ_σ)-2]) / 2;
- } else {
- dddC = 0;
- }
-
- I += dddC * ddf(Cτ_σ) * Cσ;
- }
-#pragma omp parallel for reduction(+:I)
- for (signed σ = -C.size() + τ + 3; σ < -1; σ++) {
- signed τ_σ = τ - σ;
-
- Real Cτ_σ = (C[abs(τ_σ)] + C[abs(τ_σ) - 1]) / 2;
- Real dddC;
- if (σ != 0) {
- dddC = -(σ / abs(σ)) * (C[abs(σ)+2] - 2 * C[abs(σ)+1] + 2 * C[abs(σ)-1] - C[abs(σ)-2]) / 2;
- } else {
- dddC = 0;
- }
-
- I += dddC * df(Cτ_σ);
- }
- return I;
-}
-
-Real integrateDelay(const std::vector<Real>& C, unsigned τ, Real Δτ, Real τ₀) {
- Real I = 0;
-#pragma omp parallel for reduction(+:I)
- for (signed σ = 2; σ < C.size() - τ - 2; σ++) {
- unsigned τ_σ = τ + σ;
- Real dC = -(C[σ+1] - C[σ-1]) / 2;
- Real dddC = -(C[σ+2] - 2 * C[σ+1] + 2 * C[σ-1] - C[σ-2]) / 2;
-
- I += (dC - pow(τ₀ / Δτ, 2) * dddC) * exp(-(τ_σ * Δτ / τ₀));
- }
- return I / τ₀;
-}
-
-Real energy(const std::vector<Real>& C, Real Δτ, Real τ₀) {
- Real I = 0;
- for (unsigned σ = 0; σ < C.size() - 1; σ++) {
- Real Cσ = (C[σ] + C[σ + 1]) / 2;
- Real dC = (C[σ + 1] - C[σ]) / Δτ;
-
- Real dddC = 0;
- if (σ > 1 && σ < C.size() - 2 && C.size() > 3) {
- dddC = (C[σ+1] - 3 * C[σ] + 3 * C[σ-1] - C[σ-2]) / pow(Δτ, 3);
- }
- I += Δτ * df(Cσ) * (dC - pow(τ₀, 2) * dddC);
- }
- return I;
-}
-
-int main(int argc, char* argv[]) {
- Real Δτ = 1e-3;
- Real τₘₐₓ = 1e3;
- Real τ₀ = 0;
- Real y = 0.5;
- unsigned iterations = 10;
-
- int opt;
-
- while ((opt = getopt(argc, argv, "d:T:t:y:I:")) != -1) {
- switch (opt) {
- case 'd':
- Δτ = atof(optarg);
- break;
- case 'T':
- τₘₐₓ = atof(optarg);
- break;
- case 't':
- τ₀ = atof(optarg);
- break;
- case 'y':
- y = atof(optarg);
- break;
- case 'I':
- iterations = (unsigned)atof(optarg);
- break;
- default:
- exit(1);
- }
- }
-
- Real z = 0.4794707565634420155347;
- Real Γ₀ = 1;
-
- Real τ = 0;
- std::vector<Real> C;
- C.reserve(τₘₐₓ / Δτ + 1);
-
- C.push_back(1);
-
-// while (std::cout << τ << " " << C.back() << std::endl, τ < τₘₐₓ) {
- while (τ < τₘₐₓ) {
- τ += Δτ;
- Real dC = -(z - 2 * pow(y, 2)) * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C);
- C.push_back(C.back() + Δτ * dC);
- }
-
-
- for (unsigned it = 0; it < iterations; it++) {
-
- τ = 0;
- std::vector<Real> C2;
- C2.reserve(τₘₐₓ / Δτ + 1);
- C2.push_back(1);
- while (τ < τₘₐₓ) {
- τ += Δτ;
- Real dC = -(z - 2 * pow(y, 2)) * C2.back() + integrateDelay(C, C2.size() - 1, Δτ, τ₀) - 2 / Γ₀ * pow(y, 2) * (integrate(C2) - pow(τ₀ / Δτ, 2) * integratePast(C, C2.size()-1));
- C2.push_back(C2.back() + Δτ * dC);
- }
-
- Real error = 0;
-
- for (unsigned i = 0; i < std::min(C.size(), C2.size()); i++) {
- error += pow(C[i] - C2[i], 2);
- }
-
- std::cerr << "Iteration " << it << ": " << sqrt(error / C.size()) << " " << z << std::endl;
-
- C = C2;
- }
- /*
- Real zNew = (2.0 * ((C[2] - 2 * C[1] + C[0]) / pow(Δτ, 2) - pow(τ₀, 2) * (C[4] - 4 * C[3] + 6 * C[2] - 4 * C[1] + C[0]) / pow(Δτ, 4)));
- Real zNew = (2.0 * ((C[2] - 2 * C[1] + C[0]) / pow(Δτ, 2) - pow(τ₀, 2) * (C[4] - 4 * C[3] + 6 * C[2] - 4 * C[1] + C[0]) / pow(Δτ, 4)));
-// Real zNew = (2.0 * ((83 * C[6] - 245 * C[5] + 101 * C[4] + 254 * C[3] - 31 * C[2] - 377 * C[1] + 215 * C[0]) / (132 * pow(Δτ, 2)) - pow(τ₀, 2) * (3 * C[6] - 7 * C[5] + C[4] + 6 * C[3] + C[2] - 7 * C[1] + 3 * C[0]) / (11 * pow(Δτ, 4))));
- z = z / zNew;
- τ = 0;
- C.clear();
- C.reserve(τₘₐₓ / Δτ + 1);
-
- C.push_back(1);
-
-// while (std::cout << τ << " " << C.back() << std::endl, τ < τₘₐₓ) {
- while (τ < τₘₐₓ) {
- τ += Δτ;
- Real dC = -z * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C);
- C.push_back(C.back() + Δτ * dC);
- }
- */
-
- τ = 0;
- for (Real Ci : C) {
- std::cout << τ << " " << Ci << std::endl;
- τ += Δτ;
- }
-
- std::cerr << - 2 * y / Γ₀ * energy(C, Δτ, τ₀) << std::endl;
-
- return 0;
-}
diff --git a/log-fourier.cpp b/log-fourier.cpp
new file mode 100644
index 0000000..16a7f3f
--- /dev/null
+++ b/log-fourier.cpp
@@ -0,0 +1,278 @@
+#include "log-fourier.hpp"
+#include "p-spin.hpp"
+
+#include <fstream>
+
+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), exp1kω(N), exp1kτ(N), expkω(N), expkτ(N), shift(shift) {
+ τₛ = -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("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");
+ 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));
+ }
+}
+
+LogarithmicFourierTransform::~LogarithmicFourierTransform() {
+ FFTW_DESTROY_PLAN(a_to_â);
+ FFTW_DESTROY_PLAN(â_to_a);
+ FFTW_FREE(a);
+ FFTW_FREE(â);
+ FFTW_CLEANUP();
+}
+
+Real LogarithmicFourierTransform::τ(unsigned n) const {
+ return Δτ * (n + τₛ);
+}
+
+Real LogarithmicFourierTransform::ω(unsigned n) const {
+ return Δτ * (n + ωₛ);
+}
+
+Real LogarithmicFourierTransform::s(unsigned n) const {
+ return (n + sₛ) * 2*M_PI / (pad * N * Δτ);
+}
+
+Real LogarithmicFourierTransform::t(unsigned n) const {
+ return ts[n];
+}
+
+Real LogarithmicFourierTransform::ν(unsigned n) const {
+ return νs[n];
+}
+
+std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric) {
+ std::vector<Complex> ĉ(N);
+ std::vector<Real> σs = {1};
+ /* c is either even or zero for negative arguments */
+ if (symmetric){
+ σs.push_back(-1);
+ }
+ for (Real σ : σs) {
+ for (unsigned n = 0; n < pad*N; n++) {
+ if (n < N) {
+ a[n] = c[n] * exp1kτ[n];
+ } else if (n >= (pad - 1) * N) {
+ a[n] = c[pad*N-n-1] * exp1kτ[pad*N-n-1];
+ } else {
+ a[n] = 0;
+ }
+ }
+ 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) * Γs[n];
+ }
+ FFTW_EXECUTE(â_to_a);
+ for (unsigned n = 0; n < N; n++) {
+ ĉ[n] += expkω[n] * a[(pad - 1)*N+n] / (Real)(pad*N);
+ }
+ }
+
+ for (unsigned n = 0; n < N; n++) {
+ ĉ[n] -= ĉ[N - 1];
+ if (symmetric) ĉ[n] = ĉ[n].real();
+ ĉ[n] /= shift;
+ }
+
+ return ĉ;
+}
+
+std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex>& ĉ) {
+ std::vector<Real> c(N);
+ std::vector<Real> σs = {1, -1};
+ for (Real σ : σs) {
+ for (unsigned n = 0; n < pad * N; n++) {
+ if (n < 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()) * exp1kω[pad*N-n-1];
+ } else {
+ a[n] = 0;
+ }
+ }
+ 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) * Γs[n];
+ }
+ FFTW_EXECUTE(â_to_a);
+ for (unsigned n = 0; n < N; n++) {
+ c[n] += expkτ[n] * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI);
+ }
+ }
+
+ for (unsigned n = 0; n < N; n++) {
+ c[n] -= c[N - 1];
+ c[n] *= shift;
+ }
+
+ return c;
+}
+
+std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift) {
+ return prefix + "_" + std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(β) + "_" + std::to_string(log2n) + "_" + std::to_string(Δτ) + "_" + std::to_string(shift) + ".dat";
+}
+
+void logFourierSave(const std::vector<Real>& C, const std::vector<Real>& R, const std::vector<Complex>& Ĉ, const std::vector<Complex>& Ȓ, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k) {
+ 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();
+
+ std::ofstream outfileCt(logFourierFile("Ct", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::out | std::ios::binary);
+ outfileCt.write((const char*)(Ĉ.data()), N * sizeof(Complex));
+ outfileCt.close();
+
+ std::ofstream outfileR(logFourierFile("R", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::out | std::ios::binary);
+ outfileR.write((const char*)(R.data()), N * sizeof(Real));
+ outfileR.close();
+
+ std::ofstream outfileRt(logFourierFile("Rt", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::out | std::ios::binary);
+ outfileRt.write((const char*)(Ȓ.data()), N * sizeof(Complex));
+ outfileRt.close();
+}
+
+bool logFourierLoad(std::vector<Real>& C, std::vector<Real>& R, std::vector<Complex>& Ĉ, std::vector<Complex>& Ȓ, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k) {
+ std::ifstream cfile(logFourierFile("C", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::binary);
+ std::ifstream rfile(logFourierFile("R", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::binary);
+ std::ifstream ctfile(logFourierFile("Ct", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::binary);
+ std::ifstream rtfile(logFourierFile("Rt", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::binary);
+
+ if ((!cfile.is_open() || !rfile.is_open()) || (!ctfile.is_open() || !rtfile.is_open())) {
+ return false;
+ }
+
+ unsigned N = std::pow(2, log2n);
+
+ cfile.read((char*)(C.data()), N * sizeof(Real));
+ cfile.close();
+
+ rfile.read((char*)(R.data()), N * sizeof(Real));
+ rfile.close();
+
+ ctfile.read((char*)(Ĉ.data()), N * sizeof(Complex));
+ ctfile.close();
+
+ rtfile.read((char*)(Ȓ.data()), N * sizeof(Complex));
+ rtfile.close();
+
+ return true;
+}
+
+std::tuple<std::vector<Complex>, std::vector<Complex>> ΣD(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, Real β, unsigned p, unsigned s, Real λ) {
+ std::vector<Real> D(C.size());
+ std::vector<Real> Σ(C.size());
+ Real β² = std::pow(β, 2);
+ for (unsigned n = 0; n < C.size(); n++) {
+ D[n] = β² * ∂f(λ, p, s, C[n]);
+ Σ[n] = β² * R[n] * ∂∂f(λ, p, s, C[n]);
+ }
+ std::vector<Complex> Σhat = fft.fourier(Σ, false);
+ std::vector<Complex> Dhat = fft.fourier(D, true);
+
+ return {Σhat, Dhat};
+}
+
+Real estimateZ(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ĉ, const std::vector<Real>& R, const std::vector<Complex>& Ȓ, unsigned p, unsigned s, Real λ, Real τ₀, Real β) {
+ auto [Σhat, Dhat] = ΣD(fft, C, R, β, p, s, λ);
+ Real Γ₀ = 1.0;
+
+ return (((2 * Γ₀ + Dhat[0]) * std::conj(Ȓ[0]) + Σhat[0] * Ĉ[0]) / Ĉ[0]).real();
+}
+
+Real energy(const LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real β) {
+ unsigned n₀ = 0;
+ for (unsigned n = 0; n < C.size(); n++) {
+ if (C[n] > 1 || R[n] > 1) n₀ = n % 2 == 0 ? n / 2 : (n + 1) / 2;
+ }
+ Real E = fft.t(2*n₀) * ∂f(λ, p, s, 1);
+ for (unsigned n = n₀; n < C.size()/2-1; n++) {
+ Real R₂ₙ = R[2*n];
+ Real R₂ₙ₊₁ = R[2*n+1];
+ Real R₂ₙ₊₂ = R[2*n+2];
+ Real C₂ₙ = C[2*n];
+ Real C₂ₙ₊₁ = C[2*n+1];
+ Real C₂ₙ₊₂ = C[2*n+2];
+
+ if (C₂ₙ₊₂ < 0 || R₂ₙ₊₂ < 0) break;
+
+ Real h₂ₙ = fft.t(2*n+1) - fft.t(2*n);
+ Real h₂ₙ₊₁ = fft.t(2*n+2) - fft.t(2*n+1);
+ Real f₂ₙ = R₂ₙ * ∂f(λ, p, s, C₂ₙ);
+ Real f₂ₙ₊₁ = R₂ₙ₊₁ * ∂f(λ, p, s, C₂ₙ₊₁);
+ Real f₂ₙ₊₂ = R₂ₙ₊₂ * ∂f(λ, p, s, C₂ₙ₊₂);
+
+ E += (h₂ₙ + h₂ₙ₊₁) / 6 * (
+ (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ
+ + std::pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
+ + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂
+ );
+ }
+ return β * E;
+}
+
+Real C0(const LogarithmicFourierTransform& fft, const std::vector<Complex>& Ĉ) {
+ Real C = 0;
+ for (unsigned n = 0; n < Ĉ.size()/2-1; n++) {
+ Real Ĉ₂ₙ = Ĉ[2*n].real();
+ Real Ĉ₂ₙ₊₁ = Ĉ[2*n+1].real();
+ Real Ĉ₂ₙ₊₂ = Ĉ[2*n+2].real();
+
+ Real h₂ₙ = fft.t(2*n+1) - fft.t(2*n);
+ Real h₂ₙ₊₁ = fft.t(2*n+2) - fft.t(2*n+1);
+ Real f₂ₙ = Ĉ₂ₙ;
+ Real f₂ₙ₊₁ = Ĉ₂ₙ₊₁;
+ Real f₂ₙ₊₂ = Ĉ₂ₙ₊₂;
+
+ C += (h₂ₙ + h₂ₙ₊₁) / 6 * (
+ (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ
+ + std::pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
+ + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂
+ );
+ }
+ return C * pow(fft.shift, 2) / M_PI;
+}
+
+void smooth(std::vector<Real>& C, Real ε) {
+ unsigned N = C.size();
+ bool trigger0 = false;
+ bool trigger1 = false;
+ for (unsigned i = 0; i < N; i++) {
+ if (C[i] < ε || trigger0) {
+ C[i] = 0;
+ trigger0 = true;
+ }
+ if (C[N-1-i] > 1 - ε || trigger1) {
+ C[N-1-i] = 1;
+ trigger1 = true;
+ }
+ }
+
+ Real Cmax = 0;
+ for (unsigned i = 0; i < N; i++) {
+ if (C[N-1-i] > Cmax) Cmax = C[N-1-i];
+ C[N-1-i] = Cmax;
+ }
+}
+
diff --git a/log-fourier.hpp b/log-fourier.hpp
new file mode 100644
index 0000000..755f7e9
--- /dev/null
+++ b/log-fourier.hpp
@@ -0,0 +1,58 @@
+#pragma once
+
+#include "types.hpp"
+
+#include <vector>
+#include <tuple>
+
+#include <fftw3.h>
+#include <gsl/gsl_sf_gamma.h>
+
+class LogarithmicFourierTransform {
+private:
+ Complex* a;
+ Complex* â;
+ FFTW_PLAN a_to_â;
+ FFTW_PLAN â_to_a;
+ unsigned N;
+ unsigned pad;
+ Real k;
+ Real Δτ;
+ Real τₛ;
+ Real ωₛ;
+ Real sₛ;
+ std::vector<Real> ts;
+ std::vector<Real> νs;
+ std::vector<Complex> Γs;
+ std::vector<Real> exp1kω;
+ std::vector<Real> exp1kτ;
+ std::vector<Real> expkω;
+ std::vector<Real> expkτ;
+public:
+ Real shift;
+ LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad = 4, Real shift = 0.5);
+ ~LogarithmicFourierTransform();
+ Real τ(unsigned n) const;
+ Real ω(unsigned n) const;
+ Real t(unsigned n) const;
+ Real ν(unsigned n) const;
+ Real s(unsigned n) const;
+ std::vector<Complex> fourier(const std::vector<Real>& c, bool symmetric);
+ std::vector<Real> inverse(const std::vector<Complex>& ĉ);
+};
+
+std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift);
+
+void logFourierSave(const std::vector<Real>& C, const std::vector<Real>& R, const std::vector<Complex>& Ĉ, const std::vector<Complex>& Ȓ, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift);
+
+bool logFourierLoad(std::vector<Real>& C, std::vector<Real>& R, std::vector<Complex>& Ĉ, std::vector<Complex>& Ȓ, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift);
+
+std::tuple<std::vector<Complex>, std::vector<Complex>> ΣD(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, Real β, unsigned p, unsigned s, Real λ);
+
+Real estimateZ(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ĉ, const std::vector<Real>& Ȓ, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β);
+
+Real energy(const LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real β);
+
+Real C0(const LogarithmicFourierTransform& fft, const std::vector<Complex>& Ĉ);
+
+void smooth(std::vector<Real>& C, Real ε);
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
new file mode 100644
index 0000000..e13cf1b
--- /dev/null
+++ b/log-fourier_integrator.cpp
@@ -0,0 +1,238 @@
+#include "log-fourier.hpp"
+
+#include <getopt.h>
+#include <iostream>
+#include <iomanip>
+
+int main(int argc, char* argv[]) {
+ /* Model parameters */
+ unsigned p = 2;
+ unsigned s = 2;
+ Real λ = 0.5;
+ Real τ₀ = 0;
+
+ /* Log-Fourier parameters */
+ unsigned log2n = 8;
+ Real Δτ = 0.1;
+ Real k = -0.01;
+ Real logShift = 0;
+ bool shiftSquare = false;
+
+ /* Iteration parameters */
+ Real ε = 1e-15;
+ Real γ₀ = 1;
+ Real x = 1;
+ Real β₀ = 0;
+ Real βₘₐₓ = 20;
+ Real Δβ = 0.01;
+ bool loadData = false;
+ unsigned pad = 2;
+
+ int opt;
+
+ while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:k:D:e:0:lg:x:P:h:S")) != -1) {
+ switch (opt) {
+ case 'p':
+ p = atoi(optarg);
+ break;
+ case 's':
+ s = atoi(optarg);
+ break;
+ case '2':
+ log2n = atoi(optarg);
+ break;
+ case 't':
+ τ₀ = atof(optarg);
+ break;
+ case 'b':
+ βₘₐₓ = atof(optarg);
+ break;
+ case 'd':
+ Δβ = atof(optarg);
+ break;
+ case 'g':
+ γ₀ = atof(optarg);
+ break;
+ case 'k':
+ k = atof(optarg);
+ break;
+ case 'h':
+ logShift = atof(optarg);
+ break;
+ case 'D':
+ Δτ = atof(optarg);
+ break;
+ case 'e':
+ ε = atof(optarg);
+ break;
+ case '0':
+ β₀ = atof(optarg);
+ break;
+ case 'x':
+ x = atof(optarg);
+ break;
+ case 'P':
+ pad = atoi(optarg);
+ break;
+ case 'l':
+ loadData = true;
+ break;
+ case 'S':
+ shiftSquare = true;
+ break;
+ default:
+ exit(1);
+ }
+ }
+
+ unsigned N = pow(2, log2n);
+
+ Real Γ₀ = 1;
+ Real μ₀ = τ₀ > 0 ? (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀) : Γ₀;
+
+ Real shift = μ₀ * pow(10, logShift);
+ if (shiftSquare) shift *= μ₀;
+ LogarithmicFourierTransform fft(N, k, Δτ, pad, shift);
+
+ std::cerr << "Starting, μ₀ = " << μ₀ << ", range " << fft.t(0) << " " << fft.t(N-1) << std::endl;
+
+ Real μₜ₋₁ = μ₀;
+
+ std::vector<Real> Cₜ₋₁(N);
+ std::vector<Real> Rₜ₋₁(N);
+ std::vector<Complex> Ĉₜ₋₁(N);
+ std::vector<Complex> Ȓₜ₋₁(N);
+ std::vector<Real> Γ(N);
+
+ for (unsigned n = 0; n < N; n++) {
+ Γ[n] = Γ₀ / (1 + std::pow(τ₀ * fft.ν(n), 2));
+ }
+
+ if (!loadData) {
+ /* Start from the exact solution for β = 0 */
+ for (unsigned n = 0; n < N; n++) {
+ if (τ₀ > 0) {
+ if (τ₀ == 2) {
+ Cₜ₋₁[n] = Γ₀ * std::exp(-fft.t(n) / 2) * (1 + fft.t(n) / 2);
+ } else {
+ Cₜ₋₁[n] = Γ₀ * (std::exp(-μ₀ * fft.t(n)) / μ₀ - τ₀ * std::exp(-fft.t(n) / τ₀)) / (1 - pow(μ₀ * τ₀, 2));
+ }
+ } else {
+ Cₜ₋₁[n] = Γ₀ * std::exp(-μ₀ * fft.t(n)) / μ₀;
+ }
+ Rₜ₋₁[n] = std::exp(-μ₀ * fft.t(n));
+
+ Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μ₀, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
+ Ȓₜ₋₁[n] = (Real)1 / (μ₀ + II * fft.ν(n));
+ }
+ } else {
+ if (!logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, logShift)) {
+ return 1;
+ }
+ μₜ₋₁ = estimateZ(fft, Cₜ₋₁, Ĉₜ₋₁, Rₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀);
+ }
+
+ std::vector<Real> Cₜ = Cₜ₋₁;
+ std::vector<Real> Rₜ = Rₜ₋₁;
+ std::vector<Complex> Ĉₜ = Ĉₜ₋₁;
+ std::vector<Complex> Ȓₜ = Ȓₜ₋₁;
+ Real μₜ = μₜ₋₁;
+
+ Real β = β₀ + Δβ;
+ while (β < βₘₐₓ) {
+ Real γ = γ₀;
+ Real ΔCₜ = 100;
+ Real ΔCₜ₋₁ = 101;
+ while (ΔCₜ > ε) {
+ auto [Σ, D] = ΣD(fft, Cₜ, Rₜ, β, p, s, λ);
+
+ std::vector<Complex> Ĉₜ₊₁(N);
+ std::vector<Complex> Ȓₜ₊₁(N);
+
+ Real C₀ = 0;
+ Real μ₊ = 0;
+ Real μ₋ = 0;
+
+ while (std::abs(C₀ - 1) > ε) {
+ for (unsigned n = 0; n < N; n++) {
+ Ĉₜ₊₁[n] = (((2 * Γ[n] + D[n]) * std::conj(Ȓₜ[n]) + Σ[n] * Ĉₜ[n]) / (μₜ + II * fft.ν(n))).real();
+ }
+ C₀ = C0(fft, Ĉₜ₊₁);
+ if (C₀ > 1) {
+ μ₋ = μₜ;
+ } else {
+ μ₊ = μₜ;
+ }
+ if (μ₋ > 0 && μ₊ > 0) {
+ μₜ = (μ₊ + μ₋) / 2;
+ } else {
+ μₜ *= pow(tanh(C₀-1)+1, x);
+ }
+ }
+
+ ΔCₜ = 0;
+ for (unsigned n = 0; n < N; n++) {
+ ΔCₜ += std::norm(((2 * Γ[n] + D[n]) * std::conj(Ȓₜ[n]) + Σ[n] * Ĉₜ[n]) - Ĉₜ[n] * (μₜ + II * fft.ν(n)));
+ ΔCₜ += std::norm(((Real)1 + Σ[n] * Ȓₜ[n]) - Ȓₜ[n] * (μₜ + II * fft.ν(n)));
+ }
+ ΔCₜ = sqrt(ΔCₜ) / (2*N);
+
+ for (unsigned n = 0; n < N; n++) {
+ Ȓₜ₊₁[n] = ((Real)1 + Σ[n] * Ȓₜ[n]) / (μₜ + II * fft.ν(n));
+ }
+
+ std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
+ std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
+
+ smooth(Cₜ₊₁, ε);
+ smooth(Rₜ₊₁, ε);
+
+ for (unsigned i = 0; i < N; i++) {
+ Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
+ Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]);
+ Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]);
+ Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
+ }
+
+ if (ΔCₜ > ΔCₜ₋₁ * 2 && ΔCₜ < 1e-10) {
+ γ = std::max(γ / 2, (Real)1e-2);
+ }
+
+ ΔCₜ₋₁ = ΔCₜ;
+
+ std::cerr << "\x1b[2K" << "\r";
+ std::cerr << std::setprecision(7) << β << " " << Δβ << " " << μₜ << " " << ΔCₜ << " " << γ;
+ }
+
+ if (std::isnan(ΔCₜ)) {
+ β -= Δβ;
+ Δβ *= 0.1;
+ β += Δβ;
+ Cₜ = Cₜ₋₁;
+ Rₜ = Rₜ₋₁;
+ Ĉₜ = Ĉₜ₋₁;
+ Ȓₜ = Ȓₜ₋₁;
+ μₜ = μₜ₋₁;
+ } else {
+ Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β);
+
+ std::cerr << "\x1b[2K" << "\r";
+ std::cerr << std::setprecision(7) << β << " " << Δβ << " " << μₜ << " " << Ĉₜ[0].real() << " " << E << std::endl;
+
+ logFourierSave(Cₜ, Rₜ, Ĉₜ, Ȓₜ, p, s, λ, τ₀, β, log2n, Δτ, logShift);
+
+ if (Ĉₜ[0].real() / Ĉₜ₋₁[0].real() > 1.25) {
+ Δβ = std::round(1e6 * Δβ / 2) / 1e6;
+ }
+
+ β = std::round(1e6 * (β + Δβ)) / 1e6;
+ Cₜ₋₁ = Cₜ;
+ Rₜ₋₁ = Rₜ;
+ Ĉₜ₋₁ = Ĉₜ;
+ Ȓₜ₋₁ = Ȓₜ;
+ μₜ₋₁ = μₜ;
+ }
+ }
+
+ return 0;
+}
diff --git a/log_get_energy.cpp b/log_get_energy.cpp
new file mode 100644
index 0000000..d156fd4
--- /dev/null
+++ b/log_get_energy.cpp
@@ -0,0 +1,100 @@
+#include "log-fourier.hpp"
+
+#include <getopt.h>
+#include <iostream>
+#include <iomanip>
+#include <filesystem>
+
+int main(int argc, char* argv[]) {
+ /* Model parameters */
+ unsigned p = 2;
+ unsigned s = 2;
+ Real λ = 0.5;
+ Real τ₀ = 0;
+
+ /* Log-Fourier parameters */
+ unsigned log2n = 8;
+ Real Δτ = 0.1;
+ Real k = 0.1;
+ unsigned pad = 2;
+ Real logShift = 0;
+ bool shiftSquare = false;
+
+ /* Iteration parameters */
+ Real β₀ = 0;
+ Real βₘₐₓ = 0.7;
+ Real Δβ = 0.01;
+
+ int opt;
+
+ while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:k:h:D:0:S")) != -1) {
+ switch (opt) {
+ case 'p':
+ p = atoi(optarg);
+ break;
+ case 's':
+ s = atoi(optarg);
+ break;
+ case '2':
+ log2n = atoi(optarg);
+ break;
+ case 't':
+ τ₀ = atof(optarg);
+ break;
+ case 'b':
+ βₘₐₓ = atof(optarg);
+ break;
+ case 'd':
+ Δβ = atof(optarg);
+ break;
+ case 'k':
+ k = atof(optarg);
+ break;
+ case 'h':
+ logShift = atof(optarg);
+ break;
+ case 'D':
+ Δτ = atof(optarg);
+ break;
+ case '0':
+ β₀ = atof(optarg);
+ break;
+ case 'S':
+ shiftSquare = true;
+ break;
+ default:
+ exit(1);
+ }
+ }
+
+ unsigned N = pow(2, log2n);
+ Real Γ₀ = 1;
+ Real μ₀ = τ₀ > 0 ? (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀) : Γ₀;
+
+ Real shift = μ₀ * pow(10, logShift);
+ if (shiftSquare) shift *= μ₀;
+ LogarithmicFourierTransform fft(N, k, Δτ, pad, shift);
+
+ std::vector<Real> C(N);
+ std::vector<Real> R(N);
+ std::vector<Complex> Ĉ(N);
+ std::vector<Complex> Ȓ(N);
+
+ Real β = β₀;
+
+ std::cout << std::setprecision(16);
+
+ while (β = std::round(1e6 * (β + Δβ)) / 1e6, β <= βₘₐₓ) {
+ if (std::filesystem::exists(logFourierFile("C", p, s, λ, τ₀, β, log2n, Δτ, logShift))) {
+ logFourierLoad(C, R, Ĉ, Ȓ, p, s, λ, τ₀, β, log2n, Δτ, logShift);
+
+ Real e = energy(fft, C, R, p, s, λ, β);
+
+ Real μ = estimateZ(fft, C, Ĉ, R, Ȓ, p, s, λ, τ₀, β);
+
+ std::cout << p << " " << s << " " << λ << " " << τ₀ << " " << β << " " << μ << " " << Ĉ[0].real() << " " << e << " " << Ȓ[0].real() << std::endl;
+ }
+ }
+
+ return 0;
+}
diff --git a/p-spin.cpp b/p-spin.cpp
new file mode 100644
index 0000000..ba18c6a
--- /dev/null
+++ b/p-spin.cpp
@@ -0,0 +1,25 @@
+#include "p-spin.hpp"
+
+inline Real fₚ(unsigned p, Real q) {
+ return 0.5 * pow(q, p);
+}
+
+inline Real ∂fₚ(unsigned p, Real q) {
+ return 0.5 * p * pow(q, p - 1);
+}
+
+inline Real ∂∂fₚ(unsigned p, Real q) {
+ return 0.5 * p * (p - 1) * pow(q, p - 2);
+}
+
+Real f(Real λ, unsigned p, unsigned s, Real q) {
+ return (1 - λ) * fₚ(p, q) + λ * fₚ(s, q);
+}
+
+Real ∂f(Real λ, unsigned p, unsigned s, Real q) {
+ return (1 - λ) * ∂fₚ(p, q) + λ * ∂fₚ(s, q);
+}
+
+Real ∂∂f(Real λ, unsigned p, unsigned s, Real q) {
+ return (1 - λ) * ∂∂fₚ(p, q) + λ * ∂∂fₚ(s, q);
+}
diff --git a/p-spin.hpp b/p-spin.hpp
new file mode 100644
index 0000000..484bc17
--- /dev/null
+++ b/p-spin.hpp
@@ -0,0 +1,7 @@
+#pragma once
+
+#include "types.hpp"
+
+Real f(Real λ, unsigned p, unsigned s, Real q);
+Real ∂f(Real λ, unsigned p, unsigned s, Real q);
+Real ∂∂f(Real λ, unsigned p, unsigned s, Real q);
diff --git a/print_correlations.cpp b/print_correlations.cpp
new file mode 100644
index 0000000..4627119
--- /dev/null
+++ b/print_correlations.cpp
@@ -0,0 +1,86 @@
+#include "log-fourier.hpp"
+
+#include <getopt.h>
+#include <iostream>
+#include <iomanip>
+#include <filesystem>
+
+int main(int argc, char* argv[]) {
+ /* Model parameters */
+ unsigned p = 2;
+ unsigned s = 2;
+ Real λ = 0.5;
+ Real τ₀ = 0;
+
+ /* Log-Fourier parameters */
+ unsigned log2n = 8;
+ Real Δτ = 0.1;
+ Real k = 0.1;
+ unsigned pad = 2;
+ Real logShift = 0;
+ bool shiftSquare = false;
+
+ /* Iteration parameters */
+ Real β = 0;
+
+ int opt;
+
+ while ((opt = getopt(argc, argv, "p:s:2:T:t:k:h:D:0:S")) != -1) {
+ switch (opt) {
+ case 'p':
+ p = atoi(optarg);
+ break;
+ case 's':
+ s = atoi(optarg);
+ break;
+ case '2':
+ log2n = atoi(optarg);
+ break;
+ case 't':
+ τ₀ = atof(optarg);
+ break;
+ case 'k':
+ k = atof(optarg);
+ break;
+ case 'h':
+ logShift = atof(optarg);
+ break;
+ case 'D':
+ Δτ = atof(optarg);
+ break;
+ case '0':
+ β = atof(optarg);
+ break;
+ case 'S':
+ shiftSquare = true;
+ break;
+ default:
+ exit(1);
+ }
+ }
+
+ unsigned N = pow(2, log2n);
+ Real Γ₀ = 1;
+ Real μ₀ = τ₀ > 0 ? (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀) : Γ₀;
+
+ Real shift = μ₀ * pow(10, logShift);
+ if (shiftSquare) shift *= μ₀;
+ LogarithmicFourierTransform fft(N, k, Δτ, pad, shift);
+
+ std::vector<Real> C(N);
+ std::vector<Real> R(N);
+ std::vector<Complex> Ĉ(N);
+ std::vector<Complex> Ȓ(N);
+
+ std::cout << std::setprecision(16);
+
+ if (std::filesystem::exists(logFourierFile("C", p, s, λ, τ₀, β, log2n, Δτ, logShift))) {
+ logFourierLoad(C, R, Ĉ, Ȓ, p, s, λ, τ₀, β, log2n, Δτ, logShift);
+
+ for (unsigned n = 0; n < N; n++) {
+ std::cout << fft.t(n) << " " << C[n] << " " << R[n] << " " << fft.ν(n) << " " << Ĉ[n].real() << " " << Ȓ[n].real() << " " << Ȓ[n].imag() << std::endl;
+ }
+ }
+
+ return 0;
+}
diff --git a/types.hpp b/types.hpp
new file mode 100644
index 0000000..71556fd
--- /dev/null
+++ b/types.hpp
@@ -0,0 +1,43 @@
+#pragma once
+#include <complex>
+#include <fftw3.h>
+
+using namespace std::complex_literals;
+
+#ifndef QUAD_PRECISION
+
+#define FFTW_PLAN fftw_plan
+#define FFTW_ALLOC_COMPLEX fftw_alloc_complex
+#define FFTW_IMPORT_WISDOM fftw_import_wisdom_from_filename
+#define FFTW_EXPORT_WISDOM fftw_export_wisdom_to_filename
+#define FFTW_PLAN_DFT_1D fftw_plan_dft_1d
+#define FFTW_DESTROY_PLAN fftw_destroy_plan
+#define FFTW_FREE fftw_free
+#define FFTW_CLEANUP fftw_cleanup
+#define FFTW_EXECUTE fftw_execute
+#define FFTW_COMPLEX fftw_complex
+
+using Real = double;
+
+#define II 1i
+
+#else
+
+#define FFTW_PLAN fftwl_plan
+#define FFTW_ALLOC_COMPLEX fftwl_alloc_complex
+#define FFTW_IMPORT_WISDOM fftwl_import_wisdom_from_filename
+#define FFTW_EXPORT_WISDOM fftwl_export_wisdom_to_filename
+#define FFTW_PLAN_DFT_1D fftwl_plan_dft_1d
+#define FFTW_DESTROY_PLAN fftwl_destroy_plan
+#define FFTW_FREE fftwl_free
+#define FFTW_CLEANUP fftwl_cleanup
+#define FFTW_EXECUTE fftwl_execute
+#define FFTW_COMPLEX fftwl_complex
+
+using Real = long double;
+
+#define II 1il
+
+#endif
+
+using Complex = std::complex<Real>;