From 5901a46c01bbadac922d077c75b5cc910b1b473b Mon Sep 17 00:00:00 2001
From: Francesco Biscani <bluescarni@gmail.com>
Date: Wed, 6 Jan 2021 15:18:05 +0100
Subject: [PATCH 1/4] Fix benchmark code.

---
 benchmark/n_body_creation.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmark/n_body_creation.cpp b/benchmark/n_body_creation.cpp
index 1c2de46ae..b2d55192d 100644
--- a/benchmark/n_body_creation.cpp
+++ b/benchmark/n_body_creation.cpp
@@ -66,7 +66,7 @@ int main(int argc, char *argv[])
 
     auto counter = 0u;
     for (const auto &ex : ta.get_decomposition()) {
-        std::cout << "u_" << counter++ << " = " << ex << '\n';
+        std::cout << "u_" << counter++ << " = " << ex.first << '\n';
     }
 
     std::cout << "Construction time: " << elapsed << "ms\n";

From 2854e899f582706a02015a4debdde8900474029e Mon Sep 17 00:00:00 2001
From: Francesco Biscani <bluescarni@gmail.com>
Date: Wed, 6 Jan 2021 15:18:28 +0100
Subject: [PATCH 2/4] Use operator== instead of operator= for the prime(x)
 syntax.

---
 README.md                     |  2 +-
 benchmark/stiff_equation.cpp  |  2 +-
 benchmark/taylor_jl_01.cpp    |  2 +-
 doc/tut_adaptive.rst          |  2 +-
 include/heyoka/expression.hpp |  2 +-
 src/expression.cpp            |  4 +--
 src/mascon.cpp                | 12 +++----
 src/nbody.cpp                 | 60 +++++++++++++++++------------------
 test/back_and_forth.cpp       |  2 +-
 test/e3bp.cpp                 |  5 +--
 test/llvm_state.cpp           |  2 +-
 test/one_body.cpp             |  4 +--
 test/taylor_const_sys.cpp     | 34 ++++++++++----------
 test/taylor_no_decomp_sys.cpp | 12 +++----
 tutorial/adaptive_basic.cpp   |  2 +-
 tutorial/adaptive_opt.cpp     |  2 +-
 tutorial/pendulum.cpp         |  2 +-
 17 files changed, 76 insertions(+), 75 deletions(-)

diff --git a/README.md b/README.md
index e3f4f7722..94f557a16 100644
--- a/README.md
+++ b/README.md
@@ -52,7 +52,7 @@ int main()
     auto ta = taylor_adaptive<double>{// Definition of the ODE system:
                                       // x' = v
                                       // v' = -9.8 * sin(x)
-                                      {prime(x) = v, prime(v) = -9.8 * sin(x)},
+                                      {prime(x) == v, prime(v) == -9.8 * sin(x)},
                                       // Initial conditions
                                       // for x and v.
                                       {0.05, 0.025}};
diff --git a/benchmark/stiff_equation.cpp b/benchmark/stiff_equation.cpp
index f501fbb85..0a2441529 100644
--- a/benchmark/stiff_equation.cpp
+++ b/benchmark/stiff_equation.cpp
@@ -23,7 +23,7 @@ int main()
     for (unsigned i = 1u; i < 6; ++i) {
         // Van der pool stiff equation
         auto mu = expression{number{std::pow(10, i)}};
-        taylor_adaptive<double> taylor{{prime(y1) = y2, prime(y2) = mu * (1_dbl - y1 * y1) * y2 - y1}, {2., 0.}};
+        taylor_adaptive<double> taylor{{prime(y1) == y2, prime(y2) == mu * (1_dbl - y1 * y1) * y2 - y1}, {2., 0.}};
         std::cout << "Stiffness param " << mu << ": ";
 
         auto start = high_resolution_clock::now();
diff --git a/benchmark/taylor_jl_01.cpp b/benchmark/taylor_jl_01.cpp
index 37d55af59..9aa2b74a9 100644
--- a/benchmark/taylor_jl_01.cpp
+++ b/benchmark/taylor_jl_01.cpp
@@ -32,7 +32,7 @@ int main()
 
     const auto x0 = 3.0;
 
-    taylor_adaptive<double> ta{{prime(x) = x * x}, {x0}, kw::high_accuracy = true, kw::tol = 1e-18};
+    taylor_adaptive<double> ta{{prime(x) == x * x}, {x0}, kw::high_accuracy = true, kw::tol = 1e-18};
 
     const auto final_time = 0.3333333329479479;
 
diff --git a/doc/tut_adaptive.rst b/doc/tut_adaptive.rst
index 07b5a013d..d485c8a5f 100644
--- a/doc/tut_adaptive.rst
+++ b/doc/tut_adaptive.rst
@@ -46,7 +46,7 @@ In this case, we use the ``double`` type, meaning that the integration
 will be carried out in double precision.
 
 As (mandatory) construction arguments, we pass in the system of
-differential equations using the syntax ``prime(x) = ...``, and a set
+differential equations using the syntax ``prime(x) == ...``, and a set
 of initial conditions for ``x`` and ``v`` respectively.
 
 Let's try to print to screen the integrator object:
diff --git a/include/heyoka/expression.hpp b/include/heyoka/expression.hpp
index 5ce31cfcf..0fa413df9 100644
--- a/include/heyoka/expression.hpp
+++ b/include/heyoka/expression.hpp
@@ -144,7 +144,7 @@ struct HEYOKA_DLL_PUBLIC prime_wrapper {
     prime_wrapper &operator=(prime_wrapper &&) noexcept;
     ~prime_wrapper();
 
-    std::pair<expression, expression> operator=(expression) &&;
+    std::pair<expression, expression> operator==(expression) const;
 };
 
 } // namespace detail
diff --git a/src/expression.cpp b/src/expression.cpp
index 3d769ad36..f478e5a33 100644
--- a/src/expression.cpp
+++ b/src/expression.cpp
@@ -133,9 +133,9 @@ prime_wrapper &prime_wrapper::operator=(prime_wrapper &&) noexcept = default;
 
 prime_wrapper::~prime_wrapper() = default;
 
-std::pair<expression, expression> prime_wrapper::operator=(expression e) &&
+std::pair<expression, expression> prime_wrapper::operator==(expression e) const
 {
-    return std::pair{expression{variable{std::move(m_str)}}, std::move(e)};
+    return std::pair{expression{variable{m_str}}, std::move(e)};
 }
 
 } // namespace detail
