diff options
| -rw-r--r-- | log-fourier_integrator.cpp | 22 |
1 files changed, 20 insertions, 2 deletions
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index 3fb9d62..67ec7b4 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -87,6 +87,9 @@ int main(int argc, char* argv[]) { Real β = 0; while (β < βₘₐₓ) { + β += Δβ; + Real μ₁ = μ; + Real μ₂ = 0; while (true) { Real ΔC = 100; Real ΔC₀ = 100; @@ -142,8 +145,24 @@ int main(int argc, char* argv[]) { } if (std::abs(Cₜ[0] - 1) < ε) { break; + } + if (μ₂ < μ₁) { + /* Cₜ[0] will generally start too big. Until we find a value too small, we + * scale down proportional too Cₜ[0] */ + if (Cₜ[0] > 1) { + μ *= sqrt(sqrt(std::abs(Cₜ[0]))); + } else { /* Otherwise, we've found an upper bound for μ */ + μ₂ = μ; + } } else { - μ *= sqrt(sqrt(sqrt(std::abs(Cₜ[0])))); + /* Once μ₂ is set, we can use bisection */ + if (Cₜ[0] > 1) { + μ₁ = μ; + μ = μ₁ + (μ₂ - μ₁) / 2; + } else { + μ₂ = μ; + μ = μ₁ + (μ₂ - μ₁) / 2; + } } } @@ -164,7 +183,6 @@ int main(int argc, char* argv[]) { E *= β; std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl; - β += Δβ; } for (unsigned i = 0; i < N; i++) { |