From eb64a0681627b561478d37886d50d5a9e143a28e Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sat, 20 Jun 2026 15:26:34 -0700
Subject: [PATCH 01/22] Initial symplectic integrator

---
 .../quadrature/symplectic_integration.hpp     | 72 +++++++++++++++++++
 1 file changed, 72 insertions(+)
 create mode 100644 include/boost/math/quadrature/symplectic_integration.hpp

diff --git a/include/boost/math/quadrature/symplectic_integration.hpp b/include/boost/math/quadrature/symplectic_integration.hpp
new file mode 100644
index 0000000000..51ba121b6b
--- /dev/null
+++ b/include/boost/math/quadrature/symplectic_integration.hpp
@@ -0,0 +1,72 @@
+/* Based on "Construction of higher order symplectic integrators" by Haruo Yoshida*/
+
+#include <vector>
+#include <cmath>
+#include <valarray>
+
+typedef std::valarray<double> derivative(std::valarray<double>);
+
+template <class RealType>
+std::pair<std::valarray<RealType>, std::valarray<RealType> > second_order_yoshida(const std::valarray<RealType> p0, 
+                                                                                  const std::valarray<RealType> q0, 
+                                                                                  RealType dt, 
+                                                                                  derivative dHdp, 
+                                                                                  derivative dHdq)
+{
+    std::valarray<RealType> p = p0;
+    std::valarray<RealType> q = q0;
+
+    q = q + dt / 2 * dHdp(p);
+    p = p - dt * dHdq(q);
+    q = q + dt / 2 * dHdp(p);
+    
+    return std::make_pair(p, q);
+}
+
+template <class RealType>
+std::pair<std::valarray<RealType>, std::valarray<RealType> > fourth_order_yoshida(const std::valarray<RealType> p0, 
+                                                                                  const std::valarray<RealType> q0, 
+                                                                                  RealType dt, 
+                                                                                  derivative dHdp, 
+                                                                                  derivative dHdq)
+{
+    std::valarray<RealType> p = p0;
+    std::valarray<RealType> q = q0;
+
+    RealType x0 = -(std::pow(2, 1/3) / (2 - std::pow(2, 1/3)));
+    RealType x1 = 1 / (2 - std::pow(2, 1/3));
+
+    std::vector<RealType> weights = { x1, x0, x1 };
+
+    for (unsigned i=0; i < weights.size(); i++)
+    {
+        std::tie(p, q) = second_order_yoshida(p, q, weights[i] * dt, dHdp, dHdq);
+    }
+    
+    return std::make_pair(p, q);
+}
+
+template <class RealType>
+std::pair<std::valarray<RealType>, std::valarray<RealType> > sixth_order_yoshida(const std::valarray<RealType> p0, 
+                                                  const std::valarray<RealType> q0, 
+                                                  RealType dt, 
+                                                  derivative dHdp, 
+                                                  derivative dHdq)
+{
+    std::valarray<RealType> p = p0;
+    std::valarray<RealType> q = q0;
+    
+    // Choosing "System A" solution
+    RealType w3 = 0.784513610477560;
+    RealType w2 = 0.235573213359357;
+    RealType w1 = -1.17767998417887;
+    RealType w0 = 1 - 2 * (w1 + w2 + w3);
+    std::vector<RealType> weights = { w3, w2, w1, w0, w1, w2, w3};
+
+    for (unsigned i=0; i < weights.size(); i++)
+    {
+        std::tie(p, q) = second_order_yoshida(p, q, weights[i] * dt, dHdp, dHdq);
+    }
+    
+    return std::make_pair(p, q);
+}
\ No newline at end of file

From 62c80b7078c7e5a4595394e09f5d72b51b78cce2 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 21 Jun 2026 12:27:16 -0700
Subject: [PATCH 02/22] Generalized to floats

---
 .../quadrature/symplectic_integration.hpp     | 48 +++++++++----------
 1 file changed, 24 insertions(+), 24 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic_integration.hpp b/include/boost/math/quadrature/symplectic_integration.hpp
index 51ba121b6b..5bc50515cf 100644
--- a/include/boost/math/quadrature/symplectic_integration.hpp
+++ b/include/boost/math/quadrature/symplectic_integration.hpp
@@ -6,15 +6,15 @@
 
 typedef std::valarray<double> derivative(std::valarray<double>);
 
-template <class RealType>
-std::pair<std::valarray<RealType>, std::valarray<RealType> > second_order_yoshida(const std::valarray<RealType> p0, 
-                                                                                  const std::valarray<RealType> q0, 
-                                                                                  RealType dt, 
-                                                                                  derivative dHdp, 
-                                                                                  derivative dHdq)
+template <class RealType, class ReturnType, class Func>
+std::pair<ReturnType, ReturnType> second_order_yoshida(const ReturnType p0, 
+                                                       const ReturnType q0, 
+                                                       RealType dt, 
+                                                       Func dHdp, 
+                                                       Func dHdq)
 {
-    std::valarray<RealType> p = p0;
-    std::valarray<RealType> q = q0;
+    ReturnType p = p0;
+    ReturnType q = q0;
 
     q = q + dt / 2 * dHdp(p);
     p = p - dt * dHdq(q);
@@ -23,15 +23,15 @@ std::pair<std::valarray<RealType>, std::valarray<RealType> > second_order_yoshid
     return std::make_pair(p, q);
 }
 
-template <class RealType>
-std::pair<std::valarray<RealType>, std::valarray<RealType> > fourth_order_yoshida(const std::valarray<RealType> p0, 
-                                                                                  const std::valarray<RealType> q0, 
-                                                                                  RealType dt, 
-                                                                                  derivative dHdp, 
-                                                                                  derivative dHdq)
+template <class RealType, class ReturnType, class Func>
+std::pair<ReturnType, ReturnType> fourth_order_yoshida(const ReturnType p0, 
+                                                       const ReturnType q0, 
+                                                       const RealType dt, 
+                                                       Func dHdp, 
+                                                       Func dHdq)
 {
-    std::valarray<RealType> p = p0;
-    std::valarray<RealType> q = q0;
+    ReturnType p = p0;
+    ReturnType q = q0;
 
     RealType x0 = -(std::pow(2, 1/3) / (2 - std::pow(2, 1/3)));
     RealType x1 = 1 / (2 - std::pow(2, 1/3));
@@ -46,15 +46,15 @@ std::pair<std::valarray<RealType>, std::valarray<RealType> > fourth_order_yoshid
     return std::make_pair(p, q);
 }
 
-template <class RealType>
-std::pair<std::valarray<RealType>, std::valarray<RealType> > sixth_order_yoshida(const std::valarray<RealType> p0, 
-                                                  const std::valarray<RealType> q0, 
-                                                  RealType dt, 
-                                                  derivative dHdp, 
-                                                  derivative dHdq)
+template <class RealType, class ReturnType, class Func>
+std::pair<ReturnType, ReturnType> sixth_order_yoshida(const ReturnType p0, 
+                                                       const ReturnType q0, 
+                                                       RealType dt, 
+                                                       Func dHdp, 
+                                                       Func dHdq)
 {
-    std::valarray<RealType> p = p0;
-    std::valarray<RealType> q = q0;
+    ReturnType p = p0;
+    ReturnType q = q0;
     
     // Choosing "System A" solution
     RealType w3 = 0.784513610477560;

From 45ebaff28077aceea69c4f48c5fd75e33a420d37 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 21 Jun 2026 14:17:54 -0700
Subject: [PATCH 03/22] Added driver function

---
 include/boost/math/quadrature/symplectic.hpp  | 133 ++++++++++++++++++
 .../quadrature/symplectic_integration.hpp     |  72 ----------
 2 files changed, 133 insertions(+), 72 deletions(-)
 create mode 100644 include/boost/math/quadrature/symplectic.hpp
 delete mode 100644 include/boost/math/quadrature/symplectic_integration.hpp

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
new file mode 100644
index 0000000000..300c757713
--- /dev/null
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -0,0 +1,133 @@
+/* Based on "Construction of higher order symplectic integrators" by Haruo Yoshida*/
+
+#ifndef BOOST_MATH_QUADRATURE_SYMPLECTIC_HPP
+#define BOOST_MATH_QUADRATURE_SYMPLECTIC_HPP
+
+#include <vector>
+#include <cmath>
+#include <string>
+#include <algorithm>
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+
+namespace boost{ namespace math {namespace quadrature { 
+
+template <class RealType, class ReturnType, class Func>
+std::pair<ReturnType, ReturnType> second_order_yoshida(const ReturnType p0, 
+                                                       const ReturnType q0, 
+                                                       RealType dt, 
+                                                       Func dHdp, 
+                                                       Func dHdq)
+{
+    ReturnType p = p0;
+    ReturnType q = q0;
+
+    q = q + dt / 2 * dHdp(p);
+    p = p - dt * dHdq(q);
+    q = q + dt / 2 * dHdp(p);
+    
+    return std::make_pair(p, q);
+}
+
+template <class RealType, class ReturnType, class Func>
+std::pair<ReturnType, ReturnType> fourth_order_yoshida(const ReturnType p0, 
+                                                       const ReturnType q0, 
+                                                       const RealType dt, 
+                                                       Func dHdp, 
+                                                       Func dHdq)
+{
+    ReturnType p = p0;
+    ReturnType q = q0;
+
+    RealType x0 = -(std::pow(2, 1/3) / (2 - std::pow(2, 1/3)));
+    RealType x1 = 1 / (2 - std::pow(2, 1/3));
+
+    std::vector<RealType> weights = { x1, x0, x1 };
+
+    for (unsigned i=0; i < weights.size(); i++)
+    {
+        std::tie(p, q) = second_order_yoshida(p, q, weights[i] * dt, dHdp, dHdq);
+    }
+    
+    return std::make_pair(p, q);
+}
+
+template <class RealType, class ReturnType, class Func>
+std::pair<ReturnType, ReturnType> sixth_order_yoshida(const ReturnType p0, 
+                                                      const ReturnType q0, 
+                                                      RealType dt, 
+                                                      Func dHdp, 
+                                                      Func dHdq)
+{
+    ReturnType p = p0;
+    ReturnType q = q0;
+    
+    // Choosing "System A" solution
+    RealType w3 = 0.784513610477560;
+    RealType w2 = 0.235573213359357;
+    RealType w1 = -1.17767998417887;
+    RealType w0 = 1 - 2 * (w1 + w2 + w3);
+    std::vector<RealType> weights = { w3, w2, w1, w0, w1, w2, w3};
+
+    for (unsigned i=0; i < weights.size(); i++)
+    {
+        std::tie(p, q) = second_order_yoshida(p, q, weights[i] * dt, dHdp, dHdq);
+    }
+    
+    return std::make_pair(p, q);
+}
+
+template <class ReturnType, class RealType, class Func, class Policy>
+std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
+                                                                                   const ReturnType q0,
+                                                                                   const RealType dt,
+                                                                                   const unsigned steps,
+                                                                                   Func dHdp,
+                                                                                   Func dHdq,
+                                                                                   const Policy& pol,
+                                                                                   std::string method="Y6")
+{
+    // Not sure how to make this function string nicer
+    static const char* function = "boost::math::quadrature::integrate_hamiltonian(p0, q0, %1%, steps, dHdp, dHdq)";
+
+    if ((dt <= 0) || !(boost::math::isfinite)(dt))
+    {   // Need to figure out what the correct return type is here
+        RealType val = (boost::math::policies::raise_domain_error(function, "Time step must be positive and finite but got: dt = %1%.\n", dt, pol));
+        std::vector<ReturnType> nan_vec = {ReturnType(val)};
+        return std::make_pair(nan_vec, nan_vec);
+    }
+
+    // Check if method is available
+    std::vector<std::string> available_methods = {"Y6", "Y4", "Y2"};
+
+    std::pair<ReturnType, ReturnType> (*stepper)(ReturnType, ReturnType, RealType, Func, Func);
+
+    if (method == "Y6"){
+        stepper = sixth_order_yoshida;
+    } 
+    else if (method == "Y4"){
+        stepper = fourth_order_yoshida;
+    }
+    else if (method == "Y2")
+    {
+        stepper = second_order_yoshida;
+    }
+
+    std::vector<ReturnType> p(steps);
+    std::vector<ReturnType> q(steps);
+    p[0] = p0;
+    q[0] = q0;
+
+    ReturnType p_current = p0;
+    ReturnType q_current = q0;
+    for (unsigned i=1; i < steps; i++)
+    {
+        std::tie(p_current, q_current) = stepper(p_current, q_current, dt, dHdp, dHdq);
+        p[i] = p_current;
+        q[i] = q_current;
+    }
+    return std::make_pair(p, q);
+}
+
+}}}
+#endif
diff --git a/include/boost/math/quadrature/symplectic_integration.hpp b/include/boost/math/quadrature/symplectic_integration.hpp
deleted file mode 100644
index 5bc50515cf..0000000000
--- a/include/boost/math/quadrature/symplectic_integration.hpp
+++ /dev/null
@@ -1,72 +0,0 @@
-/* Based on "Construction of higher order symplectic integrators" by Haruo Yoshida*/
-
-#include <vector>
-#include <cmath>
-#include <valarray>
-
-typedef std::valarray<double> derivative(std::valarray<double>);
-
-template <class RealType, class ReturnType, class Func>
-std::pair<ReturnType, ReturnType> second_order_yoshida(const ReturnType p0, 
-                                                       const ReturnType q0, 
-                                                       RealType dt, 
-                                                       Func dHdp, 
-                                                       Func dHdq)
-{
-    ReturnType p = p0;
-    ReturnType q = q0;
-
-    q = q + dt / 2 * dHdp(p);
-    p = p - dt * dHdq(q);
-    q = q + dt / 2 * dHdp(p);
-    
-    return std::make_pair(p, q);
-}
-
-template <class RealType, class ReturnType, class Func>
-std::pair<ReturnType, ReturnType> fourth_order_yoshida(const ReturnType p0, 
-                                                       const ReturnType q0, 
-                                                       const RealType dt, 
-                                                       Func dHdp, 
-                                                       Func dHdq)
-{
-    ReturnType p = p0;
-    ReturnType q = q0;
-
-    RealType x0 = -(std::pow(2, 1/3) / (2 - std::pow(2, 1/3)));
-    RealType x1 = 1 / (2 - std::pow(2, 1/3));
-
-    std::vector<RealType> weights = { x1, x0, x1 };
-
-    for (unsigned i=0; i < weights.size(); i++)
-    {
-        std::tie(p, q) = second_order_yoshida(p, q, weights[i] * dt, dHdp, dHdq);
-    }
-    
-    return std::make_pair(p, q);
-}
-
-template <class RealType, class ReturnType, class Func>
-std::pair<ReturnType, ReturnType> sixth_order_yoshida(const ReturnType p0, 
-                                                       const ReturnType q0, 
-                                                       RealType dt, 
-                                                       Func dHdp, 
-                                                       Func dHdq)
-{
-    ReturnType p = p0;
-    ReturnType q = q0;
-    
-    // Choosing "System A" solution
-    RealType w3 = 0.784513610477560;
-    RealType w2 = 0.235573213359357;
-    RealType w1 = -1.17767998417887;
-    RealType w0 = 1 - 2 * (w1 + w2 + w3);
-    std::vector<RealType> weights = { w3, w2, w1, w0, w1, w2, w3};
-
-    for (unsigned i=0; i < weights.size(); i++)
-    {
-        std::tie(p, q) = second_order_yoshida(p, q, weights[i] * dt, dHdp, dHdq);
-    }
-    
-    return std::make_pair(p, q);
-}
\ No newline at end of file

From 34d3cd81d4d06ae9b38603d4fd1c25a6f8d2d3da Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 21 Jun 2026 17:12:35 -0700
Subject: [PATCH 04/22] Added error handling and default arguments