diff --git a/src/mascon.cpp b/src/mascon.cpp
index b0d83780f..61dcf4b88 100644
--- a/src/mascon.cpp
+++ b/src/mascon.cpp
@@ -58,12 +58,12 @@ make_mascon_system_impl(expression Gconst, std::vector<std::vector<expression>>
 
     // Assembling the return vector containing l.h.s. and r.h.s. (note the fundamental use of pairwise_sum for
     // efficiency and to allow compact mode to do his job)
-    retval.push_back(prime(x) = vx);
-    retval.push_back(prime(y) = vy);
-    retval.push_back(prime(z) = vz);
-    retval.push_back(prime(vx) = pairwise_sum(x_acc) - centripetal_x - coriolis_x);
-    retval.push_back(prime(vy) = pairwise_sum(y_acc) - centripetal_y - coriolis_y);
-    retval.push_back(prime(vz) = pairwise_sum(z_acc) - centripetal_z - coriolis_z);
+    retval.push_back(prime(x) == vx);
+    retval.push_back(prime(y) == vy);
+    retval.push_back(prime(z) == vz);
+    retval.push_back(prime(vx) == pairwise_sum(x_acc) - centripetal_x - coriolis_x);
+    retval.push_back(prime(vy) == pairwise_sum(y_acc) - centripetal_y - coriolis_y);
+    retval.push_back(prime(vz) == pairwise_sum(z_acc) - centripetal_z - coriolis_z);
 
     return retval;
 }
diff --git a/src/nbody.cpp b/src/nbody.cpp
index 2bb312f5d..9b4fe3dcd 100644
--- a/src/nbody.cpp
+++ b/src/nbody.cpp
@@ -106,9 +106,9 @@ std::vector<std::pair<expression, expression>> make_nbody_sys_fixed_masses(std::
         // on all particles.
         for (std::uint32_t i = 0; i < n_fc_massive; ++i) {
             // r' = v.
-            retval.push_back(prime(x_vars[i]) = vx_vars[i]);
-            retval.push_back(prime(y_vars[i]) = vy_vars[i]);
-            retval.push_back(prime(z_vars[i]) = vz_vars[i]);
+            retval.push_back(prime(x_vars[i]) == vx_vars[i]);
+            retval.push_back(prime(y_vars[i]) == vy_vars[i]);
+            retval.push_back(prime(z_vars[i]) == vz_vars[i]);
 
             // Compute the total acceleration on body i,
             // accumulating the result in x/y/z_acc. Part of the
@@ -151,29 +151,29 @@ std::vector<std::pair<expression, expression>> make_nbody_sys_fixed_masses(std::
             }
 
             // Add the expressions of the accelerations to the system.
-            retval.push_back(prime(vx_vars[i]) = pairwise_sum(x_acc[i]));
-            retval.push_back(prime(vy_vars[i]) = pairwise_sum(y_acc[i]));
-            retval.push_back(prime(vz_vars[i]) = pairwise_sum(z_acc[i]));
+            retval.push_back(prime(vx_vars[i]) == pairwise_sum(x_acc[i]));
+            retval.push_back(prime(vy_vars[i]) == pairwise_sum(y_acc[i]));
+            retval.push_back(prime(vz_vars[i]) == pairwise_sum(z_acc[i]));
         }
 
         // All the accelerations on the massless particles
         // have already been accumulated in the loop above.
         // We just need to perform the pairwise sums on x/y/z_acc.
         for (std::uint32_t i = n_fc_massive; i < n; ++i) {
-            retval.push_back(prime(x_vars[i]) = vx_vars[i]);
-            retval.push_back(prime(y_vars[i]) = vy_vars[i]);
-            retval.push_back(prime(z_vars[i]) = vz_vars[i]);
+            retval.push_back(prime(x_vars[i]) == vx_vars[i]);
+            retval.push_back(prime(y_vars[i]) == vy_vars[i]);
+            retval.push_back(prime(z_vars[i]) == vz_vars[i]);
 
-            retval.push_back(prime(vx_vars[i]) = pairwise_sum(x_acc[i]));
-            retval.push_back(prime(vy_vars[i]) = pairwise_sum(y_acc[i]));
-            retval.push_back(prime(vz_vars[i]) = pairwise_sum(z_acc[i]));
+            retval.push_back(prime(vx_vars[i]) == pairwise_sum(x_acc[i]));
+            retval.push_back(prime(vy_vars[i]) == pairwise_sum(y_acc[i]));
+            retval.push_back(prime(vz_vars[i]) == pairwise_sum(z_acc[i]));
         }
     } else {
         for (std::uint32_t i = 0; i < n; ++i) {
             // r' = v.
-            retval.push_back(prime(x_vars[i]) = vx_vars[i]);
-            retval.push_back(prime(y_vars[i]) = vy_vars[i]);
-            retval.push_back(prime(z_vars[i]) = vz_vars[i]);
+            retval.push_back(prime(x_vars[i]) == vx_vars[i]);
+            retval.push_back(prime(y_vars[i]) == vy_vars[i]);
+            retval.push_back(prime(z_vars[i]) == vz_vars[i]);
 
             // Compute the total acceleration on body i,
             // accumulating the result in x/y/z_acc. Part of the
@@ -215,9 +215,9 @@ std::vector<std::pair<expression, expression>> make_nbody_sys_fixed_masses(std::
             }
 
             // Add the expressions of the accelerations to the system.
-            retval.push_back(prime(vx_vars[i]) = pairwise_sum(x_acc[i]));
-            retval.push_back(prime(vy_vars[i]) = pairwise_sum(y_acc[i]));
-            retval.push_back(prime(vz_vars[i]) = pairwise_sum(z_acc[i]));
+            retval.push_back(prime(vx_vars[i]) == pairwise_sum(x_acc[i]));
+            retval.push_back(prime(vy_vars[i]) == pairwise_sum(y_acc[i]));
+            retval.push_back(prime(vz_vars[i]) == pairwise_sum(z_acc[i]));
         }
     }
 
