Slightly better treatement of the global constant.

1 files changed, 9 insertions, 7 deletions

diff --git a/langevin.cpp b/langevin.cpp
index d704266..119a64b 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -13,13 +13,13 @@
 #include "complex_normal.hpp"
 #include "tensor.hpp"
 
-#define P_SPIN_P 3
-const unsigned p = P_SPIN_P;
+#define PSPIN_P 3
+const unsigned p = PSPIN_P; // polynomial degree of Hamiltonian
 
 using Scalar = std::complex<double>;
 using Vector = Eigen::VectorXcd;
 using Matrix = Eigen::MatrixXcd;
-using Tensor = Eigen::Tensor<Scalar, p>;
+using Tensor = Eigen::Tensor<Scalar, PSPIN_P>;
 
 using Rng = randutils::random_generator<pcg32>;
 
@@ -40,9 +40,11 @@ std::tuple<Scalar, Vector, Matrix> hamGradHess(const Tensor& J, const Vector& z)
   Matrix Jz = contractDown(J, z);
   Vector Jzz = Jz * z;
 
-  Matrix hessian = ((p - 1) * (double)p / factorial(p)) * Jz;
-  Vector gradient = ((double)p / factorial(p)) * Jzz;
-  Scalar hamiltonian = (1.0 / factorial(p)) * Jzz.dot(z);
+  double f = factorial(p);
+
+  Matrix hessian = ((p - 1) * p / f) * Jz;
+  Vector gradient = (p / f) * Jzz;
+  Scalar hamiltonian = (1 / f) * Jzz.dot(z);
 
   return {hamiltonian, gradient, hessian};
 }
@@ -142,7 +144,7 @@ int main(int argc, char* argv[]) {
   Rng r;
   complex_normal_distribution<> d(0, σ, κ);
 
-  Tensor J = generateCouplings<Scalar, p>(N, d, r.engine());
+  Tensor J = generateCouplings<Scalar, PSPIN_P>(N, d, r.engine());
   Vector z = initializeVector(N, r);
 
   std::function<bool(double, unsigned)> findSaddle = [δ](double W, unsigned) {