---
 include/boost/math/quadrature/symplectic.hpp | 47 ++++++++++++++++++--
 1 file changed, 43 insertions(+), 4 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index 300c757713..75e5ea614d 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -6,7 +6,8 @@
 #include <vector>
 #include <cmath>
 #include <string>
-#include <algorithm>
+#include <stdexcept>
+#include <boost/math/policies/policy.hpp>
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
 
@@ -84,14 +85,14 @@ std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltoni
                                                                                    const unsigned steps,
                                                                                    Func dHdp,
                                                                                    Func dHdq,
-                                                                                   const Policy& pol,
-                                                                                   std::string method="Y6")
+                                                                                   std::string method,
+                                                                                   const Policy& pol)
 {
     // Not sure how to make this function string nicer
     static const char* function = "boost::math::quadrature::integrate_hamiltonian(p0, q0, %1%, steps, dHdp, dHdq)";
 
     if ((dt <= 0) || !(boost::math::isfinite)(dt))
-    {   // Need to figure out what the correct return type is here
+    {
         RealType val = (boost::math::policies::raise_domain_error(function, "Time step must be positive and finite but got: dt = %1%.\n", dt, pol));
         std::vector<ReturnType> nan_vec = {ReturnType(val)};
         return std::make_pair(nan_vec, nan_vec);
@@ -112,6 +113,10 @@ std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltoni
     {
         stepper = second_order_yoshida;
     }
+    else
+    {
+        throw std::domain_error("Method must be in {'Y6', 'Y4', 'Y2'} but got: method = " + method);
+    }
 
     std::vector<ReturnType> p(steps);
     std::vector<ReturnType> q(steps);
@@ -129,5 +134,39 @@ std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltoni
     return std::make_pair(p, q);
 }
 
+template <class ReturnType, class RealType, class Func, class Policy>
+std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
+                                                                                   const ReturnType q0,
+                                                                                   const RealType dt,
+                                                                                   const unsigned steps,
+                                                                                   Func dHdp,
+                                                                                   Func dHdq,
+                                                                                   const Policy& pol)
+{
+    return integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, "Y6", pol);
+}
+
+template <class ReturnType, class RealType, class Func>
+std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
+                                                                                   const ReturnType q0,
+                                                                                   const RealType dt,
+                                                                                   const unsigned steps,
+                                                                                   Func dHdp,
+                                                                                   Func dHdq,
+                                                                                   std::string method)
+{
+    return integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, method, boost::math::policies::policy<>()); 
+}
+
+template <class ReturnType, class RealType, class Func>
+std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
+                                                                                   const ReturnType q0,
+                                                                                   const RealType dt,
+                                                                                   const unsigned steps,
+                                                                                   Func dHdp,
+                                                                                   Func dHdq)
+{
+    return integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, "Y6", boost::math::policies::policy<>()); 
+}
 }}}
 #endif

From 098ab7ddf424718c48641942b0c446ea62db1d59 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 21 Jun 2026 18:05:33 -0700
Subject: [PATCH 05/22] Initial testing

---
 test/Jamfile.v2          |  1 +
 test/test_symplectic.cpp | 66 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 67 insertions(+)
 create mode 100644 test/test_symplectic.cpp

diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 5b67e474df..7426c26fd6 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -1341,6 +1341,7 @@ test-suite quadrature :
    [ run test_trapezoidal.cpp /boost/test//boost_unit_test_framework : : :
       release [ requires cxx11_lambdas cxx11_auto_declarations cxx11_decltype cxx11_unified_initialization_syntax cxx11_variadic_templates ]
       $(float128_type) ]
+   [ run test_symplectic.cpp /boost/test//boost_unit_test_framework ]
 ;
 
 test-suite autodiff :
diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
new file mode 100644
index 0000000000..c3f2e06a12
--- /dev/null
+++ b/test/test_symplectic.cpp
@@ -0,0 +1,66 @@
+/*
+ * Copyright Nick Thompson, 2017
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+#define BOOST_TEST_MODULE symplectic_quadrature
+
+#include <valarray>
+#include <boost/math/tools/config.hpp>
+#include <boost/test/included/unit_test.hpp>
+#include <boost/test/tools/floating_point_comparison.hpp>
+#include <boost/math/tools/test_value.hpp>
+#include <boost/math/concepts/real_concept.hpp>
+#include <boost/math/quadrature/symplectic.hpp>
+#ifdef BOOST_HAS_FLOAT128
+#include <boost/multiprecision/complex128.hpp>
+#endif
+
+#if __has_include(<stdfloat>)
+#  include <stdfloat>
+#endif
+
+using boost::math::quadrature::integrate_hamiltonian;
+
+template <class Real>
+Real singleton_dHdp(Real p) { return p; }
+template <class Real>
+Real singleton_dHdq(Real q) { return q; }
+
+template <class Real>
+std::valarray<Real> vector_dHdp(std::valarray<Real> p) { return p; }
+template <class Real>
+std::valarray<Real> vector_dHdq(std::valarray<Real> q) { return q; }
+
+template<class RealType>
+void test_input_singleton_vs_vector()
+{
+    RealType dt = 0.05;
+    unsigned int steps = 3;
+    RealType q0 = 1;
+    RealType p0 = 0;
+
+
+    std::vector<RealType> p;
+    std::vector<RealType> q;
+
+    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, singleton_dHdp<RealType>, singleton_dHdq<RealType>);
+
+    std::valarray<RealType> q0_vector = {1};
+    std::valarray<RealType> p0_vector = {0};
+    std::vector<std::valarray<RealType> > p_vector;
+    std::vector<std::valarray<RealType> > q_vector;
+
+    std::tie(p_vector, q_vector) = boost::math::quadrature::integrate_hamiltonian(p0_vector, q0_vector, dt, steps, vector_dHdp<RealType>, vector_dHdq<RealType>);
+    
+    BOOST_CHECK_EQUAL(p[0], p_vector[0][0]);
+    BOOST_CHECK_EQUAL(p[1], p_vector[1][0]);
+    BOOST_CHECK_EQUAL(p[2], p_vector[2][0]);
+
+}
+
+BOOST_AUTO_TEST_CASE(symplectic_quadrature)
+{
+    test_input_singleton_vs_vector<double>();
+}

From fd343b732123f36e72e262de421620821a184b63 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Mon, 22 Jun 2026 17:18:16 -0700
Subject: [PATCH 06/22] Added SHO test

---
 include/boost/math/quadrature/symplectic.hpp | 12 ------
 test/test_symplectic.cpp                     | 41 ++++++++++++++++++++
 2 files changed, 41 insertions(+), 12 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index 75e5ea614d..12edb742a4 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -134,18 +134,6 @@ std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltoni
     return std::make_pair(p, q);
 }
 
-template <class ReturnType, class RealType, class Func, class Policy>
-std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
-                                                                                   const ReturnType q0,
-                                                                                   const RealType dt,
-                                                                                   const unsigned steps,
-                                                                                   Func dHdp,
-                                                                                   Func dHdq,
-                                                                                   const Policy& pol)
-{
-    return integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, "Y6", pol);
-}
-
 template <class ReturnType, class RealType, class Func>
 std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
                                                                                    const ReturnType q0,
diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index c3f2e06a12..252ddf5405 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -6,6 +6,7 @@
  */
 #define BOOST_TEST_MODULE symplectic_quadrature
 
+#include <algorithm>
 #include <valarray>
 #include <boost/math/tools/config.hpp>
 #include <boost/test/included/unit_test.hpp>
@@ -23,6 +24,7 @@
 
 using boost::math::quadrature::integrate_hamiltonian;
 
+// Equations of motion for simple harmonic oscillator 
 template <class Real>
 Real singleton_dHdp(Real p) { return p; }
 template <class Real>
@@ -33,6 +35,7 @@ std::valarray<Real> vector_dHdp(std::valarray<Real> p) { return p; }
 template <class Real>
 std::valarray<Real> vector_dHdq(std::valarray<Real> q) { return q; }
 
+// Sanity check that passing as a 1D vector does not effect results
 template<class RealType>
 void test_input_singleton_vs_vector()
 {
@@ -60,7 +63,45 @@ void test_input_singleton_vs_vector()
 
 }
 
+template <class RealType>
+void test_invalid_parameters()
+{
+    // Negative timestep
+    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(1, 0, -0.1, 10, singleton_dHdp<RealType>, singleton_dHdq<RealType>), std::domain_error);
+
+    // Method not in {'Y6', 'Y4', 'Y2'}
+    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(1, 0, 0.1, 10, singleton_dHdp<RealType>, singleton_dHdq<RealType>, "InvalidMethod"), std::domain_error);
+}
+
+/* Test if SHO energy fluctuations are below a given tolerance*/
+template <class RealType>
+void test_harmonic_oscillator(RealType tol)
+{
+    RealType dt = 0.05;
+    RealType t_end = 100;
+    unsigned int steps = t_end / dt;
+
+    RealType q0 = 1;
+    RealType p0 = 0;
+
+    std::vector<RealType> p;
+    std::vector<RealType> q;
+
+    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, singleton_dHdp<RealType>, singleton_dHdq<RealType>);
+
+    std::valarray<RealType> p_val(p.data(), p.size());
+    std::valarray<RealType> q_val(q.data(), q.size());
+    std::valarray<RealType> energy_error = std::pow(p_val, 2) + std::pow(q_val, 2) - 1;
+    std::valarray<RealType> abs_energy_error = std::abs(energy_error);
+    RealType max_error = abs_energy_error.max();
+    std::cout << max_error << std::endl;
+    BOOST_CHECK_LE(max_error, tol);
+
+}
+
 BOOST_AUTO_TEST_CASE(symplectic_quadrature)
 {
     test_input_singleton_vs_vector<double>();
+    test_invalid_parameters<double>();
+    test_harmonic_oscillator<double>(1e-10);
 }

From 63eac980d285c4c37654342e1ca41e603393c08d Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Mon, 22 Jun 2026 18:11:21 -0700
Subject: [PATCH 07/22] Switched to map functionality

---
 include/boost/math/quadrature/symplectic.hpp | 23 ++++++--------------
 test/test_symplectic.cpp                     |  2 +-
 2 files changed, 8 insertions(+), 17 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index 12edb742a4..1b270ed3c5 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -1,10 +1,11 @@
-/* Based on "Construction of higher order symplectic integrators" by Haruo Yoshida*/
+/* Copyright Jacob Hass, 2026*/
 
 #ifndef BOOST_MATH_QUADRATURE_SYMPLECTIC_HPP
 #define BOOST_MATH_QUADRATURE_SYMPLECTIC_HPP
 
 #include <vector>
 #include <cmath>
+#include <map>
 #include <string>
 #include <stdexcept>
 #include <boost/math/policies/policy.hpp>
@@ -101,22 +102,12 @@ std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltoni
     // Check if method is available
     std::vector<std::string> available_methods = {"Y6", "Y4", "Y2"};
 
-    std::pair<ReturnType, ReturnType> (*stepper)(ReturnType, ReturnType, RealType, Func, Func);
+    typedef std::pair<ReturnType, ReturnType> (*stepperType)(ReturnType, ReturnType, RealType, Func, Func);
 
-    if (method == "Y6"){
-        stepper = sixth_order_yoshida;
-    } 
-    else if (method == "Y4"){
-        stepper = fourth_order_yoshida;
-    }
-    else if (method == "Y2")
-    {
-        stepper = second_order_yoshida;
-    }
-    else
-    {
-        throw std::domain_error("Method must be in {'Y6', 'Y4', 'Y2'} but got: method = " + method);
-    }
+    std::map<std::string, stepperType> m{{"Y6", sixth_order_yoshida}, 
+                                         {"Y4", fourth_order_yoshida}, 
+                                         {"Y2", second_order_yoshida}};
+    stepperType stepper = m.at(method);
 
     std::vector<ReturnType> p(steps);
     std::vector<ReturnType> q(steps);
diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index 252ddf5405..f2608dd29c 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -70,7 +70,7 @@ void test_invalid_parameters()
     BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(1, 0, -0.1, 10, singleton_dHdp<RealType>, singleton_dHdq<RealType>), std::domain_error);
 
     // Method not in {'Y6', 'Y4', 'Y2'}
-    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(1, 0, 0.1, 10, singleton_dHdp<RealType>, singleton_dHdq<RealType>, "InvalidMethod"), std::domain_error);
+    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(1, 0, 0.1, 10, singleton_dHdp<RealType>, singleton_dHdq<RealType>, "InvalidMethod"), std::out_of_range);
 }
 
 /* Test if SHO energy fluctuations are below a given tolerance*/

From a3a916f117c930de4b3189377002af076354a785 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Wed, 24 Jun 2026 20:15:20 -0700
Subject: [PATCH 08/22] Added documentation

---
 doc/equations/harmonic_oscillator.png        | Bin 0 -> 1689 bytes
 doc/equations/harmonic_oscillator.svg        |  32 ++++++
 doc/html/quadrature.html                     |   1 +
 doc/math.qbk                                 |   1 +
 doc/quadrature/symplectic.qbk                | 103 +++++++++++++++++++
 include/boost/math/quadrature/symplectic.hpp |  12 ++-
 test/Jamfile.v2                              |   2 +-
 7 files changed, 147 insertions(+), 4 deletions(-)
 create mode 100644 doc/equations/harmonic_oscillator.png
 create mode 100644 doc/equations/harmonic_oscillator.svg
 create mode 100644 doc/quadrature/symplectic.qbk