@@ -266,9 +266,9 @@ std::vector<std::pair<expression, expression>> make_nbody_sys_par_masses(std::ui
     // on all particles.
     for (std::uint32_t i = 0; i < n_massive; ++i) {
         // r' = v.
-        retval.push_back(prime(x_vars[i]) = vx_vars[i]);
-        retval.push_back(prime(y_vars[i]) = vy_vars[i]);
-        retval.push_back(prime(z_vars[i]) = vz_vars[i]);
+        retval.push_back(prime(x_vars[i]) == vx_vars[i]);
+        retval.push_back(prime(y_vars[i]) == vy_vars[i]);
+        retval.push_back(prime(z_vars[i]) == vz_vars[i]);
 
         // Compute the total acceleration on body i,
         // accumulating the result in x/y/z_acc. Part of the
@@ -311,22 +311,22 @@ std::vector<std::pair<expression, expression>> make_nbody_sys_par_masses(std::ui
         }
 
         // Add the expressions of the accelerations to the system.
-        retval.push_back(prime(vx_vars[i]) = pairwise_sum(x_acc[i]));
-        retval.push_back(prime(vy_vars[i]) = pairwise_sum(y_acc[i]));
-        retval.push_back(prime(vz_vars[i]) = pairwise_sum(z_acc[i]));
+        retval.push_back(prime(vx_vars[i]) == pairwise_sum(x_acc[i]));
+        retval.push_back(prime(vy_vars[i]) == pairwise_sum(y_acc[i]));
+        retval.push_back(prime(vz_vars[i]) == pairwise_sum(z_acc[i]));
     }
 
     // All the accelerations on the massless particles
     // have already been accumulated in the loop above.
     // We just need to perform the pairwise sums on x/y/z_acc.
     for (std::uint32_t i = n_massive; i < n; ++i) {
-        retval.push_back(prime(x_vars[i]) = vx_vars[i]);
-        retval.push_back(prime(y_vars[i]) = vy_vars[i]);
-        retval.push_back(prime(z_vars[i]) = vz_vars[i]);
+        retval.push_back(prime(x_vars[i]) == vx_vars[i]);
+        retval.push_back(prime(y_vars[i]) == vy_vars[i]);
+        retval.push_back(prime(z_vars[i]) == vz_vars[i]);
 
-        retval.push_back(prime(vx_vars[i]) = pairwise_sum(x_acc[i]));
-        retval.push_back(prime(vy_vars[i]) = pairwise_sum(y_acc[i]));
-        retval.push_back(prime(vz_vars[i]) = pairwise_sum(z_acc[i]));
+        retval.push_back(prime(vx_vars[i]) == pairwise_sum(x_acc[i]));
+        retval.push_back(prime(vy_vars[i]) == pairwise_sum(y_acc[i]));
+        retval.push_back(prime(vz_vars[i]) == pairwise_sum(z_acc[i]));
     }
 
     return retval;
diff --git a/test/back_and_forth.cpp b/test/back_and_forth.cpp
index d424c81ec..14a892493 100644
--- a/test/back_and_forth.cpp
+++ b/test/back_and_forth.cpp
@@ -43,7 +43,7 @@ TEST_CASE("pendulum")
     auto tester = [](auto fp_x, bool high_accuracy, bool compact_mode) {
         using fp_t = decltype(fp_x);
 
-        taylor_adaptive<fp_t> ta{{prime("th"_var) = "v"_var, prime("v"_var) = -9.8_dbl / 1.5_dbl * sin("th"_var)},
+        taylor_adaptive<fp_t> ta{{prime("th"_var) == "v"_var, prime("v"_var) == -9.8_dbl / 1.5_dbl * sin("th"_var)},
                                  {fp_t{0.05}, fp_t{0.025}},
                                  kw::high_accuracy = high_accuracy,
                                  kw::compact_mode = compact_mode};
diff --git a/test/e3bp.cpp b/test/e3bp.cpp
index 60229c119..ff0691f8b 100644
--- a/test/e3bp.cpp
+++ b/test/e3bp.cpp
@@ -69,8 +69,9 @@ TEST_CASE("e3bp")
         = cart_to_ham_ec(std::vector{1.20793759666736, -0.493320558636725, 1.19760678594565, -0.498435147674914,
                                      0.548228167205306, 0.49662691628363});
 
-    taylor_adaptive_dbl tad{{prime(xi) = diff(ham, pxi), prime(eta) = diff(ham, peta), prime(phi) = diff(ham, pphi),
-                             prime(pxi) = -diff(ham, xi), prime(peta) = -diff(ham, eta), prime(pphi) = -diff(ham, phi)},
+    taylor_adaptive_dbl tad{{prime(xi) == diff(ham, pxi), prime(eta) == diff(ham, peta), prime(phi) == diff(ham, pphi),
+                             prime(pxi) == -diff(ham, xi), prime(peta) == -diff(ham, eta),
+                             prime(pphi) == -diff(ham, phi)},
                             init_state};
 
     const auto &st = tad.get_state();
diff --git a/test/llvm_state.cpp b/test/llvm_state.cpp
index 1adb4bffb..8e976cf75 100644
--- a/test/llvm_state.cpp
+++ b/test/llvm_state.cpp
@@ -22,7 +22,7 @@ TEST_CASE("basic")
 
     llvm_state s;
     auto [x, y] = make_vars("x", "y");
-    taylor_add_jet_dbl(s, "foo", {prime(x) = y, prime(y) = (1_dbl - x * x) * y - x}, 21, 1, true, false);
+    taylor_add_jet_dbl(s, "foo", {prime(x) == y, prime(y) == (1_dbl - x * x) * y - x}, 21, 1, true, false);
 
     std::cout << s.get_ir() << '\n';
 }
diff --git a/test/one_body.cpp b/test/one_body.cpp
index 604ca61b4..216c58559 100644
--- a/test/one_body.cpp
+++ b/test/one_body.cpp
@@ -29,8 +29,8 @@ TEST_CASE("one_body")
 
     auto [s_x, s_v] = kep_to_cart({1., .5, .0001, 0., 0., 0.}, 1.);
 
-    taylor_adaptive<double> ta{{prime(x) = vx, prime(y) = vy, prime(z) = vz, prime(vx) = -x * r_m3,
-                                prime(vy) = -y * r_m3, prime(vz) = -z * r_m3},
+    taylor_adaptive<double> ta{{prime(x) == vx, prime(y) == vy, prime(z) == vz, prime(vx) == -x * r_m3,
+                                prime(vy) == -y * r_m3, prime(vz) == -z * r_m3},
                                {s_x[0], s_x[1], s_x[2], s_v[0], s_v[1], s_v[2]},
                                kw::high_accuracy = true,
                                kw::tol = 1e-18};
diff --git a/test/taylor_const_sys.cpp b/test/taylor_const_sys.cpp
index 01813d0fb..e0753781b 100644
--- a/test/taylor_const_sys.cpp
+++ b/test/taylor_const_sys.cpp
@@ -92,8 +92,8 @@ TEST_CASE("taylor const sys")
             llvm_state s{kw::opt_level = opt_level};
 
             taylor_add_jet<fp_t>(s, "jet",
-                                 {prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
-                                  prime(z) = expression{number{fp_t{0}}}},
+                                 {prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
+                                  prime(z) == expression{number{fp_t{0}}}},
                                  1, 1, high_accuracy, compact_mode);
 
             s.compile();
@@ -118,8 +118,8 @@ TEST_CASE("taylor const sys")
 
             taylor_add_jet<fp_t>(
                 s, "jet",
-                {prime(x) = par[0], prime(y) = expression{number{fp_t{-2}}}, prime(z) = expression{number{fp_t{0}}}}, 1,
-                1, high_accuracy, compact_mode);
+                {prime(x) == par[0], prime(y) == expression{number{fp_t{-2}}}, prime(z) == expression{number{fp_t{0}}}},
+                1, 1, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -144,8 +144,8 @@ TEST_CASE("taylor const sys")
             llvm_state s{kw::opt_level = opt_level};
 
             taylor_add_jet<fp_t>(s, "jet",
-                                 {prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
-                                  prime(z) = expression{number{fp_t{0}}}},
+                                 {prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
+                                  prime(z) == expression{number{fp_t{0}}}},
                                  1, 2, high_accuracy, compact_mode);
 
             s.compile();
@@ -181,8 +181,8 @@ TEST_CASE("taylor const sys")
 
             taylor_add_jet<fp_t>(
                 s, "jet",
-                {prime(x) = expression{number{fp_t{1}}}, prime(y) = par[1], prime(z) = expression{number{fp_t{0}}}}, 1,
-                2, high_accuracy, compact_mode);
+                {prime(x) == expression{number{fp_t{1}}}, prime(y) == par[1], prime(z) == expression{number{fp_t{0}}}},
+                1, 2, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -218,8 +218,8 @@ TEST_CASE("taylor const sys")
             llvm_state s{kw::opt_level = opt_level};
 
             taylor_add_jet<fp_t>(s, "jet",
-                                 {prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
-                                  prime(z) = expression{number{fp_t{0}}}},
+                                 {prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
+                                  prime(z) == expression{number{fp_t{0}}}},
                                  2, 1, high_accuracy, compact_mode);
 
             s.compile();
@@ -246,8 +246,8 @@ TEST_CASE("taylor const sys")
             llvm_state s{kw::opt_level = opt_level};
 
             taylor_add_jet<fp_t>(s, "jet",
-                                 {prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
-                                  prime(z) = expression{number{fp_t{0}}}},
+                                 {prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
+                                  prime(z) == expression{number{fp_t{0}}}},
                                  2, 2, high_accuracy, compact_mode);
 
             s.compile();
@@ -291,8 +291,8 @@ TEST_CASE("taylor const sys")
             llvm_state s{kw::opt_level = opt_level};
 
             taylor_add_jet<fp_t>(s, "jet",
-                                 {prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
-                                  prime(z) = expression{number{fp_t{0}}}},
+                                 {prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
+                                  prime(z) == expression{number{fp_t{0}}}},
                                  3, 3, high_accuracy, compact_mode);
 
             s.compile();
@@ -336,7 +336,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) = par[0], prime(y) = par[1], prime(z) = par[2]}, 3, 3,
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) == par[0], prime(y) == par[1], prime(z) == par[2]}, 3, 3,
                                  high_accuracy, compact_mode);
 
             s.compile();
