Implemented binary search for μ

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++) {