diff --git a/doc/equations/harmonic_oscillator.png b/doc/equations/harmonic_oscillator.png
new file mode 100644
index 0000000000000000000000000000000000000000..427a3ea50877449cde9658f2b6f46ed8c86016c8
GIT binary patch
literal 1689
zcmV;K24?w*P)<h;3K|Lk000e1NJLTq004~u001ip1^@s6MH>rN00009a7bBm000ie
z000ie0hKEb8vp<XOG!jQRCwC$Txo1oMHn5mi>8(>!P3f7mtc*GCa5SBv{<7CRMfIa
zDppG=7)z}vF<McO5EY9C6$>gLxT1m@)QC%>Kh(rMMsSIW)#wjVBQX+}qTk6qH=R!J
zy?w9mwRc|MoaCgPx%a($XXcx4zGdbqm7ORQHVUs2whEUE4^j>e5x7lQtd-vsUhlfa
zAquYvZ_>&ugfChAz5|2@yJ3qL`N9LV@^iwK7QW|T;Y{I|!ijF!?7;-#hr&UYz1QW!
zTZL7^|Adp=P}zae!tKJtEPlt4PH5~wiEx$h5SGt!LX*9~;4c==r`-beekE1CuTk*#
zdEqSKbY5TPgeE(I_{-JGi@E$m*LoDMjT62@L0T$&)nJtqn(V}E3P}k&x?H<9kab6A
zSeu;C^fG*bm?8YpxONgde}WSlhajvGK5o>_Vds}Pp>YVpPGP-K_q6c63~WoBLNkgE
zMYC|4>KqX5S11zh7nT_Y@*l#fPG}OwaN$A<N)9JoOMev_gx&u(jNp5PcW1Dj2BE3!
z<LAN)eL*Y37ThPy3;V$F!m$Pi^4gGqYh@O7@&lX|O%$Lzl(!8c(^Tt4!qvi8gkLK!
z6OZ!Tko(cTiEQ%)X?T9Ap_wZwG`V3PK+NJHO88SWB8xQ1ucvqND^=KJi#mH3XNm7-
zADAU}95euel-~*;>4ooqQh5Q=NHkDr4iAy)pM~oz>$qatGbdQ^E;%th_j}<Ly|7f>
zmx#80R`^B8wJ{W^*_H`S4Yi*Q+65D8XiWCc5oul~uhpE4(B;Y3Bx^&iH%4foxB%A>
z8>x2iy+_3SJ{wiv$An<-w59R;DpW^$1v0TDYC?9WNw!63A&2;cX0u_>6jS@28nAR;
z$ZP&SqSJ-7bn$`l!oAeoV7v169fe{cg$7>KV1w|U_C#nQhxmj>w`az4(pRwF-vWYN
zi;Jm=5xo0=9@lBg3*hbMeHiN_4GxQX_*UT$!rKG*Sd#V()Y&ttYgRdmE1b|IjiWSs
z1~HYhS&!`>8#E)fkrUqsyx`^hf1c`L0MEB))^wwV;?Cg|1ZctD;DA0Bn*0Rc%RF*Z
zOD!>_n-GY;1!O44a$Kr|aU#EHW`tMuxdk(M8I9iFW_K}<1tHoJGx)85`t8O{H*s--
z?{#c|J%bp{7LHwR0)2?!-S6oL4>U$-5?>3oCi$JmsjgAvD3-eADQUFPu+5{rxg&wk
z-~`5t>3BFiqg~y2RhK65eXwbkhS*18oacn5NA`?A*-mD;faP;TK8LLTEn0abSDK`N
z0sbYzGip$oxhJK+3r#JTCLgeAe$xbDhUQUi<#Ghi>0jjq?|dMLm++yyeERJxC=}bc
zw&rSf4)RZ*DKCRu!zIVT4xz!b^pp-Z&VNTuvxwsh+`wh4fNnQ<BZa6sAf6Km`0$lf
zYh!6mPuAQt2uQhRbDTlH-4?)`IFGT0!hR$)uump&Y>Fc8s|c_Ij}g`xu@wl%CCWR9
zcPQ)5QC=i@L^jRnIgB$*rFv~s9l4S?(UhV9V*=k#jStQKiSl+H6e@?5ForI#4ur$^
z2h?X0p!r)==d3c~ML#KTuaH9uSRLYNAvC@-gWCn;^uYkFau+8=u^alH1GBC@pmXsn
z5PLJYUGwRXoNv_o(fCu`(Dxj}xb~a)h}y=-G<pjD(c9@<zF>FUqr8xO8SCb9mXD6H
zUHLR_%~U}PC&>^u^yx(zRd_Ya=-_&K4(;rCJzVgv@`8_bTHOL_E1Qr*H~L{W#OVie
zJL@$gvW6O|jH-Q#asQ{v+Xt=V8sXcPwQCVQs6lRs)8DaMRObYoAE&rX_)UP9@|W`B
zPC9Tpv(=(@A^!BSGl0_%9NXBa!ZWrTsDX5xDg0A)L`(rh!S~P(vH=L9zO4Fwzj40@
z&ikx1;u__~J-d|`3DtpQ-RnjjBFo6zAP1AD;{Aa{qysh^OH!%wPCq)N02sB+s_XHJ
zW}~C^EU=F^gf|;?>okpoY<^tHS%mpl1%wmol-)n~W<7?}I9zM+3E?K~zG~I6_f>am
z_o=SGnh9@e7Tw*ZnA(LU68l^pHQKeF?R7{2cM9JRVS{$x65gu=XN0~Auz`@vLeA&9
j7(JD$h<0s?yVL&%h;d{n+q~X_00000NkvXXu0mjfKGipv

literal 0
HcmV?d00001

diff --git a/doc/equations/harmonic_oscillator.svg b/doc/equations/harmonic_oscillator.svg
new file mode 100644
index 0000000000..7320acb737
--- /dev/null
+++ b/doc/equations/harmonic_oscillator.svg
@@ -0,0 +1,32 @@
+<?xml version='1.0'?>
+<!-- This file was generated by dvisvgm 1.14.1 -->
+<svg height='32.5582pt' version='1.1' viewBox='238.04 98.1763 105.382 32.5582' width='105.382pt' xmlns='http://www.w3.org/2000/svg' xmlns:xlink='http://www.w3.org/1999/xlink'>
+<defs>
+<path d='M4.77011 -2.76164H8.06974C8.23711 -2.76164 8.4523 -2.76164 8.4523 -2.97684C8.4523 -3.20399 8.24907 -3.20399 8.06974 -3.20399H4.77011V-6.50361C4.77011 -6.67098 4.77011 -6.88618 4.55492 -6.88618C4.32777 -6.88618 4.32777 -6.68294 4.32777 -6.50361V-3.20399H1.02814C0.860772 -3.20399 0.645579 -3.20399 0.645579 -2.98879C0.645579 -2.76164 0.848817 -2.76164 1.02814 -2.76164H4.32777V0.537983C4.32777 0.705355 4.32777 0.920548 4.54296 0.920548C4.77011 0.920548 4.77011 0.71731 4.77011 0.537983V-2.76164Z' id='g2-43'/>
+<path d='M3.44309 -7.66326C3.44309 -7.93823 3.44309 -7.95019 3.20399 -7.95019C2.91706 -7.6274 2.3193 -7.18506 1.08792 -7.18506V-6.83836C1.36289 -6.83836 1.96065 -6.83836 2.61818 -7.14919V-0.920548C2.61818 -0.490162 2.58232 -0.3467 1.53026 -0.3467H1.15965V0C1.48244 -0.0239103 2.64209 -0.0239103 3.03661 -0.0239103S4.57883 -0.0239103 4.90162 0V-0.3467H4.53101C3.47895 -0.3467 3.44309 -0.490162 3.44309 -0.920548V-7.66326Z' id='g2-49'/>
+<path d='M5.26027 -2.00847H4.99726C4.96139 -1.80523 4.86575 -1.1477 4.7462 -0.956413C4.66252 -0.848817 3.98107 -0.848817 3.62242 -0.848817H1.41071C1.7335 -1.12379 2.46276 -1.88892 2.7736 -2.17584C4.59078 -3.84956 5.26027 -4.47123 5.26027 -5.65479C5.26027 -7.02964 4.17235 -7.95019 2.78555 -7.95019S0.585803 -6.76663 0.585803 -5.73848C0.585803 -5.12877 1.11183 -5.12877 1.1477 -5.12877C1.39875 -5.12877 1.70959 -5.30809 1.70959 -5.69066C1.70959 -6.0254 1.48244 -6.25255 1.1477 -6.25255C1.0401 -6.25255 1.01619 -6.25255 0.980324 -6.2406C1.20747 -7.05355 1.85305 -7.60349 2.63014 -7.60349C3.64633 -7.60349 4.268 -6.75467 4.268 -5.65479C4.268 -4.63861 3.68219 -3.75392 3.00075 -2.98879L0.585803 -0.286924V0H4.94944L5.26027 -2.00847Z' id='g2-50'/>
+<path d='M8.06974 -3.87347C8.23711 -3.87347 8.4523 -3.87347 8.4523 -4.08867C8.4523 -4.31582 8.24907 -4.31582 8.06974 -4.31582H1.02814C0.860772 -4.31582 0.645579 -4.31582 0.645579 -4.10062C0.645579 -3.87347 0.848817 -3.87347 1.02814 -3.87347H8.06974ZM8.06974 -1.64981C8.23711 -1.64981 8.4523 -1.64981 8.4523 -1.86501C8.4523 -2.09215 8.24907 -2.09215 8.06974 -2.09215H1.02814C0.860772 -2.09215 0.645579 -2.09215 0.645579 -1.87696C0.645579 -1.64981 0.848817 -1.64981 1.02814 -1.64981H8.06974Z' id='g2-61'/>
+<path d='M8.94247 -7.29265C9.05006 -7.69913 9.07397 -7.81868 9.92279 -7.81868C10.138 -7.81868 10.2575 -7.81868 10.2575 -8.03387C10.2575 -8.16538 10.1499 -8.16538 10.0782 -8.16538C9.86301 -8.16538 9.61196 -8.14147 9.38481 -8.14147H7.9741C7.74695 -8.14147 7.49589 -8.16538 7.26874 -8.16538C7.18506 -8.16538 7.04159 -8.16538 7.04159 -7.93823C7.04159 -7.81868 7.12528 -7.81868 7.35243 -7.81868C8.06974 -7.81868 8.06974 -7.72304 8.06974 -7.59153C8.06974 -7.56762 8.06974 -7.49589 8.02192 -7.31656L7.29265 -4.42341H3.68219L4.3995 -7.29265C4.5071 -7.69913 4.53101 -7.81868 5.37983 -7.81868C5.59502 -7.81868 5.71457 -7.81868 5.71457 -8.03387C5.71457 -8.16538 5.60697 -8.16538 5.53524 -8.16538C5.32005 -8.16538 5.06899 -8.14147 4.84184 -8.14147H3.43113C3.20399 -8.14147 2.95293 -8.16538 2.72578 -8.16538C2.64209 -8.16538 2.49863 -8.16538 2.49863 -7.93823C2.49863 -7.81868 2.58232 -7.81868 2.80946 -7.81868C3.52677 -7.81868 3.52677 -7.72304 3.52677 -7.59153C3.52677 -7.56762 3.52677 -7.49589 3.47895 -7.31656L1.86501 -0.884682C1.75741 -0.466252 1.7335 -0.3467 0.908593 -0.3467C0.633624 -0.3467 0.549938 -0.3467 0.549938 -0.119552C0.549938 0 0.681445 0 0.71731 0C0.932503 0 1.18356 -0.0239103 1.41071 -0.0239103H2.82142C3.04857 -0.0239103 3.29963 0 3.52677 0C3.62242 0 3.75392 0 3.75392 -0.227148C3.75392 -0.3467 3.64633 -0.3467 3.467 -0.3467C2.73773 -0.3467 2.73773 -0.442341 2.73773 -0.561893C2.73773 -0.573848 2.73773 -0.657534 2.76164 -0.753176L3.58655 -4.07671H7.20897C7.00573 -3.28767 6.39601 -0.789041 6.3721 -0.71731C6.2406 -0.358655 6.04932 -0.358655 5.34396 -0.3467C5.2005 -0.3467 5.0929 -0.3467 5.0929 -0.119552C5.0929 0 5.22441 0 5.26027 0C5.47547 0 5.72653 -0.0239103 5.95367 -0.0239103H7.36438C7.59153 -0.0239103 7.84259 0 8.06974 0C8.16538 0 8.29689 0 8.29689 -0.227148C8.29689 -0.3467 8.18929 -0.3467 8.00996 -0.3467C7.2807 -0.3467 7.2807 -0.442341 7.2807 -0.561893C7.2807 -0.573848 7.2807 -0.657534 7.30461 -0.753176L8.94247 -7.29265Z' id='g0-72'/>
+<path d='M3.3594 -7.99801C3.37136 -8.04583 3.39527 -8.11756 3.39527 -8.17733C3.39527 -8.29689 3.27572 -8.29689 3.25181 -8.29689C3.23985 -8.29689 2.80946 -8.26102 2.59427 -8.23711C2.39103 -8.22516 2.21171 -8.20125 1.99651 -8.18929C1.70959 -8.16538 1.6259 -8.15342 1.6259 -7.93823C1.6259 -7.81868 1.74545 -7.81868 1.86501 -7.81868C2.47472 -7.81868 2.47472 -7.71108 2.47472 -7.59153C2.47472 -7.54371 2.47472 -7.5198 2.41494 -7.30461L0.705355 -0.466252C0.657534 -0.286924 0.657534 -0.263014 0.657534 -0.191283C0.657534 0.071731 0.860772 0.119552 0.980324 0.119552C1.31507 0.119552 1.3868 -0.143462 1.48244 -0.514072L2.04433 -2.74969C2.90511 -2.65405 3.41918 -2.29539 3.41918 -1.72154C3.41918 -1.64981 3.41918 -1.60199 3.38331 -1.42267C3.33549 -1.24334 3.33549 -1.09988 3.33549 -1.0401C3.33549 -0.3467 3.78979 0.119552 4.3995 0.119552C4.94944 0.119552 5.23636 -0.382565 5.332 -0.549938C5.58306 -0.992279 5.73848 -1.66177 5.73848 -1.70959C5.73848 -1.76936 5.69066 -1.81719 5.61893 -1.81719C5.51133 -1.81719 5.49938 -1.76936 5.45156 -1.57808C5.28418 -0.956413 5.03313 -0.119552 4.42341 -0.119552C4.18431 -0.119552 4.02889 -0.239103 4.02889 -0.6934C4.02889 -0.920548 4.07671 -1.18356 4.12453 -1.36289C4.17235 -1.57808 4.17235 -1.59004 4.17235 -1.7335C4.17235 -2.43885 3.53873 -2.83337 2.43885 -2.97684C2.86924 -3.23985 3.29963 -3.7061 3.467 -3.88543C4.14844 -4.65056 4.61469 -5.03313 5.16463 -5.03313C5.4396 -5.03313 5.51133 -4.96139 5.59502 -4.88966C5.15268 -4.84184 4.98531 -4.53101 4.98531 -4.29191C4.98531 -4.00498 5.21245 -3.90934 5.37983 -3.90934C5.70262 -3.90934 5.98954 -4.18431 5.98954 -4.56687C5.98954 -4.91357 5.71457 -5.27223 5.17659 -5.27223C4.51905 -5.27223 3.98107 -4.80598 3.13225 -3.84956C3.0127 -3.7061 2.57036 -3.25181 2.12802 -3.08443L3.3594 -7.99801Z' id='g0-107'/>
+<path d='M2.46276 -3.50286C2.48667 -3.5746 2.78555 -4.17235 3.2279 -4.55492C3.53873 -4.84184 3.94521 -5.03313 4.41146 -5.03313C4.88966 -5.03313 5.05704 -4.67447 5.05704 -4.19626C5.05704 -4.12453 5.05704 -3.88543 4.91357 -3.32354L4.61469 -2.09215C4.51905 -1.7335 4.29191 -0.848817 4.268 -0.71731C4.22017 -0.537983 4.14844 -0.227148 4.14844 -0.179328C4.14844 -0.0119552 4.27995 0.119552 4.45928 0.119552C4.81793 0.119552 4.87771 -0.155417 4.98531 -0.585803L5.70262 -3.44309C5.72653 -3.53873 6.34819 -5.03313 7.66326 -5.03313C8.14147 -5.03313 8.30884 -4.67447 8.30884 -4.19626C8.30884 -3.52677 7.84259 -2.22366 7.57958 -1.50635C7.47198 -1.21943 7.4122 -1.06401 7.4122 -0.848817C7.4122 -0.310834 7.78281 0.119552 8.35666 0.119552C9.46849 0.119552 9.88692 -1.63786 9.88692 -1.70959C9.88692 -1.76936 9.8391 -1.81719 9.76737 -1.81719C9.65978 -1.81719 9.64782 -1.78132 9.58804 -1.57808C9.31308 -0.621669 8.87073 -0.119552 8.39253 -0.119552C8.27298 -0.119552 8.08169 -0.131507 8.08169 -0.514072C8.08169 -0.824907 8.22516 -1.20747 8.27298 -1.33898C8.48817 -1.91283 9.02615 -3.32354 9.02615 -4.01694C9.02615 -4.73425 8.60772 -5.27223 7.69913 -5.27223C6.89813 -5.27223 6.25255 -4.81793 5.77435 -4.11258C5.73848 -4.75816 5.34396 -5.27223 4.44732 -5.27223C3.38331 -5.27223 2.82142 -4.51905 2.60623 -4.22017C2.57036 -4.90162 2.0802 -5.27223 1.55417 -5.27223C1.20747 -5.27223 0.932503 -5.10486 0.705355 -4.65056C0.490162 -4.22017 0.32279 -3.49091 0.32279 -3.44309S0.37061 -3.33549 0.454296 -3.33549C0.549938 -3.33549 0.561893 -3.34745 0.633624 -3.62242C0.812951 -4.32777 1.0401 -5.03313 1.51831 -5.03313C1.79328 -5.03313 1.88892 -4.84184 1.88892 -4.48319C1.88892 -4.22017 1.76936 -3.75392 1.68568 -3.38331L1.35093 -2.09215C1.30311 -1.86501 1.17161 -1.32702 1.11183 -1.11183C1.02814 -0.800996 0.896638 -0.239103 0.896638 -0.179328C0.896638 -0.0119552 1.02814 0.119552 1.20747 0.119552C1.35093 0.119552 1.51831 0.0478207 1.61395 -0.131507C1.63786 -0.191283 1.74545 -0.609714 1.80523 -0.848817L2.06824 -1.92478L2.46276 -3.50286Z' id='g0-109'/>
+<path d='M0.514072 1.51831C0.430386 1.87696 0.382565 1.9726 -0.107597 1.9726C-0.251059 1.9726 -0.37061 1.9726 -0.37061 2.19975C-0.37061 2.22366 -0.358655 2.3193 -0.227148 2.3193C-0.071731 2.3193 0.0956413 2.29539 0.251059 2.29539H0.765131C1.01619 2.29539 1.6259 2.3193 1.87696 2.3193C1.94869 2.3193 2.09215 2.3193 2.09215 2.10411C2.09215 1.9726 2.00847 1.9726 1.80523 1.9726C1.25529 1.9726 1.21943 1.88892 1.21943 1.79328C1.21943 1.64981 1.75741 -0.406476 1.82914 -0.681445C1.96065 -0.3467 2.28344 0.119552 2.90511 0.119552C4.25604 0.119552 5.71457 -1.63786 5.71457 -3.39527C5.71457 -4.49514 5.0929 -5.27223 4.19626 -5.27223C3.43113 -5.27223 2.78555 -4.53101 2.65405 -4.36364C2.55841 -4.96139 2.09215 -5.27223 1.61395 -5.27223C1.26725 -5.27223 0.992279 -5.10486 0.765131 -4.65056C0.549938 -4.22017 0.382565 -3.49091 0.382565 -3.44309S0.430386 -3.33549 0.514072 -3.33549C0.609714 -3.33549 0.621669 -3.34745 0.6934 -3.62242C0.872727 -4.32777 1.09988 -5.03313 1.57808 -5.03313C1.85305 -5.03313 1.94869 -4.84184 1.94869 -4.48319C1.94869 -4.19626 1.91283 -4.07671 1.86501 -3.86152L0.514072 1.51831ZM2.58232 -3.73001C2.666 -4.06476 3.00075 -4.41146 3.19203 -4.57883C3.32354 -4.69838 3.71806 -5.03313 4.17235 -5.03313C4.69838 -5.03313 4.93748 -4.5071 4.93748 -3.88543C4.93748 -3.31158 4.60274 -1.96065 4.30386 -1.33898C4.00498 -0.6934 3.45504 -0.119552 2.90511 -0.119552C2.09215 -0.119552 1.96065 -1.1477 1.96065 -1.19552C1.96065 -1.23138 1.98456 -1.32702 1.99651 -1.3868L2.58232 -3.73001Z' id='g0-112'/>
+<path d='M5.66675 -4.87771C5.28418 -4.80598 5.14072 -4.51905 5.14072 -4.29191C5.14072 -4.00498 5.36787 -3.90934 5.53524 -3.90934C5.8939 -3.90934 6.14496 -4.22017 6.14496 -4.54296C6.14496 -5.04508 5.57111 -5.27223 5.06899 -5.27223C4.33973 -5.27223 3.93325 -4.55492 3.82565 -4.32777C3.55068 -5.22441 2.80946 -5.27223 2.59427 -5.27223C1.37484 -5.27223 0.729265 -3.7061 0.729265 -3.44309C0.729265 -3.39527 0.777086 -3.33549 0.860772 -3.33549C0.956413 -3.33549 0.980324 -3.40722 1.00423 -3.45504C1.41071 -4.78207 2.21171 -5.03313 2.55841 -5.03313C3.09639 -5.03313 3.20399 -4.53101 3.20399 -4.24408C3.20399 -3.98107 3.13225 -3.7061 2.98879 -3.13225L2.58232 -1.4944C2.40299 -0.777086 2.05629 -0.119552 1.42267 -0.119552C1.36289 -0.119552 1.06401 -0.119552 0.812951 -0.274969C1.24334 -0.358655 1.33898 -0.71731 1.33898 -0.860772C1.33898 -1.09988 1.15965 -1.24334 0.932503 -1.24334C0.645579 -1.24334 0.334745 -0.992279 0.334745 -0.609714C0.334745 -0.107597 0.896638 0.119552 1.41071 0.119552C1.98456 0.119552 2.39103 -0.334745 2.64209 -0.824907C2.83337 -0.119552 3.43113 0.119552 3.87347 0.119552C5.0929 0.119552 5.73848 -1.44658 5.73848 -1.70959C5.73848 -1.76936 5.69066 -1.81719 5.61893 -1.81719C5.51133 -1.81719 5.49938 -1.75741 5.46351 -1.66177C5.14072 -0.609714 4.44732 -0.119552 3.90934 -0.119552C3.49091 -0.119552 3.26376 -0.430386 3.26376 -0.920548C3.26376 -1.18356 3.31158 -1.37484 3.50286 -2.16389L3.9213 -3.78979C4.10062 -4.5071 4.5071 -5.03313 5.05704 -5.03313C5.08095 -5.03313 5.41569 -5.03313 5.66675 -4.87771Z' id='g0-120'/>
+<path d='M2.24757 -1.6259C2.37509 -1.74545 2.70984 -2.00847 2.83736 -2.12005C3.33151 -2.57435 3.80174 -3.0127 3.80174 -3.73798C3.80174 -4.68643 3.00473 -5.30012 2.00847 -5.30012C1.05205 -5.30012 0.422416 -4.57484 0.422416 -3.8655C0.422416 -3.47497 0.73325 -3.41918 0.844832 -3.41918C1.0122 -3.41918 1.25928 -3.53873 1.25928 -3.84159C1.25928 -4.25604 0.860772 -4.25604 0.765131 -4.25604C0.996264 -4.83786 1.53026 -5.03711 1.9208 -5.03711C2.66202 -5.03711 3.04458 -4.40747 3.04458 -3.73798C3.04458 -2.90909 2.46276 -2.30336 1.52229 -1.33898L0.518057 -0.302864C0.422416 -0.215193 0.422416 -0.199253 0.422416 0H3.57061L3.80174 -1.42665H3.55467C3.53076 -1.26725 3.467 -0.868742 3.37136 -0.71731C3.32354 -0.653549 2.71781 -0.653549 2.59029 -0.653549H1.17161L2.24757 -1.6259Z' id='g1-50'/>
+</defs>
+<g id='page1' transform='matrix(1.25 0 0 1.25 0 0)'>
+<use x='190.432' xlink:href='#g0-72' y='96.2674'/>
+<use x='204.339' xlink:href='#g2-61' y='96.2674'/>
+<use x='220.668' xlink:href='#g0-112' y='88.1796'/>
+<use x='226.522' xlink:href='#g1-50' y='83.8412'/>
+<rect height='0.478187' width='16.0923' x='217.926' y='93.0395'/>
+<use x='217.926' xlink:href='#g2-50' y='104.468'/>
+<use x='223.757' xlink:href='#g0-109' y='104.468'/>
+<use x='237.87' xlink:href='#g2-43' y='96.2674'/>
+<use x='250.793' xlink:href='#g2-49' y='88.1796'/>
+<rect height='0.478187' width='5.85299' x='250.793' y='93.0395'/>
+<use x='250.793' xlink:href='#g2-50' y='104.468'/>
+<use x='257.842' xlink:href='#g0-107' y='96.2674'/>
+<use x='264.308' xlink:href='#g0-120' y='96.2674'/>
+<use x='270.936' xlink:href='#g1-50' y='91.3312'/>
+</g>
+</svg>
diff --git a/doc/html/quadrature.html b/doc/html/quadrature.html
index 9041b488ae..dadad954c4 100644
--- a/doc/html/quadrature.html
+++ b/doc/html/quadrature.html
@@ -51,6 +51,7 @@
 </dl></dd>
 <dt><span class="section"><a href="math_toolkit/fourier_integrals.html">Fourier Integrals</a></span></dt>
 <dt><span class="section"><a href="math_toolkit/naive_monte_carlo.html">Naive Monte Carlo Integration</a></span></dt>
+<dt><span class="section"><a href="math_toolkit/symplectic.html">Symplectic Integration</a></span></dt>
 <dt><span class="section"><a href="math_toolkit/wavelet_transforms.html">Wavelet Transforms</a></span></dt>
 <dt><span class="section"><a href="math_toolkit/diff.html">Numerical Differentiation</a></span></dt>
 <dt><span class="section"><a href="math_toolkit/autodiff.html">Automatic Differentiation</a></span></dt>
diff --git a/doc/math.qbk b/doc/math.qbk
index f48c87b88c..d2575db293 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -772,6 +772,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [include quadrature/double_exponential.qbk]
 [include quadrature/ooura_fourier_integrals.qbk]
 [include quadrature/naive_monte_carlo.qbk]
+[include quadrature/symplectic.qpk]
 [include quadrature/wavelet_transforms.qbk]
 [include differentiation/numerical_differentiation.qbk]
 [include differentiation/autodiff.qbk]
diff --git a/doc/quadrature/symplectic.qbk b/doc/quadrature/symplectic.qbk
new file mode 100644
index 0000000000..b52ee70467
--- /dev/null
+++ b/doc/quadrature/symplectic.qbk
@@ -0,0 +1,103 @@
+[/
+Copyright (c) 2026 Jacob Hass
+Use, modification and distribution are subject to the
+Boost Software License, Version 1.0. (See accompanying file
+LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+]
+
+[section:symplectic Symplectic Integration]
+
+[heading Synopsis]
+
+    #include <boost/math/quadrature/symmetric.hpp>
+    namespace boost{ namespace math{ namespace quadrature {
+
+    template <class ReturnType, class RealType, class Func>
+    std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
+                                                                                       const ReturnType q0,
+                                                                                       const RealType dt,
+                                                                                       const unsigned steps,
+                                                                                       Func dHdp,
+                                                                                       Func dHdq,
+                                                                                       std::string method = "Y6")
+
+    template <class ReturnType, class RealType, class Func, class ``__Policy``>
+    std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
+                                                                                       const ReturnType q0,
+                                                                                       const RealType dt,
+                                                                                       const unsigned steps,
+                                                                                       Func dHdp,
+                                                                                       Func dHdq,
+                                                                                       std::string method,
+                                                                                       const ``__Policy``& pol)
+
+    }}} // namespaces
+
+[heading Description]
+
+The functional `integrate_hamiltonian` calculates the phase space trajectory for a given Hamiltonian.
+The trajectories are calculated using symplectic integration which preserves the energy of 
+a system. Even higher order traditional numerical integration algorithms will gain or lose 
+energy at long times. The functional implements the methods in [@https://www.sciencedirect.com/science/article/abs/pii/0375960190900923 Yoshida's] 
+landmark paper. We assume that the Hamiltonian is separable so that it can be written as `H = T(p) + V(q)`. 
+
+We now give an example for a simple harmonic oscillator with the Hamiltonian
+
+[equation harmonic_oscillator]
+
+For simplicity we will assume `k = m = 1`. Then the partial derivatives of the Hamiltonian
+with respect to `p` and `q` are 
+
+    double dHdp(double p)
+    {
+        return p;
+    }
+
+    double dHdq(double q)
+    {
+        return q;
+    }
+
+Note that the functional can be readily generalized to multiple coordinates by changing the 
+signature of `dHdp` to 
+
+    std::valarray<double> dHdp(std::valarray<double> p)
+    {
+        // calculate the partial derivatives with respect to each p_i
+    }
+
+The function must return a `valarray` type, as opposed to `std::vector`, because we must be 
+able to perform arithmetic operations on the return type. We then define the timestep size 
+and number of steps of the algorithm to go from `t=0` to `t=100`
+
+    RealType dt = 0.05;
+    RealType t_end = 100;
+    unsigned int steps = t_end / dt;
+
+Lastly, we define the initial conditions so that the oscillator starts from rest at `x=1` 
+    
+    RealType q0 = 1;
+    RealType p0 = 0;
+
+We now evolve the system using the following 
+
+    std::vector<RealType> p;
+    std::vector<RealType> q;
+
+    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, "Y6");
+
+The ouput vectors `q, p` are the position and momentum of the system at each time. In 
+higher dimensions the output will be a vector of vectors. 
+
+The last argument that we pass to `integrate_hamiltonian`, "Y6" here, sets the integration 
+method to use. The string "Y6" stands for Yoshida's 6th order integrator. Yoshida's 2nd and 
+4th order methods are also available by passing the string "Y2", or "Y4" respectively.  
+
+[optional_policy]
+
+References:
+
+Yoshida, Haruo. [`Construction of higher order symplectic integrators], Physics Letters A, 150.5-7 (1990): 262-268.
+
+[endsect] [/section:symplectic Symplectic Integration]
+
diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index 1b270ed3c5..da5ae6547d 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -65,9 +65,15 @@ std::pair<ReturnType, ReturnType> sixth_order_yoshida(const ReturnType p0,
     ReturnType q = q0;
     
     // Choosing "System A" solution
-    RealType w3 = 0.784513610477560;
-    RealType w2 = 0.235573213359357;
-    RealType w1 = -1.17767998417887;
+    // The following Mathematica command can calculate these coefficients to arbitrary precision
+    // FindRoot[{w0+2(w1+w2+w3)==1,
+    //           w0^3 + 2(w1^3 + w2^3+w3^3)==0, 
+    //           w0^5 + 2(w1^5 + w2^5+w3^5)==0, 
+    //           1/6(w0^2w1^3-w0^4*w1) - 1/6(w0^3w1^2-w0*w1^4)+1/6((w0+2w1)^2w2^3-(w0+2w1)(w0^3+2w1^3)w2) - 1/6((w0^3+2w1^3)w2^2-(w0+2w1)w2^4) +1/6((w0+2w1+2w2)^2w3^3-(w0+2w1+2w2)(w0^3+2w1^3+2w2^3)w3)-1/6((w0^3+2w1^3+2w2^3)w3^2-(w0+2w1+2w2)w3^4)==0}, {{w0,1.3151863206839063}, {w1,-1.17767998417887}, {w2,0.235573213359357}, {w3,0.784513610477560}}, WorkingPrecision->64]
+    RealType w1 = static_cast<RealType>(-1.17767998417887100694641568096431573463926925263459848447536851379674155618156L);
+    RealType w2 = static_cast<RealType>(0.23557321335935813368479318297853460168646808210340111900349313095621471215223L);
+    RealType w3 = static_cast<RealType>(0.78451361047755726381949763386634987577682441745149338456794779895125997479548L);
+    // w0 = 1.31518632068391121888424972823886251435195350615940796180785516777853373846773
     RealType w0 = 1 - 2 * (w1 + w2 + w3);
     std::vector<RealType> weights = { w3, w2, w1, w0, w1, w2, w3};
 
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 7426c26fd6..694ccefdde 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -1341,7 +1341,7 @@ test-suite quadrature :
    [ run test_trapezoidal.cpp /boost/test//boost_unit_test_framework : : :
       release [ requires cxx11_lambdas cxx11_auto_declarations cxx11_decltype cxx11_unified_initialization_syntax cxx11_variadic_templates ]
       $(float128_type) ]
-   [ run test_symplectic.cpp /boost/test//boost_unit_test_framework ]
+   [ run test_symplectic.cpp ]
 ;
 
 test-suite autodiff :

From 12818850d14e1dc9531c795a2c5a14ca714c257a Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Wed, 24 Jun 2026 20:23:52 -0700
Subject: [PATCH 09/22] Fixed typo

---
 doc/math.qbk | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/math.qbk b/doc/math.qbk
index d2575db293..7e46490c61 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -772,7 +772,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [include quadrature/double_exponential.qbk]
 [include quadrature/ooura_fourier_integrals.qbk]
 [include quadrature/naive_monte_carlo.qbk]
-[include quadrature/symplectic.qpk]
+[include quadrature/symplectic.qbk]
 [include quadrature/wavelet_transforms.qbk]
 [include differentiation/numerical_differentiation.qbk]
 [include differentiation/autodiff.qbk]

From 5e7705b913aeac205b31ccfd900cacf68c1baf71 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Fri, 26 Jun 2026 15:44:03 -0700
Subject: [PATCH 10/22] Updated docs and max function

---
 doc/quadrature/symplectic.qbk                | 2 +-
 include/boost/math/quadrature/symplectic.hpp | 6 +++++-
 test/test_symplectic.cpp                     | 2 +-
 3 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/doc/quadrature/symplectic.qbk b/doc/quadrature/symplectic.qbk
index b52ee70467..1765a3482d 100644
--- a/doc/quadrature/symplectic.qbk
+++ b/doc/quadrature/symplectic.qbk
@@ -42,7 +42,7 @@ energy at long times. The functional implements the methods in [@https://www.sci
 landmark paper. We assume that the Hamiltonian is separable so that it can be written as `H = T(p) + V(q)`. 
 
 We now give an example for a simple harmonic oscillator with the Hamiltonian
-
+[/ $H = \frac{p^2}{2m} + \frac{1}{2}kx^2$ ]
 [equation harmonic_oscillator]
 
 For simplicity we will assume `k = m = 1`. Then the partial derivatives of the Hamiltonian
diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index da5ae6547d..f07486ab59 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -1,4 +1,8 @@
-/* Copyright Jacob Hass, 2026*/
+// Copyright Jacob Hass, 2026
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
 
 #ifndef BOOST_MATH_QUADRATURE_SYMPLECTIC_HPP
 #define BOOST_MATH_QUADRATURE_SYMPLECTIC_HPP
diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index f2608dd29c..efdf7f79dd 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -93,7 +93,7 @@ void test_harmonic_oscillator(RealType tol)
     std::valarray<RealType> q_val(q.data(), q.size());
     std::valarray<RealType> energy_error = std::pow(p_val, 2) + std::pow(q_val, 2) - 1;
     std::valarray<RealType> abs_energy_error = std::abs(energy_error);
-    RealType max_error = abs_energy_error.max();
+    RealType max_error = *std::max_element(std::begin(abs_energy_error), std::end(abs_energy_error));
     std::cout << max_error << std::endl;
     BOOST_CHECK_LE(max_error, tol);
 

From e31880b78252b5850d6c99600cd52cee980c9ace Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Fri, 26 Jun 2026 16:29:13 -0700
Subject: [PATCH 11/22] Changed functions to only accept vectors

---
 include/boost/math/quadrature/symplectic.hpp | 22 +++++--
 test/test_symplectic.cpp                     | 67 ++++++--------------
 2 files changed, 38 insertions(+), 51 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index f07486ab59..b71558b501 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -28,10 +28,24 @@ std::pair<ReturnType, ReturnType> second_order_yoshida(const ReturnType p0,
     ReturnType p = p0;
     ReturnType q = q0;
 
-    q = q + dt / 2 * dHdp(p);
-    p = p - dt * dHdq(q);
-    q = q + dt / 2 * dHdp(p);
-    
+    // Half step in q
+    auto dHdp_val = dHdp(p);
+    for (unsigned i=0; i < q.size(); i++){
+        q[i] = q[i] + dt / 2 * dHdp_val[i];
+    }
+
+    // Full step in p
+    auto dHdq_val = dHdq(q);
+    for (unsigned i=0; i < p.size(); i++){
+        p[i] = p[i] - dt * dHdq_val[i];
+    }
+
+    // Half step in q
+    dHdp_val = dHdp(p);
+    for (unsigned i=0; i < q.size(); i++){
+        q[i] = q[i] + dt / 2 * dHdp_val[i];
+    }
+
     return std::make_pair(p, q);
 }
 
diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index efdf7f79dd..3748579792 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -1,5 +1,5 @@
 /*
- * Copyright Nick Thompson, 2017
+ * Copyright Jacob Hass, 2026
  * Use, modification and distribution are subject to the
  * Boost Software License, Version 1.0. (See accompanying file
  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
@@ -24,53 +24,22 @@
 
 using boost::math::quadrature::integrate_hamiltonian;
 
-// Equations of motion for simple harmonic oscillator 
+// Equations of motion for simple harmonic oscillator
 template <class Real>
-Real singleton_dHdp(Real p) { return p; }
+std::vector<Real> vector_dHdp(std::vector<Real> p) { return p; }
 template <class Real>
-Real singleton_dHdq(Real q) { return q; }
-
-template <class Real>
-std::valarray<Real> vector_dHdp(std::valarray<Real> p) { return p; }
-template <class Real>
-std::valarray<Real> vector_dHdq(std::valarray<Real> q) { return q; }
-
-// Sanity check that passing as a 1D vector does not effect results
-template<class RealType>
-void test_input_singleton_vs_vector()
-{
-    RealType dt = 0.05;
-    unsigned int steps = 3;
-    RealType q0 = 1;
-    RealType p0 = 0;
-
-
-    std::vector<RealType> p;
-    std::vector<RealType> q;
-
-    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, singleton_dHdp<RealType>, singleton_dHdq<RealType>);
-
-    std::valarray<RealType> q0_vector = {1};
-    std::valarray<RealType> p0_vector = {0};
-    std::vector<std::valarray<RealType> > p_vector;
-    std::vector<std::valarray<RealType> > q_vector;
-
-    std::tie(p_vector, q_vector) = boost::math::quadrature::integrate_hamiltonian(p0_vector, q0_vector, dt, steps, vector_dHdp<RealType>, vector_dHdq<RealType>);
-    
-    BOOST_CHECK_EQUAL(p[0], p_vector[0][0]);
-    BOOST_CHECK_EQUAL(p[1], p_vector[1][0]);
-    BOOST_CHECK_EQUAL(p[2], p_vector[2][0]);
-
-}
+std::vector<Real> vector_dHdq(std::vector<Real> q) { return q; }
 
 template <class RealType>
 void test_invalid_parameters()
 {
+    std::vector<RealType> q0 = {1};
+    std::vector<RealType> p0 = {0};
     // Negative timestep
-    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(1, 0, -0.1, 10, singleton_dHdp<RealType>, singleton_dHdq<RealType>), std::domain_error);
+    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(q0, p0, -0.1, 10, vector_dHdp<RealType>, vector_dHdq<RealType>), std::domain_error);
 
     // Method not in {'Y6', 'Y4', 'Y2'}
-    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(1, 0, 0.1, 10, singleton_dHdp<RealType>, singleton_dHdq<RealType>, "InvalidMethod"), std::out_of_range);
+    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(q0, p0, 0.1, 10, vector_dHdp<RealType>, vector_dHdq<RealType>, "InvalidMethod"), std::out_of_range);
 }
 
 /* Test if SHO energy fluctuations are below a given tolerance*/
@@ -81,16 +50,21 @@ void test_harmonic_oscillator(RealType tol)
     RealType t_end = 100;
     unsigned int steps = t_end / dt;
 
-    RealType q0 = 1;
-    RealType p0 = 0;
+    std::vector<RealType> q0 = {1};
+    std::vector<RealType> p0 = {0};
 
-    std::vector<RealType> p;
-    std::vector<RealType> q;
+    std::vector<std::vector<RealType> > p;
+    std::vector<std::vector<RealType> > q;
 
-    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, singleton_dHdp<RealType>, singleton_dHdq<RealType>);
+    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, vector_dHdp<RealType>, vector_dHdp<RealType>);
 
-    std::valarray<RealType> p_val(p.data(), p.size());
-    std::valarray<RealType> q_val(q.data(), q.size());
+    std::valarray<RealType> p_val(p.size());
+    std::valarray<RealType> q_val(q.size());
+    for (unsigned i=0; i < p.size(); i++)
+    {
+        p_val[i] = p[i][0];
+        q_val[i] = q[i][0];
+    }
     std::valarray<RealType> energy_error = std::pow(p_val, 2) + std::pow(q_val, 2) - 1;
     std::valarray<RealType> abs_energy_error = std::abs(energy_error);
     RealType max_error = *std::max_element(std::begin(abs_energy_error), std::end(abs_energy_error));
@@ -101,7 +75,6 @@ void test_harmonic_oscillator(RealType tol)
 
 BOOST_AUTO_TEST_CASE(symplectic_quadrature)
 {
-    test_input_singleton_vs_vector<double>();
     test_invalid_parameters<double>();
     test_harmonic_oscillator<double>(1e-10);
 }

From 645aa1a8981c80b2173e5a5a1924815bd3d5771b Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Fri, 26 Jun 2026 19:21:07 -0700
Subject: [PATCH 12/22] Completely eliminated valarrays

---
 test/test_symplectic.cpp | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index 3748579792..a1088ef2b9 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -7,7 +7,6 @@
 #define BOOST_TEST_MODULE symplectic_quadrature
 
 #include <algorithm>
-#include <valarray>
 #include <boost/math/tools/config.hpp>
 #include <boost/test/included/unit_test.hpp>
 #include <boost/test/tools/floating_point_comparison.hpp>
@@ -46,6 +45,8 @@ void test_invalid_parameters()
 template <class RealType>
 void test_harmonic_oscillator(RealType tol)
 {
+    BOOST_MATH_STD_USING;
+
     RealType dt = 0.05;
     RealType t_end = 100;
     unsigned int steps = t_end / dt;
@@ -58,15 +59,17 @@ void test_harmonic_oscillator(RealType tol)
 
     std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, vector_dHdp<RealType>, vector_dHdp<RealType>);
 
-    std::valarray<RealType> p_val(p.size());
-    std::valarray<RealType> q_val(q.size());
+    RealType p_val;
+    RealType q_val;
+    std::vector<RealType> abs_energy_error(p.size());
     for (unsigned i=0; i < p.size(); i++)
     {
-        p_val[i] = p[i][0];
-        q_val[i] = q[i][0];
+        p_val = p[i][0];
+        q_val = q[i][0];
+
+        abs_energy_error[i] = std::abs(std::pow(p_val, 2) + std::pow(q_val, 2) - 1);
     }
-    std::valarray<RealType> energy_error = std::pow(p_val, 2) + std::pow(q_val, 2) - 1;
-    std::valarray<RealType> abs_energy_error = std::abs(energy_error);
+
     RealType max_error = *std::max_element(std::begin(abs_energy_error), std::end(abs_energy_error));
     std::cout << max_error << std::endl;
     BOOST_CHECK_LE(max_error, tol);

From cdd64cf715877afa6c2495df7117b0837c4ee92d Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 5 Jul 2026 16:41:48 -0700
Subject: [PATCH 13/22] Added a couple more tests

---
 include/boost/math/quadrature/symplectic.hpp |  5 ++++-
 test/test_symplectic.cpp                     | 12 ++++++------
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index b71558b501..d2c6e3e788 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -56,11 +56,14 @@ std::pair<ReturnType, ReturnType> fourth_order_yoshida(const ReturnType p0,
                                                        Func dHdp, 
                                                        Func dHdq)
 {
+    BOOST_MATH_STD_USING
+    
     ReturnType p = p0;
     ReturnType q = q0;
 
-    RealType x0 = -(std::pow(2, 1/3) / (2 - std::pow(2, 1/3)));
+    // RealType x0 = -(std::pow(2, 1/3) / (2 - std::pow(2, 1/3)));
     RealType x1 = 1 / (2 - std::pow(2, 1/3));
+    RealType x0 = 1 - 2 * x1; 
 
     std::vector<RealType> weights = { x1, x0, x1 };
 
diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index a1088ef2b9..853597fea3 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -43,9 +43,9 @@ void test_invalid_parameters()
 
 /* Test if SHO energy fluctuations are below a given tolerance*/
 template <class RealType>
-void test_harmonic_oscillator(RealType tol)
+void test_harmonic_oscillator(const RealType tol, const std::string method)
 {
-    BOOST_MATH_STD_USING;
+    BOOST_MATH_STD_USING
 
     RealType dt = 0.05;
     RealType t_end = 100;
@@ -57,7 +57,7 @@ void test_harmonic_oscillator(RealType tol)
     std::vector<std::vector<RealType> > p;
     std::vector<std::vector<RealType> > q;
 
-    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, vector_dHdp<RealType>, vector_dHdp<RealType>);
+    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, vector_dHdp<RealType>, vector_dHdp<RealType>, method);
 
     RealType p_val;
     RealType q_val;
@@ -71,13 +71,13 @@ void test_harmonic_oscillator(RealType tol)
     }
 
     RealType max_error = *std::max_element(std::begin(abs_energy_error), std::end(abs_energy_error));
-    std::cout << max_error << std::endl;
     BOOST_CHECK_LE(max_error, tol);
-
 }
 
 BOOST_AUTO_TEST_CASE(symplectic_quadrature)
 {
     test_invalid_parameters<double>();
-    test_harmonic_oscillator<double>(1e-10);
+    test_harmonic_oscillator<double>(1e-10, "Y6");
+    test_harmonic_oscillator<double>(7e-4, "Y4");
+    test_harmonic_oscillator<double>(7e-4, "Y2");
 }

From 290e3266034820c389e7654d852236c6f5552d64 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 5 Jul 2026 16:58:30 -0700
Subject: [PATCH 14/22] Changed function namespaces

---
 include/boost/math/quadrature/symplectic.hpp | 32 ++++++++++++++------
 1 file changed, 23 insertions(+), 9 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index d2c6e3e788..45168d78c5 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -16,7 +16,7 @@
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
 
-namespace boost{ namespace math {namespace quadrature { 
+namespace boost{ namespace math { namespace quadrature { namespace detail {
 
 template <class RealType, class ReturnType, class Func>
 std::pair<ReturnType, ReturnType> second_order_yoshida(const ReturnType p0, 
@@ -107,14 +107,14 @@ std::pair<ReturnType, ReturnType> sixth_order_yoshida(const ReturnType p0,
 }
 
 template <class ReturnType, class RealType, class Func, class Policy>
-std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
-                                                                                   const ReturnType q0,
-                                                                                   const RealType dt,
-                                                                                   const unsigned steps,
-                                                                                   Func dHdp,
-                                                                                   Func dHdq,
-                                                                                   std::string method,
-                                                                                   const Policy& pol)
+std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian_imp(const ReturnType p0,
+                                                                                       const ReturnType q0,
+                                                                                       const RealType dt,
+                                                                                       const unsigned steps,
+                                                                                       Func dHdp,
+                                                                                       Func dHdq,
+                                                                                       std::string method,
+                                                                                       const Policy& pol)
 {
     // Not sure how to make this function string nicer
     static const char* function = "boost::math::quadrature::integrate_hamiltonian(p0, q0, %1%, steps, dHdp, dHdq)";
@@ -151,6 +151,20 @@ std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltoni
     }
     return std::make_pair(p, q);
 }
+} // namespace detail
+
+template <class ReturnType, class RealType, class Func, class Policy>
+std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
+                                                                                   const ReturnType q0,
+                                                                                   const RealType dt,
+                                                                                   const unsigned steps,
+                                                                                   Func dHdp,
+                                                                                   Func dHdq,
+                                                                                   std::string method,
+                                                                                   const Policy& pol)
+{
+    return detail::integrate_hamiltonian_imp(p0, q0, dt, steps, dHdp, dHdq, method, pol); 
+}
 
 template <class ReturnType, class RealType, class Func>
 std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,

From bb09041a0a94c2993020bc1f925ed45cce88779c Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 5 Jul 2026 18:51:04 -0700
Subject: [PATCH 15/22] Added SRKN from paper

---
 include/boost/math/quadrature/symplectic.hpp | 53 +++++++++++++++++++-
 test/test_symplectic.cpp                     |  1 +
 2 files changed, 52 insertions(+), 2 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index 45168d78c5..ba7bec42b3 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -95,7 +95,7 @@ std::pair<ReturnType, ReturnType> sixth_order_yoshida(const ReturnType p0,
     RealType w2 = static_cast<RealType>(0.23557321335935813368479318297853460168646808210340111900349313095621471215223L);
     RealType w3 = static_cast<RealType>(0.78451361047755726381949763386634987577682441745149338456794779895125997479548L);
     // w0 = 1.31518632068391121888424972823886251435195350615940796180785516777853373846773
-    RealType w0 = 1 - 2 * (w1 + w2 + w3);
+    RealType w0 = 1. - 2. * (w1 + w2 + w3);
     std::vector<RealType> weights = { w3, w2, w1, w0, w1, w2, w3};
 
     for (unsigned i=0; i < weights.size(); i++)
@@ -106,6 +106,54 @@ std::pair<ReturnType, ReturnType> sixth_order_yoshida(const ReturnType p0,
     return std::make_pair(p, q);
 }
 
+template <class RealType, class ReturnType, class Func>
+std::pair<ReturnType, ReturnType> SRKN_b_order_6(const ReturnType p0,
+                                                 const ReturnType q0,
+                                                 RealType dt, 
+                                                 Func dHdp,
+                                                 Func dHdq)
+{ // This method implements SRKN_b^6 in Table 3 here 
+  // https://www.sciencedirect.com/science/article/pii/S0377042701004927
+
+    ReturnType p = p0;
+    ReturnType q = q0;
+    
+    RealType b1 = static_cast<RealType>(0.0829844064174052);
+    RealType b2 = static_cast<RealType>(0.396309801498368);
+    RealType b3 = static_cast<RealType>(-0.0390563049223486);
+    RealType b4 = 1. - 2. * (b1 + b2 + b3);
+
+    RealType a1 = static_cast<RealType>(0.245298957184271);
+    RealType a2 = static_cast<RealType>(0.604872665711080);
+    RealType a3 = 0.5 - (a1 + a2);
+
+    std::vector<RealType> b_weights = {b1, b2, b3, b4, b3, b2};
+    std::vector<RealType> a_weights = {a1, a2, a3, a3, a2, a1};
+    
+    RealType a, b;
+    for (unsigned int i=0; i < b_weights.size(); i++)
+    {
+        b = b_weights[i];
+        a = a_weights[i];
+        
+        auto dHdp_val = dHdp(p);
+        for (unsigned j=0; j < q.size(); j++){
+            q[j] = q[j] + b * dt * dHdp_val[j];
+        }
+        
+        auto dHdq_val = dHdq(q);
+        for (unsigned j=0; j < p.size(); j++){
+            p[j] = p[j] - a * dt * dHdq_val[j];
+        }
+    }
+    // Have to do one more step in p
+    auto dHdp_val = dHdp(p);
+    for (unsigned j=0; j < q.size(); j++){
+        q[j] = q[j] + b1 * dt * dHdp_val[j];
+    }
+    return std::make_pair(p, q);
+}
+
 template <class ReturnType, class RealType, class Func, class Policy>
 std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian_imp(const ReturnType p0,
                                                                                        const ReturnType q0,
@@ -133,7 +181,8 @@ std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltoni
 
     std::map<std::string, stepperType> m{{"Y6", sixth_order_yoshida}, 
                                          {"Y4", fourth_order_yoshida}, 
-                                         {"Y2", second_order_yoshida}};
+                                         {"Y2", second_order_yoshida},
+                                         {"SRKNB6", SRKN_b_order_6}};
     stepperType stepper = m.at(method);
 
     std::vector<ReturnType> p(steps);
diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index 853597fea3..6dc60098a9 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -80,4 +80,5 @@ BOOST_AUTO_TEST_CASE(symplectic_quadrature)
     test_harmonic_oscillator<double>(1e-10, "Y6");
     test_harmonic_oscillator<double>(7e-4, "Y4");
     test_harmonic_oscillator<double>(7e-4, "Y2");
+    test_harmonic_oscillator<double>(1e-11, "SRKNB6");
 }

From e6d67b015c1df2c9d63d7181404782315b07a352 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 5 Jul 2026 20:22:24 -0700
Subject: [PATCH 16/22] Added simple pendulum test

---
 test/test_symplectic.cpp | 60 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 59 insertions(+), 1 deletion(-)

diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index 6dc60098a9..e4124e4aaf 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -13,6 +13,7 @@
 #include <boost/math/tools/test_value.hpp>
 #include <boost/math/concepts/real_concept.hpp>
 #include <boost/math/quadrature/symplectic.hpp>
+#include <boost/math/constants/constants.hpp>
 #ifdef BOOST_HAS_FLOAT128
 #include <boost/multiprecision/complex128.hpp>
 #endif
@@ -29,6 +30,21 @@ std::vector<Real> vector_dHdp(std::vector<Real> p) { return p; }
 template <class Real>
 std::vector<Real> vector_dHdq(std::vector<Real> q) { return q; }
 
+// Equations of motion for simple pendulum
+template <class Real>
+std::vector<Real> pendulum_vector_dHdp(std::vector<Real> p) { return p; }
+template <class Real>
+std::vector<Real> pendulum_vector_dHdq(std::vector<Real> q) 
+{ 
+    BOOST_MATH_STD_USING
+    std::vector<Real> q_return(q.size());
+    for (unsigned i=0; i<q.size(); i++)
+    {
+        q_return[i] = std::sin(q[i]);
+    }
+    return q_return;
+}
+
 template <class RealType>
 void test_invalid_parameters()
 {
@@ -57,7 +73,7 @@ void test_harmonic_oscillator(const RealType tol, const std::string method)
     std::vector<std::vector<RealType> > p;
     std::vector<std::vector<RealType> > q;
 
-    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, vector_dHdp<RealType>, vector_dHdp<RealType>, method);
+    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, vector_dHdp<RealType>, vector_dHdq<RealType>, method);
 
     RealType p_val;
     RealType q_val;
@@ -74,11 +90,53 @@ void test_harmonic_oscillator(const RealType tol, const std::string method)
     BOOST_CHECK_LE(max_error, tol);
 }
 
+/* Test if SHO energy fluctuations are below a given tolerance*/
+template <class RealType>
+void test_pendulum(const RealType tol, const std::string method)
+{
+    BOOST_MATH_STD_USING
+
+    RealType dt = 0.05;
+    RealType t_end = 100;
+    unsigned int steps = t_end / dt;
+
+    std::vector<RealType> q0 = {boost::math::constants::pi<RealType>() / 2};
+    std::vector<RealType> p0 = {0};
+
+    std::vector<std::vector<RealType> > p;
+    std::vector<std::vector<RealType> > q;
+
+    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, pendulum_vector_dHdp<RealType>, pendulum_vector_dHdq<RealType>, method);
+
+    RealType p_val;
+    RealType q_val;
+    std::vector<RealType> abs_energy_error(p.size());
+    for (unsigned i=0; i < p.size(); i++)
+    {
+        p_val = p[i][0];
+        q_val = q[i][0];
+        // The energy of the system is H = p^2 / 2 + (1-cos(q)). With our initial condition
+        // this yields H = 1. Thus, we just want p^2 / 2 - cos(q) = 0
+        abs_energy_error[i] = std::abs(std::pow(p_val, 2) / 2 - std::cos(q_val));
+    }
+
+    RealType max_error = *std::max_element(std::begin(abs_energy_error), std::end(abs_energy_error));
+    BOOST_CHECK_LE(max_error, tol);
+}
+
 BOOST_AUTO_TEST_CASE(symplectic_quadrature)
 {
     test_invalid_parameters<double>();
+    
+    // SHO Tests
     test_harmonic_oscillator<double>(1e-10, "Y6");
     test_harmonic_oscillator<double>(7e-4, "Y4");
     test_harmonic_oscillator<double>(7e-4, "Y2");
     test_harmonic_oscillator<double>(1e-11, "SRKNB6");
+
+    // Pendulum Tests
+    test_pendulum<double>(1e-10, "Y6");
+    test_pendulum<double>(5e-4, "Y4");
+    test_pendulum<double>(5e-4, "Y2");
+    test_pendulum<double>(1e-8, "SRKNB6");
 }

From bc972a6ff358ed42c4e86a7653cad7a0b3c8039d Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Mon, 6 Jul 2026 20:49:31 -0700
Subject: [PATCH 17/22] Changed template to RandomAccessContainer

---
 include/boost/math/quadrature/symplectic.hpp | 150 ++++++++++---------
 1 file changed, 78 insertions(+), 72 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index ba7bec42b3..3e584564d6 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -18,15 +18,15 @@
 
 namespace boost{ namespace math { namespace quadrature { namespace detail {
 
-template <class RealType, class ReturnType, class Func>
-std::pair<ReturnType, ReturnType> second_order_yoshida(const ReturnType p0, 
-                                                       const ReturnType q0, 
-                                                       RealType dt, 
-                                                       Func dHdp, 
-                                                       Func dHdq)
+template <typename RandomAccessContainer, class Func>
+std::pair<RandomAccessContainer, RandomAccessContainer> second_order_yoshida(const RandomAccessContainer p0, 
+                                                                             const RandomAccessContainer q0, 
+                                                                             const typename RandomAccessContainer::value_type dt, 
+                                                                             Func dHdp, 
+                                                                             Func dHdq)
 {
-    ReturnType p = p0;
-    ReturnType q = q0;
+    RandomAccessContainer p = p0;
+    RandomAccessContainer q = q0;
 
     // Half step in q
     auto dHdp_val = dHdp(p);
@@ -49,17 +49,18 @@ std::pair<ReturnType, ReturnType> second_order_yoshida(const ReturnType p0,
     return std::make_pair(p, q);
 }
 
-template <class RealType, class ReturnType, class Func>
-std::pair<ReturnType, ReturnType> fourth_order_yoshida(const ReturnType p0, 
-                                                       const ReturnType q0, 
-                                                       const RealType dt, 
-                                                       Func dHdp, 
-                                                       Func dHdq)
+template <typename RandomAccessContainer, class Func>
+std::pair<RandomAccessContainer, RandomAccessContainer> fourth_order_yoshida(const RandomAccessContainer p0, 
+                                                                             const RandomAccessContainer q0, 
+                                                                             const typename RandomAccessContainer::value_type dt, 
+                                                                             Func dHdp, 
+                                                                             Func dHdq)
 {
     BOOST_MATH_STD_USING
-    
-    ReturnType p = p0;
-    ReturnType q = q0;
+    using RealType = typename RandomAccessContainer::value_type;
+
+    RandomAccessContainer p = p0;
+    RandomAccessContainer q = q0;
 
     // RealType x0 = -(std::pow(2, 1/3) / (2 - std::pow(2, 1/3)));
     RealType x1 = 1 / (2 - std::pow(2, 1/3));
@@ -75,15 +76,17 @@ std::pair<ReturnType, ReturnType> fourth_order_yoshida(const ReturnType p0,
     return std::make_pair(p, q);
 }
 
-template <class RealType, class ReturnType, class Func>
-std::pair<ReturnType, ReturnType> sixth_order_yoshida(const ReturnType p0, 
-                                                      const ReturnType q0, 
-                                                      RealType dt, 
-                                                      Func dHdp, 
-                                                      Func dHdq)
+template <typename RandomAccessContainer, class Func>
+std::pair<RandomAccessContainer, RandomAccessContainer> sixth_order_yoshida(const RandomAccessContainer p0, 
+                                                                            const RandomAccessContainer q0, 
+                                                                            const typename RandomAccessContainer::value_type dt, 
+                                                                            Func dHdp, 
+                                                                            Func dHdq)
 {
-    ReturnType p = p0;
-    ReturnType q = q0;
+    using RealType = typename RandomAccessContainer::value_type;
+
+    RandomAccessContainer p = p0;
+    RandomAccessContainer q = q0;
     
     // Choosing "System A" solution
     // The following Mathematica command can calculate these coefficients to arbitrary precision
@@ -106,17 +109,19 @@ std::pair<ReturnType, ReturnType> sixth_order_yoshida(const ReturnType p0,
     return std::make_pair(p, q);
 }
 
-template <class RealType, class ReturnType, class Func>
-std::pair<ReturnType, ReturnType> SRKN_b_order_6(const ReturnType p0,
-                                                 const ReturnType q0,
-                                                 RealType dt, 
-                                                 Func dHdp,
-                                                 Func dHdq)
+template <typename RandomAccessContainer, class Func>
+std::pair<RandomAccessContainer, RandomAccessContainer> SRKN_b_order_6(const RandomAccessContainer p0,
+                                                                       const RandomAccessContainer q0,
+                                                                       const typename RandomAccessContainer::value_type dt, 
+                                                                       Func dHdp,
+                                                                       Func dHdq)
 { // This method implements SRKN_b^6 in Table 3 here 
   // https://www.sciencedirect.com/science/article/pii/S0377042701004927
 
-    ReturnType p = p0;
-    ReturnType q = q0;
+    using RealType = typename RandomAccessContainer::value_type;
+
+    RandomAccessContainer p = p0;
+    RandomAccessContainer q = q0;
     
     RealType b1 = static_cast<RealType>(0.0829844064174052);
     RealType b2 = static_cast<RealType>(0.396309801498368);
@@ -154,30 +159,31 @@ std::pair<ReturnType, ReturnType> SRKN_b_order_6(const ReturnType p0,
     return std::make_pair(p, q);
 }
 
-template <class ReturnType, class RealType, class Func, class Policy>
-std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian_imp(const ReturnType p0,
-                                                                                       const ReturnType q0,
-                                                                                       const RealType dt,
-                                                                                       const unsigned steps,
-                                                                                       Func dHdp,
-                                                                                       Func dHdq,
-                                                                                       std::string method,
-                                                                                       const Policy& pol)
+template <typename RandomAccessContainer, class Func, class Policy>
+std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian_imp(const RandomAccessContainer p0,
+                                                                                                             const RandomAccessContainer q0,
+                                                                                                             const typename RandomAccessContainer::value_type dt,
+                                                                                                             const unsigned steps,
+                                                                                                             Func dHdp,
+                                                                                                             Func dHdq,
+                                                                                                             std::string method,
+                                                                                                             const Policy& pol)
 {
+    using RealType = typename RandomAccessContainer::value_type;
     // Not sure how to make this function string nicer
     static const char* function = "boost::math::quadrature::integrate_hamiltonian(p0, q0, %1%, steps, dHdp, dHdq)";
 
     if ((dt <= 0) || !(boost::math::isfinite)(dt))
     {
         RealType val = (boost::math::policies::raise_domain_error(function, "Time step must be positive and finite but got: dt = %1%.\n", dt, pol));
-        std::vector<ReturnType> nan_vec = {ReturnType(val)};
+        std::vector<RandomAccessContainer> nan_vec = {RandomAccessContainer(val)};
         return std::make_pair(nan_vec, nan_vec);
     }
 
     // Check if method is available
     std::vector<std::string> available_methods = {"Y6", "Y4", "Y2"};
 
-    typedef std::pair<ReturnType, ReturnType> (*stepperType)(ReturnType, ReturnType, RealType, Func, Func);
+    typedef std::pair<RandomAccessContainer, RandomAccessContainer> (*stepperType)(RandomAccessContainer, RandomAccessContainer, RealType, Func, Func);
 
     std::map<std::string, stepperType> m{{"Y6", sixth_order_yoshida}, 
                                          {"Y4", fourth_order_yoshida}, 
@@ -185,13 +191,13 @@ std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltoni
                                          {"SRKNB6", SRKN_b_order_6}};
     stepperType stepper = m.at(method);
 
-    std::vector<ReturnType> p(steps);
-    std::vector<ReturnType> q(steps);
+    std::vector<RandomAccessContainer> p(steps);
+    std::vector<RandomAccessContainer> q(steps);
     p[0] = p0;
     q[0] = q0;
 
-    ReturnType p_current = p0;
-    ReturnType q_current = q0;
+    RandomAccessContainer p_current = p0;
+    RandomAccessContainer q_current = q0;
     for (unsigned i=1; i < steps; i++)
     {
         std::tie(p_current, q_current) = stepper(p_current, q_current, dt, dHdp, dHdq);
@@ -202,38 +208,38 @@ std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltoni
 }
 } // namespace detail
 
-template <class ReturnType, class RealType, class Func, class Policy>
-std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
-                                                                                   const ReturnType q0,
-                                                                                   const RealType dt,
-                                                                                   const unsigned steps,
-                                                                                   Func dHdp,
-                                                                                   Func dHdq,
-                                                                                   std::string method,
-                                                                                   const Policy& pol)
+template <typename RandomAccessContainer, class Func, class Policy>
+std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian(const RandomAccessContainer p0,
+                                                                                                         const RandomAccessContainer q0,
+                                                                                                         const typename RandomAccessContainer::value_type dt,
+                                                                                                         const unsigned steps,
+                                                                                                         Func dHdp,
+                                                                                                         Func dHdq,
+                                                                                                         std::string method,
+                                                                                                         const Policy& pol)
 {
     return detail::integrate_hamiltonian_imp(p0, q0, dt, steps, dHdp, dHdq, method, pol); 
 }
 
-template <class ReturnType, class RealType, class Func>
-std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
-                                                                                   const ReturnType q0,
-                                                                                   const RealType dt,
-                                                                                   const unsigned steps,
-                                                                                   Func dHdp,
-                                                                                   Func dHdq,
-                                                                                   std::string method)
+template <typename RandomAccessContainer, class Func>
+std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian(const RandomAccessContainer p0,
+                                                                                                         const RandomAccessContainer q0,
+                                                                                                         const typename RandomAccessContainer::value_type dt,
+                                                                                                         const unsigned steps,
+                                                                                                         Func dHdp,
+                                                                                                         Func dHdq,
+                                                                                                         std::string method)
 {
     return integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, method, boost::math::policies::policy<>()); 
 }
 
-template <class ReturnType, class RealType, class Func>
-std::pair<std::vector<ReturnType>, std::vector<ReturnType> > integrate_hamiltonian(const ReturnType p0,
-                                                                                   const ReturnType q0,
-                                                                                   const RealType dt,
-                                                                                   const unsigned steps,
-                                                                                   Func dHdp,
-                                                                                   Func dHdq)
+template <typename RandomAccessContainer, class Func>
+std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian(const RandomAccessContainer p0,
+                                                                                                         const RandomAccessContainer q0,
+                                                                                                         const typename RandomAccessContainer::value_type dt,
+                                                                                                         const unsigned steps,
+                                                                                                         Func dHdp,
+                                                                                                         Func dHdq)
 {
     return integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, "Y6", boost::math::policies::policy<>()); 
 }

From 40eea88c566fcde0705382e82527bcf4c8b250ea Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sat, 11 Jul 2026 19:58:55 -0700
Subject: [PATCH 18/22] Added singleton support

---
 include/boost/math/quadrature/symplectic.hpp | 137 ++++++++++++-------
 test/test_symplectic.cpp                     |  28 ++--
 2 files changed, 104 insertions(+), 61 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index 3e584564d6..310d23c46a 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -18,10 +18,64 @@
 
 namespace boost{ namespace math { namespace quadrature { namespace detail {
 
-template <typename RandomAccessContainer, class Func>
+template<typename...>
+using void_t = void;
+
+template<typename T, typename U, typename = void>
+struct has_plus : std::false_type {};
+
+template<typename T, typename U>
+struct has_plus<T, U, std::void_t<decltype(std::declval<T>() + std::declval<U>())> > : std::true_type {};
+
+template <typename T, typename U>
+typename std::enable_if<has_plus<T, U>::value, T>::type
+add(T x, U y) 
+{
+    return x + y;
+}
+
+template <typename T, typename U>
+typename std::enable_if<!has_plus<T, U>::value, T>::type
+add(T vec1, U vec2) 
+{
+    for (unsigned i=0; i < vec1.size(); i++)
+    {
+        vec1[i] = vec1[i] + vec2[i];
+    }
+    return vec1;
+}
+
+template<typename...>
+using void_t = void;
+
+template<typename T, typename U, typename = void>
+struct has_mult : std::false_type {};
+
+template<typename T, typename U>
+struct has_mult<T, U, std::void_t<decltype(std::declval<T>() * std::declval<U>())> > : std::true_type {};
+
+template <typename T, typename U>
+typename std::enable_if<has_mult<T, U>::value, T>::type
+mult_prefactor(T x, U prefactor) 
+{
+    return x * prefactor;
+}
+
+template <typename T, typename U>
+typename std::enable_if<!has_plus<T, U>::value, T>::type
+mult_prefactor(T vec1, U prefactor) 
+{
+    for (unsigned i=0; i < vec1.size(); i++)
+    {
+        vec1[i] = vec1[i] * prefactor;
+    }
+    return vec1;
+}
+
+template <typename RandomAccessContainer, typename RealType, class Func>
 std::pair<RandomAccessContainer, RandomAccessContainer> second_order_yoshida(const RandomAccessContainer p0, 
                                                                              const RandomAccessContainer q0, 
-                                                                             const typename RandomAccessContainer::value_type dt, 
+                                                                             const RealType dt, 
                                                                              Func dHdp, 
                                                                              Func dHdq)
 {
@@ -29,35 +83,31 @@ std::pair<RandomAccessContainer, RandomAccessContainer> second_order_yoshida(con
     RandomAccessContainer q = q0;
 
     // Half step in q
-    auto dHdp_val = dHdp(p);
-    for (unsigned i=0; i < q.size(); i++){
-        q[i] = q[i] + dt / 2 * dHdp_val[i];
-    }
+    RandomAccessContainer dHdp_val = dHdp(p);
+    RandomAccessContainer dq = mult_prefactor(dHdp_val, dt / 2);
+    q = add(q, dq);
 
     // Full step in p
-    auto dHdq_val = dHdq(q);
-    for (unsigned i=0; i < p.size(); i++){
-        p[i] = p[i] - dt * dHdq_val[i];
-    }
+    RandomAccessContainer dHdq_val = dHdq(q);
+    RandomAccessContainer dp = mult_prefactor(dHdq_val, -dt);
+    p = add(p, dp);
 
     // Half step in q
     dHdp_val = dHdp(p);
-    for (unsigned i=0; i < q.size(); i++){
-        q[i] = q[i] + dt / 2 * dHdp_val[i];
-    }
+    dq = mult_prefactor(dHdp_val, dt / 2);
+    q = add(q, dq);
 
     return std::make_pair(p, q);
 }
 
-template <typename RandomAccessContainer, class Func>
+template <typename RandomAccessContainer, typename RealType, class Func>
 std::pair<RandomAccessContainer, RandomAccessContainer> fourth_order_yoshida(const RandomAccessContainer p0, 
                                                                              const RandomAccessContainer q0, 
-                                                                             const typename RandomAccessContainer::value_type dt, 
+                                                                             const RealType dt, 
                                                                              Func dHdp, 
                                                                              Func dHdq)
 {
     BOOST_MATH_STD_USING
-    using RealType = typename RandomAccessContainer::value_type;
 
     RandomAccessContainer p = p0;
     RandomAccessContainer q = q0;
@@ -76,15 +126,13 @@ std::pair<RandomAccessContainer, RandomAccessContainer> fourth_order_yoshida(con
     return std::make_pair(p, q);
 }
 
-template <typename RandomAccessContainer, class Func>
+template <typename RandomAccessContainer, typename RealType, class Func>
 std::pair<RandomAccessContainer, RandomAccessContainer> sixth_order_yoshida(const RandomAccessContainer p0, 
                                                                             const RandomAccessContainer q0, 
-                                                                            const typename RandomAccessContainer::value_type dt, 
+                                                                            const RealType dt, 
                                                                             Func dHdp, 
                                                                             Func dHdq)
 {
-    using RealType = typename RandomAccessContainer::value_type;
-
     RandomAccessContainer p = p0;
     RandomAccessContainer q = q0;
     
@@ -109,17 +157,15 @@ std::pair<RandomAccessContainer, RandomAccessContainer> sixth_order_yoshida(cons
     return std::make_pair(p, q);
 }
 
-template <typename RandomAccessContainer, class Func>
+template <typename RandomAccessContainer, typename RealType, class Func>
 std::pair<RandomAccessContainer, RandomAccessContainer> SRKN_b_order_6(const RandomAccessContainer p0,
                                                                        const RandomAccessContainer q0,
-                                                                       const typename RandomAccessContainer::value_type dt, 
+                                                                       const RealType dt, 
                                                                        Func dHdp,
                                                                        Func dHdq)
 { // This method implements SRKN_b^6 in Table 3 here 
   // https://www.sciencedirect.com/science/article/pii/S0377042701004927
 
-    using RealType = typename RandomAccessContainer::value_type;
-
     RandomAccessContainer p = p0;
     RandomAccessContainer q = q0;
     
@@ -141,35 +187,32 @@ std::pair<RandomAccessContainer, RandomAccessContainer> SRKN_b_order_6(const Ran
         b = b_weights[i];
         a = a_weights[i];
         
-        auto dHdp_val = dHdp(p);
-        for (unsigned j=0; j < q.size(); j++){
-            q[j] = q[j] + b * dt * dHdp_val[j];
-        }
-        
-        auto dHdq_val = dHdq(q);
-        for (unsigned j=0; j < p.size(); j++){
-            p[j] = p[j] - a * dt * dHdq_val[j];
-        }
-    }
-    // Have to do one more step in p
-    auto dHdp_val = dHdp(p);
-    for (unsigned j=0; j < q.size(); j++){
-        q[j] = q[j] + b1 * dt * dHdp_val[j];
+        RandomAccessContainer dHdp_val = dHdp(p);
+        RandomAccessContainer dq = mult_prefactor(dHdp_val, dt * b);
+        q = add(q, dq);
+
+        RandomAccessContainer dHdq_val = dHdq(q);
+        RandomAccessContainer dp = mult_prefactor(dHdq_val, -a * dt);
+        p = add(p, dp);
     }
+    // Need to do one more step in q
+    RandomAccessContainer dHdp_val = dHdp(p);
+    RandomAccessContainer dq = mult_prefactor(dHdp_val, dt * b1);
+    q = add(q, dq);
+
     return std::make_pair(p, q);
 }
 
-template <typename RandomAccessContainer, class Func, class Policy>
+template <typename RandomAccessContainer, typename RealType, class Func, class Policy>
 std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian_imp(const RandomAccessContainer p0,
                                                                                                              const RandomAccessContainer q0,
-                                                                                                             const typename RandomAccessContainer::value_type dt,
+                                                                                                             const RealType dt,
                                                                                                              const unsigned steps,
                                                                                                              Func dHdp,
                                                                                                              Func dHdq,
                                                                                                              std::string method,
                                                                                                              const Policy& pol)
 {
-    using RealType = typename RandomAccessContainer::value_type;
     // Not sure how to make this function string nicer
     static const char* function = "boost::math::quadrature::integrate_hamiltonian(p0, q0, %1%, steps, dHdp, dHdq)";
 
@@ -208,10 +251,10 @@ std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer>
 }
 } // namespace detail
 
-template <typename RandomAccessContainer, class Func, class Policy>
+template <typename RandomAccessContainer, typename RealType, class Func, class Policy>
 std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian(const RandomAccessContainer p0,
                                                                                                          const RandomAccessContainer q0,
-                                                                                                         const typename RandomAccessContainer::value_type dt,
+                                                                                                         const RealType dt,
                                                                                                          const unsigned steps,
                                                                                                          Func dHdp,
                                                                                                          Func dHdq,
@@ -221,10 +264,10 @@ std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer>
     return detail::integrate_hamiltonian_imp(p0, q0, dt, steps, dHdp, dHdq, method, pol); 
 }
 
-template <typename RandomAccessContainer, class Func>
+template <typename RandomAccessContainer, typename RealType, class Func>
 std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian(const RandomAccessContainer p0,
                                                                                                          const RandomAccessContainer q0,
-                                                                                                         const typename RandomAccessContainer::value_type dt,
+                                                                                                         const RealType dt,
                                                                                                          const unsigned steps,
                                                                                                          Func dHdp,
                                                                                                          Func dHdq,
@@ -233,10 +276,10 @@ std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer>
     return integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, method, boost::math::policies::policy<>()); 
 }
 
-template <typename RandomAccessContainer, class Func>
+template <typename RandomAccessContainer, typename RealType, class Func>
 std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian(const RandomAccessContainer p0,
                                                                                                          const RandomAccessContainer q0,
-                                                                                                         const typename RandomAccessContainer::value_type dt,
+                                                                                                         const RealType dt,
                                                                                                          const unsigned steps,
                                                                                                          Func dHdp,
                                                                                                          Func dHdq)
diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index e4124e4aaf..121214a60e 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -26,9 +26,9 @@ using boost::math::quadrature::integrate_hamiltonian;
 
 // Equations of motion for simple harmonic oscillator
 template <class Real>
-std::vector<Real> vector_dHdp(std::vector<Real> p) { return p; }
+Real oscillator_dHdp(Real p) { return p; }
 template <class Real>
-std::vector<Real> vector_dHdq(std::vector<Real> q) { return q; }
+Real oscillator_dHdq(Real q) { return q; }
 
 // Equations of motion for simple pendulum
 template <class Real>
@@ -48,13 +48,13 @@ std::vector<Real> pendulum_vector_dHdq(std::vector<Real> q)
 template <class RealType>
 void test_invalid_parameters()
 {
-    std::vector<RealType> q0 = {1};
-    std::vector<RealType> p0 = {0};
+    RealType q0 = 1;
+    RealType p0 = 0;
     // Negative timestep
-    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(q0, p0, -0.1, 10, vector_dHdp<RealType>, vector_dHdq<RealType>), std::domain_error);
+    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(q0, p0, -0.1, 10, oscillator_dHdp<RealType>, oscillator_dHdq<RealType>), std::domain_error);
 
     // Method not in {'Y6', 'Y4', 'Y2'}
-    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(q0, p0, 0.1, 10, vector_dHdp<RealType>, vector_dHdq<RealType>, "InvalidMethod"), std::out_of_range);
+    BOOST_CHECK_THROW(boost::math::quadrature::integrate_hamiltonian(q0, p0, 0.1, 10, oscillator_dHdp<RealType>, oscillator_dHdq<RealType>, "InvalidMethod"), std::out_of_range);
 }
 
 /* Test if SHO energy fluctuations are below a given tolerance*/
@@ -67,21 +67,21 @@ void test_harmonic_oscillator(const RealType tol, const std::string method)
     RealType t_end = 100;
     unsigned int steps = t_end / dt;
 
-    std::vector<RealType> q0 = {1};
-    std::vector<RealType> p0 = {0};
+    RealType q0 = 1;
+    RealType p0 = 0;
 
-    std::vector<std::vector<RealType> > p;
-    std::vector<std::vector<RealType> > q;
+    std::vector<RealType> p;
+    std::vector<RealType> q;
 
-    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, vector_dHdp<RealType>, vector_dHdq<RealType>, method);
+    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, oscillator_dHdp<RealType>, oscillator_dHdq<RealType>, method);
 
     RealType p_val;
     RealType q_val;
     std::vector<RealType> abs_energy_error(p.size());
     for (unsigned i=0; i < p.size(); i++)
     {
-        p_val = p[i][0];
-        q_val = q[i][0];
+        p_val = p[i];
+        q_val = q[i];
 
         abs_energy_error[i] = std::abs(std::pow(p_val, 2) + std::pow(q_val, 2) - 1);
     }
@@ -128,7 +128,7 @@ BOOST_AUTO_TEST_CASE(symplectic_quadrature)
 {
     test_invalid_parameters<double>();
     
-    // SHO Tests
+    // Simple Harmonic Oscillator Tests
     test_harmonic_oscillator<double>(1e-10, "Y6");
     test_harmonic_oscillator<double>(7e-4, "Y4");
     test_harmonic_oscillator<double>(7e-4, "Y2");

From bba624144280fc373a3d4c13adb3ac396f505b91 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 12 Jul 2026 10:13:51 -0700
Subject: [PATCH 19/22] Added higher order method

---
 test/test_symplectic.cpp | 79 ++++++++++++++++++++++++++++++++++++++--
 1 file changed, 76 insertions(+), 3 deletions(-)

diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index 121214a60e..fc89c8ad0f 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -14,6 +14,7 @@
 #include <boost/math/concepts/real_concept.hpp>
 #include <boost/math/quadrature/symplectic.hpp>
 #include <boost/math/constants/constants.hpp>
+#include <iostream>
 #ifdef BOOST_HAS_FLOAT128
 #include <boost/multiprecision/complex128.hpp>
 #endif
@@ -45,6 +46,35 @@ std::vector<Real> pendulum_vector_dHdq(std::vector<Real> q)
     return q_return;
 }
 
+// Equations of motion for Henon Heiles potential
+template <class Real> 
+std::vector<Real> hh_dHdp(std::vector<Real> p) 
+{
+    return p;
+}
+
+template <class Real> 
+std::vector<Real> hh_dHdq(std::vector<Real> q)
+{
+    BOOST_MATH_STD_USING
+    std::vector<Real> dHdq(q.size());
+    dHdq.at(0) = q.at(0) + 2 * q.at(0) * q.at(1); // x + 2*xy
+    dHdq.at(1) = q.at(1) + pow(q.at(1), 2.0) - pow(q.at(1), 2.0); // y + x^2 - y^2
+    return dHdq;
+}
+
+// The energy of the system is H = 1/2 (px^2 + py^2) + 1/2(x^2 + y^2) + x^2 y - y^3/3
+template <class Real>
+Real hh_energy(std::vector<Real> p, std::vector<Real> q)
+{
+    BOOST_MATH_STD_USING
+    Real x = q[0];
+    Real y = q[1];
+    Real px = p[0];
+    Real py = p[1];
+    return 0.5 * (pow(px, 2.) + pow(py, 2.)) + 0.5 * (pow(x, 2.) + pow(y, 2.)) + pow(x, 2.) * y - pow(y, 3.) / 3.;
+}
+
 template <class RealType>
 void test_invalid_parameters()
 {
@@ -100,8 +130,8 @@ void test_pendulum(const RealType tol, const std::string method)
     RealType t_end = 100;
     unsigned int steps = t_end / dt;
 
-    std::vector<RealType> q0 = {boost::math::constants::pi<RealType>() / 2};
-    std::vector<RealType> p0 = {0};
+    std::vector<RealType> q0 = {boost::math::constants::pi<RealType>() * 0.5};
+    std::vector<RealType> p0 = {0.};
 
     std::vector<std::vector<RealType> > p;
     std::vector<std::vector<RealType> > q;
@@ -117,13 +147,51 @@ void test_pendulum(const RealType tol, const std::string method)
         q_val = q[i][0];
         // The energy of the system is H = p^2 / 2 + (1-cos(q)). With our initial condition
         // this yields H = 1. Thus, we just want p^2 / 2 - cos(q) = 0
-        abs_energy_error[i] = std::abs(std::pow(p_val, 2) / 2 - std::cos(q_val));
+        abs_energy_error[i] = std::abs(std::pow(p_val, 2.) / 2. - std::cos(q_val));
     }
 
     RealType max_error = *std::max_element(std::begin(abs_energy_error), std::end(abs_energy_error));
     BOOST_CHECK_LE(max_error, tol);
 }
 
+/* Test if SHO energy fluctuations are below a given tolerance*/
+template <class RealType>
+void test_hh_model(const RealType tol, const std::string method)
+{
+    BOOST_MATH_STD_USING
+
+    RealType dt = 0.005;
+    RealType t_end = 100.;
+    unsigned int steps = t_end / dt;
+    
+    std::vector<RealType> q0 = {0.05, 0};
+    std::vector<RealType> p0 = {0, 0.025};
+
+    RealType total_energy = hh_energy(p0, q0);
+
+    std::vector<std::vector<RealType> > p;
+    std::vector<std::vector<RealType> > q;
+
+    std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, hh_dHdp<RealType>, hh_dHdq<RealType>, method);
+
+    std::vector<RealType> p_val;
+    std::vector<RealType> q_val;
+    std::vector<RealType> abs_energy_error(p.size());
+    RealType sum = 0;
+    for (unsigned i=0; i < p.size(); i++)
+    {
+        p_val = p[i];
+        q_val = q[i];
+        
+        abs_energy_error[i] = std::abs(hh_energy(p_val, q_val) - total_energy);
+        sum += abs_energy_error[i];
+    }
+
+    std::cout << std::setprecision(16) << sum / p.size() << std::endl;
+    RealType max_error = *std::max_element(std::begin(abs_energy_error), std::end(abs_energy_error));
+    BOOST_CHECK_LE(max_error, tol);
+}
+
 BOOST_AUTO_TEST_CASE(symplectic_quadrature)
 {
     test_invalid_parameters<double>();
@@ -139,4 +207,9 @@ BOOST_AUTO_TEST_CASE(symplectic_quadrature)
     test_pendulum<double>(5e-4, "Y4");
     test_pendulum<double>(5e-4, "Y2");
     test_pendulum<double>(1e-8, "SRKNB6");
+
+    test_hh_model<double>(1e-16, "SRKNB11");
+    // test_hh_model<double>(5e-16, "Y4");
+    // test_hh_model<double>(5e-16, "Y2");
+    // test_hh_model<double>(1e-16, "SRKNB6");
 }

From 89ef4c030a75376d056d44a7aead38c5f9cfc97a Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 12 Jul 2026 10:15:03 -0700
Subject: [PATCH 20/22] Added higher order method (previous commit was wrong
 file)

---
 include/boost/math/quadrature/symplectic.hpp | 62 ++++++++++++++++++--
 test/test_symplectic.cpp                     |  2 +-
 2 files changed, 58 insertions(+), 6 deletions(-)

diff --git a/include/boost/math/quadrature/symplectic.hpp b/include/boost/math/quadrature/symplectic.hpp
index 310d23c46a..241ceeae69 100644
--- a/include/boost/math/quadrature/symplectic.hpp
+++ b/include/boost/math/quadrature/symplectic.hpp
@@ -84,7 +84,7 @@ std::pair<RandomAccessContainer, RandomAccessContainer> second_order_yoshida(con
 
     // Half step in q
     RandomAccessContainer dHdp_val = dHdp(p);
-    RandomAccessContainer dq = mult_prefactor(dHdp_val, dt / 2);
+    RandomAccessContainer dq = mult_prefactor(dHdp_val, 0.5 * dt);
     q = add(q, dq);
 
     // Full step in p
@@ -94,7 +94,7 @@ std::pair<RandomAccessContainer, RandomAccessContainer> second_order_yoshida(con
 
     // Half step in q
     dHdp_val = dHdp(p);
-    dq = mult_prefactor(dHdp_val, dt / 2);
+    dq = mult_prefactor(dHdp_val, 0.5 * dt);
     q = add(q, dq);
 
     return std::make_pair(p, q);
@@ -113,8 +113,8 @@ std::pair<RandomAccessContainer, RandomAccessContainer> fourth_order_yoshida(con
     RandomAccessContainer q = q0;
 
     // RealType x0 = -(std::pow(2, 1/3) / (2 - std::pow(2, 1/3)));
-    RealType x1 = 1 / (2 - std::pow(2, 1/3));
-    RealType x0 = 1 - 2 * x1; 
+    RealType x1 = 1. / (2. - std::cbrt(2.));
+    RealType x0 = 1. - 2. * x1; 
 
     std::vector<RealType> weights = { x1, x0, x1 };
 
@@ -203,6 +203,57 @@ std::pair<RandomAccessContainer, RandomAccessContainer> SRKN_b_order_6(const Ran
     return std::make_pair(p, q);
 }
 
+template <typename RandomAccessContainer, typename RealType, class Func>
+std::pair<RandomAccessContainer, RandomAccessContainer> SRKN_b_order_11(const RandomAccessContainer p0,
+                                                                       const RandomAccessContainer q0,
+                                                                       const RealType dt, 
+                                                                       Func dHdp,
+                                                                       Func dHdq)
+{ // This method implements SRKN_b^11 in Table 3 here 
+  // https://www.sciencedirect.com/science/article/pii/S0377042701004927
+
+    RandomAccessContainer p = p0;
+    RandomAccessContainer q = q0;
+    
+    RealType b1 = static_cast<RealType>(0.0414649985182624);
+    RealType b2 = static_cast<RealType>(0.198128671918067);
+    RealType b3 = static_cast<RealType>(-0.0400061921041533);
+    RealType b4 = static_cast<RealType>(0.0752539843015807);
+    RealType b5 = static_cast<RealType>(-0.0115113874206879);
+    RealType b6 = 0.5 - (b1 + b2 + b3 + b4 + b5);
+
+    RealType a1 = static_cast<RealType>(0.123229775946271);
+    RealType a2 = static_cast<RealType>(0.290553797799558);
+    RealType a3 = static_cast<RealType>(-0.127049212625417);
+    RealType a4 = static_cast<RealType>(-0.246331761062075);
+    RealType a5 = static_cast<RealType>(0.357208872795928);
+    RealType a6 = 1. - 2. * (a1 + a2 + a3 + a4 + a5);
+
+    std::vector<RealType> b_weights = {b1, b2, b3, b4, b5, b6, b6, b5, b4, b3, b2};
+    std::vector<RealType> a_weights = {a1, a2, a3, a4, a5, a6, a5, a4, a3, a2, a1};
+    
+    RealType a, b;
+    for (unsigned int i=0; i < b_weights.size(); i++)
+    {
+        b = b_weights[i];
+        a = a_weights[i];
+        
+        RandomAccessContainer dHdp_val = dHdp(p);
+        RandomAccessContainer dq = mult_prefactor(dHdp_val, dt * b);
+        q = add(q, dq);
+
+        RandomAccessContainer dHdq_val = dHdq(q);
+        RandomAccessContainer dp = mult_prefactor(dHdq_val, -a * dt);
+        p = add(p, dp);
+    }
+    // Need to do one more step in q
+    RandomAccessContainer dHdp_val = dHdp(p);
+    RandomAccessContainer dq = mult_prefactor(dHdp_val, dt * b1);
+    q = add(q, dq);
+
+    return std::make_pair(p, q);
+}
+
 template <typename RandomAccessContainer, typename RealType, class Func, class Policy>
 std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian_imp(const RandomAccessContainer p0,
                                                                                                              const RandomAccessContainer q0,
@@ -231,7 +282,8 @@ std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer>
     std::map<std::string, stepperType> m{{"Y6", sixth_order_yoshida}, 
                                          {"Y4", fourth_order_yoshida}, 
                                          {"Y2", second_order_yoshida},
-                                         {"SRKNB6", SRKN_b_order_6}};
+                                         {"SRKNB6", SRKN_b_order_6},
+                                         {"SRKNB11", SRKN_b_order_11}};
     stepperType stepper = m.at(method);
 
     std::vector<RandomAccessContainer> p(steps);
diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index fc89c8ad0f..99965f2cbc 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -208,7 +208,7 @@ BOOST_AUTO_TEST_CASE(symplectic_quadrature)
     test_pendulum<double>(5e-4, "Y2");
     test_pendulum<double>(1e-8, "SRKNB6");
 
-    test_hh_model<double>(1e-16, "SRKNB11");
+    // test_hh_model<double>(1e-16, "SRKNB11");
     // test_hh_model<double>(5e-16, "Y4");
     // test_hh_model<double>(5e-16, "Y2");
     // test_hh_model<double>(1e-16, "SRKNB6");

From 7ab85402f77bc7f611538b1ff03a793503d64295 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 12 Jul 2026 10:18:39 -0700
Subject: [PATCH 21/22] Added order 11 tests

---
 test/test_symplectic.cpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index 99965f2cbc..dc94c65580 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -201,12 +201,14 @@ BOOST_AUTO_TEST_CASE(symplectic_quadrature)
     test_harmonic_oscillator<double>(7e-4, "Y4");
     test_harmonic_oscillator<double>(7e-4, "Y2");
     test_harmonic_oscillator<double>(1e-11, "SRKNB6");
+    test_harmonic_oscillator<double>(2e-14, "SRKNB11");
 
     // Pendulum Tests
     test_pendulum<double>(1e-10, "Y6");
     test_pendulum<double>(5e-4, "Y4");
     test_pendulum<double>(5e-4, "Y2");
     test_pendulum<double>(1e-8, "SRKNB6");
+    test_pendulum<double>(1e-10, "SRKNB11");
 
     // test_hh_model<double>(1e-16, "SRKNB11");
     // test_hh_model<double>(5e-16, "Y4");

From f304b0ac5427c31b8d717af507fb32443c77c8e7 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 12 Jul 2026 10:37:50 -0700
Subject: [PATCH 22/22] Fixed hh test derivatives

---
 test/test_symplectic.cpp | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/test/test_symplectic.cpp b/test/test_symplectic.cpp
index dc94c65580..2f2e8c3212 100644
--- a/test/test_symplectic.cpp
+++ b/test/test_symplectic.cpp
@@ -58,8 +58,8 @@ std::vector<Real> hh_dHdq(std::vector<Real> q)
 {
     BOOST_MATH_STD_USING
     std::vector<Real> dHdq(q.size());
-    dHdq.at(0) = q.at(0) + 2 * q.at(0) * q.at(1); // x + 2*xy
-    dHdq.at(1) = q.at(1) + pow(q.at(1), 2.0) - pow(q.at(1), 2.0); // y + x^2 - y^2
+    dHdq[0] = q[0] + 2.0 * q[0] * q[1]; // x + 2*xy
+    dHdq[1] = q[1] + pow(q[0], 2.0) - pow(q[1], 2.0); // y + x^2 - y^2
     return dHdq;
 }
 
@@ -164,8 +164,8 @@ void test_hh_model(const RealType tol, const std::string method)
     RealType t_end = 100.;
     unsigned int steps = t_end / dt;
     
-    std::vector<RealType> q0 = {0.05, 0};
-    std::vector<RealType> p0 = {0, 0.025};
+    std::vector<RealType> q0 = {0.5, 0};
+    std::vector<RealType> p0 = {0, 0.25};
 
     RealType total_energy = hh_energy(p0, q0);
 
@@ -210,8 +210,9 @@ BOOST_AUTO_TEST_CASE(symplectic_quadrature)
     test_pendulum<double>(1e-8, "SRKNB6");
     test_pendulum<double>(1e-10, "SRKNB11");
 
-    // test_hh_model<double>(1e-16, "SRKNB11");
-    // test_hh_model<double>(5e-16, "Y4");
-    // test_hh_model<double>(5e-16, "Y2");
-    // test_hh_model<double>(1e-16, "SRKNB6");
+    test_hh_model<double>(1e-14, "SRKNB11");
+    test_hh_model<double>(1e-14, "Y6");
+    test_hh_model<double>(5e-11, "Y4");
+    test_hh_model<double>(5e-6, "Y2");
+    test_hh_model<double>(1e-13, "SRKNB6");
 }