@@ -380,8 +380,8 @@ TEST_CASE("taylor const sys")
         }
 
         // Do the batch/scalar comparison.
-        compare_batch_scalar<fp_t>({prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
-                                    prime(z) = expression{number{fp_t{0}}}},
+        compare_batch_scalar<fp_t>({prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
+                                    prime(z) == expression{number{fp_t{0}}}},
                                    opt_level, high_accuracy, compact_mode);
     };
 
diff --git a/test/taylor_no_decomp_sys.cpp b/test/taylor_no_decomp_sys.cpp
index 0d6fae57d..da36ba774 100644
--- a/test/taylor_no_decomp_sys.cpp
+++ b/test/taylor_no_decomp_sys.cpp
@@ -90,7 +90,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) = y, prime(y) = x}, 1, 1, high_accuracy, compact_mode);
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) == y, prime(y) == x}, 1, 1, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -110,7 +110,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) = y, prime(y) = x}, 1, 2, high_accuracy, compact_mode);
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) == y, prime(y) == x}, 1, 2, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -137,7 +137,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) = y, prime(y) = x}, 2, 1, high_accuracy, compact_mode);
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) == y, prime(y) == x}, 2, 1, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -159,7 +159,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) = y, prime(y) = x}, 2, 2, high_accuracy, compact_mode);
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) == y, prime(y) == x}, 2, 2, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -192,7 +192,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) = y, prime(y) = x}, 3, 3, high_accuracy, compact_mode);
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) == y, prime(y) == x}, 3, 3, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -237,7 +237,7 @@ TEST_CASE("taylor const sys")
         }
 
         // Do the batch/scalar comparison.
-        compare_batch_scalar<fp_t>({prime(x) = y, prime(y) = x}, opt_level, high_accuracy, compact_mode);
+        compare_batch_scalar<fp_t>({prime(x) == y, prime(y) == x}, opt_level, high_accuracy, compact_mode);
     };
 
     for (auto cm : {true, false}) {
diff --git a/tutorial/adaptive_basic.cpp b/tutorial/adaptive_basic.cpp
index 2797f0841..ce7c09216 100644
--- a/tutorial/adaptive_basic.cpp
+++ b/tutorial/adaptive_basic.cpp
@@ -22,7 +22,7 @@ int main()
     auto ta = taylor_adaptive<double>{// Definition of the ODE system:
                                       // x' = v
                                       // v' = -9.8 * sin(x)
-                                      {prime(x) = v, prime(v) = -9.8 * sin(x)},
+                                      {prime(x) == v, prime(v) == -9.8 * sin(x)},
                                       // Initial conditions
                                       // for x and v.
                                       {0.05, 0.025}};
diff --git a/tutorial/adaptive_opt.cpp b/tutorial/adaptive_opt.cpp
index 094bb91a8..afcad4516 100644
--- a/tutorial/adaptive_opt.cpp
+++ b/tutorial/adaptive_opt.cpp
@@ -25,7 +25,7 @@ int main()
     auto ta = taylor_adaptive<double>{// Definition of the ODE system:
                                       // x' = v
                                       // v' = -9.8 * sin(x)
-                                      {prime(x) = v, prime(v) = -9.8 * sin(x)},
+                                      {prime(x) == v, prime(v) == -9.8 * sin(x)},
                                       // Initial conditions
                                       // for x and v.
                                       {0.05, 0.025},
diff --git a/tutorial/pendulum.cpp b/tutorial/pendulum.cpp
index 9266a8f3e..fe973746c 100644
--- a/tutorial/pendulum.cpp
+++ b/tutorial/pendulum.cpp
@@ -22,7 +22,7 @@ int main()
     auto ta = taylor_adaptive<double>{// Definition of the ODE system:
                                       // x' = v
                                       // v' = -9.8 * sin(x)
-                                      {prime(x) = v, prime(v) = -9.8 * sin(x)},
+                                      {prime(x) == v, prime(v) == -9.8 * sin(x)},
                                       // Initial conditions
                                       // for x and v.
                                       {0.05, 0.025}};

From 2754bb60d0c05e4e6dc237ea0e4bcd97fa850e4d Mon Sep 17 00:00:00 2001
From: Francesco Biscani <bluescarni@gmail.com>
Date: Wed, 6 Jan 2021 15:30:16 +0100
Subject: [PATCH 3/4] Implement a pair of trivial simplifications expression's
 operator.

---
 src/expression.cpp  | 8 ++++++++
 test/expression.cpp | 8 ++++++++
 2 files changed, 16 insertions(+)

diff --git a/src/expression.cpp b/src/expression.cpp
index f478e5a33..72bb840e3 100644
--- a/src/expression.cpp
+++ b/src/expression.cpp
@@ -224,6 +224,10 @@ expression operator+(expression e1, expression e2)
                 // e1 + 0 = e1.
                 return expression{std::forward<decltype(v1)>(v1)};
             }
+            if (std::visit([](const auto &x) { return x < 0; }, v2.value())) {
+                // e1 + -x = e1 - x.
+                return expression{std::forward<decltype(v1)>(v1)} - expression{-std::forward<decltype(v2)>(v2)};
+            }
             // NOTE: fall through the standard case if e2 is not zero.
         }
 
