From 1e8da9cdc0efdd675ca7e8e7b4eade7a2d7aa211 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sat, 19 Apr 2025 16:38:15 -0300
Subject: Implemented binary search for μ
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 log-fourier_integrator.cpp | 22 ++++++++++++++++++++--
 1 file 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++) {
-- 
cgit v1.2.3-70-g09d2