@@ -257,6 +261,10 @@ expression operator-(expression e1, expression e2)
                 // e1 - 0 = e1.
                 return expression{std::forward<decltype(v1)>(v1)};
             }
+            if (std::visit([](const auto &x) { return x < 0; }, v2.value())) {
+                // e1 - -x = e1 + x.
+                return expression{std::forward<decltype(v1)>(v1)} + expression{-std::forward<decltype(v2)>(v2)};
+            }
             // NOTE: fall through the standard case if e2 is not zero.
         }
 
diff --git a/test/expression.cpp b/test/expression.cpp
index 1b8afc8b2..639bdc51c 100644
--- a/test/expression.cpp
+++ b/test/expression.cpp
@@ -396,3 +396,11 @@ TEST_CASE("get_param_size")
     REQUIRE(get_param_size(par[122] + par[123]) == 124u);
     REQUIRE(get_param_size(par[500] - sin(cos(par[1] + "y"_var) + par[4])) == 501u);
 }
+
+TEST_CASE("binary simpls")
+{
+    auto [x, y] = make_vars("x", "y");
+
+    REQUIRE(x + -1. == x - 1.);
+    REQUIRE(y - -1. == y + 1.);
+}

From f0e9d985eaa0c16fc5f21d34dce404e063368a05 Mon Sep 17 00:00:00 2001
From: Francesco Biscani <bluescarni@gmail.com>
Date: Wed, 6 Jan 2021 20:28:15 +0100
Subject: [PATCH 4/4] Revert "Use operator== instead of operator= for the
 prime(x) syntax."

This reverts commit 2854e899f582706a02015a4debdde8900474029e.
---
 README.md                     |  2 +-
 benchmark/stiff_equation.cpp  |  2 +-
 benchmark/taylor_jl_01.cpp    |  2 +-
 doc/tut_adaptive.rst          |  2 +-
 include/heyoka/expression.hpp |  2 +-
 src/expression.cpp            |  4 +--
 src/mascon.cpp                | 12 +++----
 src/nbody.cpp                 | 60 +++++++++++++++++------------------
 test/back_and_forth.cpp       |  2 +-
 test/e3bp.cpp                 |  5 ++-
 test/llvm_state.cpp           |  2 +-
 test/one_body.cpp             |  4 +--
 test/taylor_const_sys.cpp     | 34 ++++++++++----------
 test/taylor_no_decomp_sys.cpp | 12 +++----
 tutorial/adaptive_basic.cpp   |  2 +-
 tutorial/adaptive_opt.cpp     |  2 +-
 tutorial/pendulum.cpp         |  2 +-
 17 files changed, 75 insertions(+), 76 deletions(-)

diff --git a/README.md b/README.md
index 94f557a16..e3f4f7722 100644
--- a/README.md
+++ b/README.md
@@ -52,7 +52,7 @@ int main()
     auto ta = taylor_adaptive<double>{// Definition of the ODE system:
                                       // x' = v
                                       // v' = -9.8 * sin(x)
-                                      {prime(x) == v, prime(v) == -9.8 * sin(x)},
+                                      {prime(x) = v, prime(v) = -9.8 * sin(x)},
                                       // Initial conditions
                                       // for x and v.
                                       {0.05, 0.025}};
diff --git a/benchmark/stiff_equation.cpp b/benchmark/stiff_equation.cpp
index 0a2441529..f501fbb85 100644
--- a/benchmark/stiff_equation.cpp
+++ b/benchmark/stiff_equation.cpp
@@ -23,7 +23,7 @@ int main()
     for (unsigned i = 1u; i < 6; ++i) {
         // Van der pool stiff equation
         auto mu = expression{number{std::pow(10, i)}};
-        taylor_adaptive<double> taylor{{prime(y1) == y2, prime(y2) == mu * (1_dbl - y1 * y1) * y2 - y1}, {2., 0.}};
+        taylor_adaptive<double> taylor{{prime(y1) = y2, prime(y2) = mu * (1_dbl - y1 * y1) * y2 - y1}, {2., 0.}};
         std::cout << "Stiffness param " << mu << ": ";
 
         auto start = high_resolution_clock::now();
diff --git a/benchmark/taylor_jl_01.cpp b/benchmark/taylor_jl_01.cpp
index 9aa2b74a9..37d55af59 100644
--- a/benchmark/taylor_jl_01.cpp
+++ b/benchmark/taylor_jl_01.cpp
@@ -32,7 +32,7 @@ int main()
 
     const auto x0 = 3.0;
 
-    taylor_adaptive<double> ta{{prime(x) == x * x}, {x0}, kw::high_accuracy = true, kw::tol = 1e-18};
+    taylor_adaptive<double> ta{{prime(x) = x * x}, {x0}, kw::high_accuracy = true, kw::tol = 1e-18};
 
     const auto final_time = 0.3333333329479479;
 
diff --git a/doc/tut_adaptive.rst b/doc/tut_adaptive.rst
index d485c8a5f..07b5a013d 100644
--- a/doc/tut_adaptive.rst
+++ b/doc/tut_adaptive.rst
@@ -46,7 +46,7 @@ In this case, we use the ``double`` type, meaning that the integration
 will be carried out in double precision.
 
 As (mandatory) construction arguments, we pass in the system of
-differential equations using the syntax ``prime(x) == ...``, and a set
+differential equations using the syntax ``prime(x) = ...``, and a set
 of initial conditions for ``x`` and ``v`` respectively.
 
 Let's try to print to screen the integrator object:
diff --git a/include/heyoka/expression.hpp b/include/heyoka/expression.hpp
index 0fa413df9..5ce31cfcf 100644
--- a/include/heyoka/expression.hpp
+++ b/include/heyoka/expression.hpp
@@ -144,7 +144,7 @@ struct HEYOKA_DLL_PUBLIC prime_wrapper {
     prime_wrapper &operator=(prime_wrapper &&) noexcept;
     ~prime_wrapper();
 
-    std::pair<expression, expression> operator==(expression) const;
+    std::pair<expression, expression> operator=(expression) &&;
 };
 
 } // namespace detail
diff --git a/src/expression.cpp b/src/expression.cpp
index 72bb840e3..deb30acce 100644
--- a/src/expression.cpp
+++ b/src/expression.cpp
@@ -133,9 +133,9 @@ prime_wrapper &prime_wrapper::operator=(prime_wrapper &&) noexcept = default;
 
 prime_wrapper::~prime_wrapper() = default;
 
-std::pair<expression, expression> prime_wrapper::operator==(expression e) const
+std::pair<expression, expression> prime_wrapper::operator=(expression e) &&
 {
-    return std::pair{expression{variable{m_str}}, std::move(e)};
+    return std::pair{expression{variable{std::move(m_str)}}, std::move(e)};
 }
 
 } // namespace detail
diff --git a/src/mascon.cpp b/src/mascon.cpp
index 61dcf4b88..b0d83780f 100644
--- a/src/mascon.cpp
+++ b/src/mascon.cpp
@@ -58,12 +58,12 @@ make_mascon_system_impl(expression Gconst, std::vector<std::vector<expression>>
 
     // Assembling the return vector containing l.h.s. and r.h.s. (note the fundamental use of pairwise_sum for
     // efficiency and to allow compact mode to do his job)
-    retval.push_back(prime(x) == vx);
-    retval.push_back(prime(y) == vy);
-    retval.push_back(prime(z) == vz);
-    retval.push_back(prime(vx) == pairwise_sum(x_acc) - centripetal_x - coriolis_x);
-    retval.push_back(prime(vy) == pairwise_sum(y_acc) - centripetal_y - coriolis_y);
-    retval.push_back(prime(vz) == pairwise_sum(z_acc) - centripetal_z - coriolis_z);
+    retval.push_back(prime(x) = vx);
+    retval.push_back(prime(y) = vy);
+    retval.push_back(prime(z) = vz);
+    retval.push_back(prime(vx) = pairwise_sum(x_acc) - centripetal_x - coriolis_x);
+    retval.push_back(prime(vy) = pairwise_sum(y_acc) - centripetal_y - coriolis_y);
+    retval.push_back(prime(vz) = pairwise_sum(z_acc) - centripetal_z - coriolis_z);
 
     return retval;
 }
diff --git a/src/nbody.cpp b/src/nbody.cpp
index 9b4fe3dcd..2bb312f5d 100644
--- a/src/nbody.cpp
+++ b/src/nbody.cpp
@@ -106,9 +106,9 @@ std::vector<std::pair<expression, expression>> make_nbody_sys_fixed_masses(std::
         // on all particles.
         for (std::uint32_t i = 0; i < n_fc_massive; ++i) {
             // r' = v.
-            retval.push_back(prime(x_vars[i]) == vx_vars[i]);
-            retval.push_back(prime(y_vars[i]) == vy_vars[i]);
-            retval.push_back(prime(z_vars[i]) == vz_vars[i]);
+            retval.push_back(prime(x_vars[i]) = vx_vars[i]);
+            retval.push_back(prime(y_vars[i]) = vy_vars[i]);
+            retval.push_back(prime(z_vars[i]) = vz_vars[i]);
 
             // Compute the total acceleration on body i,
             // accumulating the result in x/y/z_acc. Part of the
@@ -151,29 +151,29 @@ std::vector<std::pair<expression, expression>> make_nbody_sys_fixed_masses(std::
             }
 
             // Add the expressions of the accelerations to the system.
-            retval.push_back(prime(vx_vars[i]) == pairwise_sum(x_acc[i]));
-            retval.push_back(prime(vy_vars[i]) == pairwise_sum(y_acc[i]));
-            retval.push_back(prime(vz_vars[i]) == pairwise_sum(z_acc[i]));
+            retval.push_back(prime(vx_vars[i]) = pairwise_sum(x_acc[i]));
+            retval.push_back(prime(vy_vars[i]) = pairwise_sum(y_acc[i]));
+            retval.push_back(prime(vz_vars[i]) = pairwise_sum(z_acc[i]));
         }
 
         // All the accelerations on the massless particles
         // have already been accumulated in the loop above.
         // We just need to perform the pairwise sums on x/y/z_acc.
         for (std::uint32_t i = n_fc_massive; i < n; ++i) {
-            retval.push_back(prime(x_vars[i]) == vx_vars[i]);
-            retval.push_back(prime(y_vars[i]) == vy_vars[i]);
-            retval.push_back(prime(z_vars[i]) == vz_vars[i]);
+            retval.push_back(prime(x_vars[i]) = vx_vars[i]);
+            retval.push_back(prime(y_vars[i]) = vy_vars[i]);
+            retval.push_back(prime(z_vars[i]) = vz_vars[i]);
 
-            retval.push_back(prime(vx_vars[i]) == pairwise_sum(x_acc[i]));
-            retval.push_back(prime(vy_vars[i]) == pairwise_sum(y_acc[i]));
-            retval.push_back(prime(vz_vars[i]) == pairwise_sum(z_acc[i]));
+            retval.push_back(prime(vx_vars[i]) = pairwise_sum(x_acc[i]));
+            retval.push_back(prime(vy_vars[i]) = pairwise_sum(y_acc[i]));
+            retval.push_back(prime(vz_vars[i]) = pairwise_sum(z_acc[i]));
         }
     } else {
         for (std::uint32_t i = 0; i < n; ++i) {
             // r' = v.
-            retval.push_back(prime(x_vars[i]) == vx_vars[i]);
-            retval.push_back(prime(y_vars[i]) == vy_vars[i]);
-            retval.push_back(prime(z_vars[i]) == vz_vars[i]);
+            retval.push_back(prime(x_vars[i]) = vx_vars[i]);
+            retval.push_back(prime(y_vars[i]) = vy_vars[i]);
+            retval.push_back(prime(z_vars[i]) = vz_vars[i]);
 
             // Compute the total acceleration on body i,
             // accumulating the result in x/y/z_acc. Part of the
@@ -215,9 +215,9 @@ std::vector<std::pair<expression, expression>> make_nbody_sys_fixed_masses(std::
             }
 
             // Add the expressions of the accelerations to the system.
-            retval.push_back(prime(vx_vars[i]) == pairwise_sum(x_acc[i]));
-            retval.push_back(prime(vy_vars[i]) == pairwise_sum(y_acc[i]));
-            retval.push_back(prime(vz_vars[i]) == pairwise_sum(z_acc[i]));
+            retval.push_back(prime(vx_vars[i]) = pairwise_sum(x_acc[i]));
+            retval.push_back(prime(vy_vars[i]) = pairwise_sum(y_acc[i]));
+            retval.push_back(prime(vz_vars[i]) = pairwise_sum(z_acc[i]));
         }
     }
 
@@ -266,9 +266,9 @@ std::vector<std::pair<expression, expression>> make_nbody_sys_par_masses(std::ui
     // on all particles.
     for (std::uint32_t i = 0; i < n_massive; ++i) {
         // r' = v.
-        retval.push_back(prime(x_vars[i]) == vx_vars[i]);
-        retval.push_back(prime(y_vars[i]) == vy_vars[i]);
-        retval.push_back(prime(z_vars[i]) == vz_vars[i]);
+        retval.push_back(prime(x_vars[i]) = vx_vars[i]);
+        retval.push_back(prime(y_vars[i]) = vy_vars[i]);
+        retval.push_back(prime(z_vars[i]) = vz_vars[i]);
 
         // Compute the total acceleration on body i,
         // accumulating the result in x/y/z_acc. Part of the
@@ -311,22 +311,22 @@ std::vector<std::pair<expression, expression>> make_nbody_sys_par_masses(std::ui
         }
 
         // Add the expressions of the accelerations to the system.
-        retval.push_back(prime(vx_vars[i]) == pairwise_sum(x_acc[i]));
-        retval.push_back(prime(vy_vars[i]) == pairwise_sum(y_acc[i]));
-        retval.push_back(prime(vz_vars[i]) == pairwise_sum(z_acc[i]));
+        retval.push_back(prime(vx_vars[i]) = pairwise_sum(x_acc[i]));
+        retval.push_back(prime(vy_vars[i]) = pairwise_sum(y_acc[i]));
+        retval.push_back(prime(vz_vars[i]) = pairwise_sum(z_acc[i]));
     }
 
     // All the accelerations on the massless particles
     // have already been accumulated in the loop above.
     // We just need to perform the pairwise sums on x/y/z_acc.
     for (std::uint32_t i = n_massive; i < n; ++i) {
-        retval.push_back(prime(x_vars[i]) == vx_vars[i]);
-        retval.push_back(prime(y_vars[i]) == vy_vars[i]);
-        retval.push_back(prime(z_vars[i]) == vz_vars[i]);
+        retval.push_back(prime(x_vars[i]) = vx_vars[i]);
+        retval.push_back(prime(y_vars[i]) = vy_vars[i]);
+        retval.push_back(prime(z_vars[i]) = vz_vars[i]);
 
-        retval.push_back(prime(vx_vars[i]) == pairwise_sum(x_acc[i]));
-        retval.push_back(prime(vy_vars[i]) == pairwise_sum(y_acc[i]));
-        retval.push_back(prime(vz_vars[i]) == pairwise_sum(z_acc[i]));
+        retval.push_back(prime(vx_vars[i]) = pairwise_sum(x_acc[i]));
+        retval.push_back(prime(vy_vars[i]) = pairwise_sum(y_acc[i]));
+        retval.push_back(prime(vz_vars[i]) = pairwise_sum(z_acc[i]));
     }
 
     return retval;
diff --git a/test/back_and_forth.cpp b/test/back_and_forth.cpp
index 14a892493..d424c81ec 100644
--- a/test/back_and_forth.cpp
+++ b/test/back_and_forth.cpp
@@ -43,7 +43,7 @@ TEST_CASE("pendulum")
     auto tester = [](auto fp_x, bool high_accuracy, bool compact_mode) {
         using fp_t = decltype(fp_x);
 
-        taylor_adaptive<fp_t> ta{{prime("th"_var) == "v"_var, prime("v"_var) == -9.8_dbl / 1.5_dbl * sin("th"_var)},
+        taylor_adaptive<fp_t> ta{{prime("th"_var) = "v"_var, prime("v"_var) = -9.8_dbl / 1.5_dbl * sin("th"_var)},
                                  {fp_t{0.05}, fp_t{0.025}},
                                  kw::high_accuracy = high_accuracy,
                                  kw::compact_mode = compact_mode};
diff --git a/test/e3bp.cpp b/test/e3bp.cpp
index ff0691f8b..60229c119 100644
--- a/test/e3bp.cpp
+++ b/test/e3bp.cpp
@@ -69,9 +69,8 @@ TEST_CASE("e3bp")
         = cart_to_ham_ec(std::vector{1.20793759666736, -0.493320558636725, 1.19760678594565, -0.498435147674914,
                                      0.548228167205306, 0.49662691628363});
 
-    taylor_adaptive_dbl tad{{prime(xi) == diff(ham, pxi), prime(eta) == diff(ham, peta), prime(phi) == diff(ham, pphi),
-                             prime(pxi) == -diff(ham, xi), prime(peta) == -diff(ham, eta),
-                             prime(pphi) == -diff(ham, phi)},
+    taylor_adaptive_dbl tad{{prime(xi) = diff(ham, pxi), prime(eta) = diff(ham, peta), prime(phi) = diff(ham, pphi),
+                             prime(pxi) = -diff(ham, xi), prime(peta) = -diff(ham, eta), prime(pphi) = -diff(ham, phi)},
                             init_state};
 
     const auto &st = tad.get_state();
diff --git a/test/llvm_state.cpp b/test/llvm_state.cpp
index 8e976cf75..1adb4bffb 100644
--- a/test/llvm_state.cpp
+++ b/test/llvm_state.cpp
@@ -22,7 +22,7 @@ TEST_CASE("basic")
 
     llvm_state s;
     auto [x, y] = make_vars("x", "y");
-    taylor_add_jet_dbl(s, "foo", {prime(x) == y, prime(y) == (1_dbl - x * x) * y - x}, 21, 1, true, false);
+    taylor_add_jet_dbl(s, "foo", {prime(x) = y, prime(y) = (1_dbl - x * x) * y - x}, 21, 1, true, false);
 
     std::cout << s.get_ir() << '\n';
 }
diff --git a/test/one_body.cpp b/test/one_body.cpp
index 216c58559..604ca61b4 100644
--- a/test/one_body.cpp
+++ b/test/one_body.cpp
@@ -29,8 +29,8 @@ TEST_CASE("one_body")
 
     auto [s_x, s_v] = kep_to_cart({1., .5, .0001, 0., 0., 0.}, 1.);
 
-    taylor_adaptive<double> ta{{prime(x) == vx, prime(y) == vy, prime(z) == vz, prime(vx) == -x * r_m3,
-                                prime(vy) == -y * r_m3, prime(vz) == -z * r_m3},
+    taylor_adaptive<double> ta{{prime(x) = vx, prime(y) = vy, prime(z) = vz, prime(vx) = -x * r_m3,
+                                prime(vy) = -y * r_m3, prime(vz) = -z * r_m3},
                                {s_x[0], s_x[1], s_x[2], s_v[0], s_v[1], s_v[2]},
                                kw::high_accuracy = true,
                                kw::tol = 1e-18};
diff --git a/test/taylor_const_sys.cpp b/test/taylor_const_sys.cpp
index e0753781b..01813d0fb 100644
--- a/test/taylor_const_sys.cpp
+++ b/test/taylor_const_sys.cpp
@@ -92,8 +92,8 @@ TEST_CASE("taylor const sys")
             llvm_state s{kw::opt_level = opt_level};
 
             taylor_add_jet<fp_t>(s, "jet",
-                                 {prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
-                                  prime(z) == expression{number{fp_t{0}}}},
+                                 {prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
+                                  prime(z) = expression{number{fp_t{0}}}},
                                  1, 1, high_accuracy, compact_mode);
 
             s.compile();
@@ -118,8 +118,8 @@ TEST_CASE("taylor const sys")
 
             taylor_add_jet<fp_t>(
                 s, "jet",
-                {prime(x) == par[0], prime(y) == expression{number{fp_t{-2}}}, prime(z) == expression{number{fp_t{0}}}},
-                1, 1, high_accuracy, compact_mode);
+                {prime(x) = par[0], prime(y) = expression{number{fp_t{-2}}}, prime(z) = expression{number{fp_t{0}}}}, 1,
+                1, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -144,8 +144,8 @@ TEST_CASE("taylor const sys")
             llvm_state s{kw::opt_level = opt_level};
 
             taylor_add_jet<fp_t>(s, "jet",
-                                 {prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
-                                  prime(z) == expression{number{fp_t{0}}}},
+                                 {prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
+                                  prime(z) = expression{number{fp_t{0}}}},
                                  1, 2, high_accuracy, compact_mode);
 
             s.compile();
@@ -181,8 +181,8 @@ TEST_CASE("taylor const sys")
 
             taylor_add_jet<fp_t>(
                 s, "jet",
-                {prime(x) == expression{number{fp_t{1}}}, prime(y) == par[1], prime(z) == expression{number{fp_t{0}}}},
-                1, 2, high_accuracy, compact_mode);
+                {prime(x) = expression{number{fp_t{1}}}, prime(y) = par[1], prime(z) = expression{number{fp_t{0}}}}, 1,
+                2, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -218,8 +218,8 @@ TEST_CASE("taylor const sys")
             llvm_state s{kw::opt_level = opt_level};
 
             taylor_add_jet<fp_t>(s, "jet",
-                                 {prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
-                                  prime(z) == expression{number{fp_t{0}}}},
+                                 {prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
+                                  prime(z) = expression{number{fp_t{0}}}},
                                  2, 1, high_accuracy, compact_mode);
 
             s.compile();
@@ -246,8 +246,8 @@ TEST_CASE("taylor const sys")
             llvm_state s{kw::opt_level = opt_level};
 
             taylor_add_jet<fp_t>(s, "jet",
-                                 {prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
-                                  prime(z) == expression{number{fp_t{0}}}},
+                                 {prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
+                                  prime(z) = expression{number{fp_t{0}}}},
                                  2, 2, high_accuracy, compact_mode);
 
             s.compile();
@@ -291,8 +291,8 @@ TEST_CASE("taylor const sys")
             llvm_state s{kw::opt_level = opt_level};
 
             taylor_add_jet<fp_t>(s, "jet",
-                                 {prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
-                                  prime(z) == expression{number{fp_t{0}}}},
+                                 {prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
+                                  prime(z) = expression{number{fp_t{0}}}},
                                  3, 3, high_accuracy, compact_mode);
 
             s.compile();
@@ -336,7 +336,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) == par[0], prime(y) == par[1], prime(z) == par[2]}, 3, 3,
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) = par[0], prime(y) = par[1], prime(z) = par[2]}, 3, 3,
                                  high_accuracy, compact_mode);
 
             s.compile();
@@ -380,8 +380,8 @@ TEST_CASE("taylor const sys")
         }
 
         // Do the batch/scalar comparison.
-        compare_batch_scalar<fp_t>({prime(x) == expression{number{fp_t{1}}}, prime(y) == expression{number{fp_t{-2}}},
-                                    prime(z) == expression{number{fp_t{0}}}},
+        compare_batch_scalar<fp_t>({prime(x) = expression{number{fp_t{1}}}, prime(y) = expression{number{fp_t{-2}}},
+                                    prime(z) = expression{number{fp_t{0}}}},
                                    opt_level, high_accuracy, compact_mode);
     };
 
diff --git a/test/taylor_no_decomp_sys.cpp b/test/taylor_no_decomp_sys.cpp
index da36ba774..0d6fae57d 100644
--- a/test/taylor_no_decomp_sys.cpp
+++ b/test/taylor_no_decomp_sys.cpp
@@ -90,7 +90,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) == y, prime(y) == x}, 1, 1, high_accuracy, compact_mode);
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) = y, prime(y) = x}, 1, 1, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -110,7 +110,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) == y, prime(y) == x}, 1, 2, high_accuracy, compact_mode);
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) = y, prime(y) = x}, 1, 2, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -137,7 +137,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) == y, prime(y) == x}, 2, 1, high_accuracy, compact_mode);
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) = y, prime(y) = x}, 2, 1, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -159,7 +159,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) == y, prime(y) == x}, 2, 2, high_accuracy, compact_mode);
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) = y, prime(y) = x}, 2, 2, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -192,7 +192,7 @@ TEST_CASE("taylor const sys")
         {
             llvm_state s{kw::opt_level = opt_level};
 
-            taylor_add_jet<fp_t>(s, "jet", {prime(x) == y, prime(y) == x}, 3, 3, high_accuracy, compact_mode);
+            taylor_add_jet<fp_t>(s, "jet", {prime(x) = y, prime(y) = x}, 3, 3, high_accuracy, compact_mode);
 
             s.compile();
 
@@ -237,7 +237,7 @@ TEST_CASE("taylor const sys")
         }
 
         // Do the batch/scalar comparison.
-        compare_batch_scalar<fp_t>({prime(x) == y, prime(y) == x}, opt_level, high_accuracy, compact_mode);
+        compare_batch_scalar<fp_t>({prime(x) = y, prime(y) = x}, opt_level, high_accuracy, compact_mode);
     };
 
     for (auto cm : {true, false}) {
diff --git a/tutorial/adaptive_basic.cpp b/tutorial/adaptive_basic.cpp
index ce7c09216..2797f0841 100644
--- a/tutorial/adaptive_basic.cpp
+++ b/tutorial/adaptive_basic.cpp
@@ -22,7 +22,7 @@ int main()
     auto ta = taylor_adaptive<double>{// Definition of the ODE system:
                                       // x' = v
                                       // v' = -9.8 * sin(x)
-                                      {prime(x) == v, prime(v) == -9.8 * sin(x)},
+                                      {prime(x) = v, prime(v) = -9.8 * sin(x)},
                                       // Initial conditions
                                       // for x and v.
                                       {0.05, 0.025}};
diff --git a/tutorial/adaptive_opt.cpp b/tutorial/adaptive_opt.cpp
index afcad4516..094bb91a8 100644
--- a/tutorial/adaptive_opt.cpp
+++ b/tutorial/adaptive_opt.cpp
@@ -25,7 +25,7 @@ int main()
     auto ta = taylor_adaptive<double>{// Definition of the ODE system:
                                       // x' = v
                                       // v' = -9.8 * sin(x)
-                                      {prime(x) == v, prime(v) == -9.8 * sin(x)},
+                                      {prime(x) = v, prime(v) = -9.8 * sin(x)},
                                       // Initial conditions
                                       // for x and v.
                                       {0.05, 0.025},
diff --git a/tutorial/pendulum.cpp b/tutorial/pendulum.cpp
index fe973746c..9266a8f3e 100644
--- a/tutorial/pendulum.cpp
+++ b/tutorial/pendulum.cpp
@@ -22,7 +22,7 @@ int main()
     auto ta = taylor_adaptive<double>{// Definition of the ODE system:
                                       // x' = v
                                       // v' = -9.8 * sin(x)
-                                      {prime(x) == v, prime(v) == -9.8 * sin(x)},
+                                      {prime(x) = v, prime(v) = -9.8 * sin(x)},
                                       // Initial conditions
                                       // for x and v.
                                       {0.05, 0.025}};