Many changes.

7 files changed, 176 insertions, 84 deletions

diff --git a/complex_normal.hpp b/complex_normal.hpp
index 8ed75d1..2ffcea1 100644
--- a/complex_normal.hpp
+++ b/complex_normal.hpp
@@ -6,8 +6,8 @@
 template <class T = double> class complex_normal_distribution {
 private:
   std::normal_distribution<T> d;
-  double δx;
-  double δy;
+  T δx;
+  T δy;
   std::complex<T> μ;
   std::complex<T> eθ;
 
diff --git a/dynamics.hpp b/dynamics.hpp
index d421d13..10b1be2 100644
--- a/dynamics.hpp
+++ b/dynamics.hpp
@@ -13,10 +13,10 @@ class gradientDescentStallException: public std::exception {
   }
 };
 
-template <class Scalar, int p>
-std::tuple<double, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double ε, double γ0 = 1, double δγ = 2) {
+template <class Real, class Scalar, int p>
+std::tuple<Real, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, Real ε, Real γ0 = 1, Real δγ = 2) {
   Vector<Scalar> z = z0;
-  double γ = γ0;
+  Real γ = γ0;
 
   auto [W, dW] = WdW(J, z);
 
@@ -42,12 +42,12 @@ std::tuple<double, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& J, c
   return {W, z};
 }
 
-template <class Scalar, int p>
-Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double ε, double δW = 2, double γ0 = 1, double δγ = 2) {
+template <class Real, class Scalar, int p>
+Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, Real ε, Real δW = 2, Real γ0 = 1, Real δγ = 2) {
   Vector<Scalar> z = z0;
   Vector<Scalar> ζ = euclideanToStereographic(z);
 
-  double W;
+  Real W;
   std::tie(W, std::ignore) = WdW(J, z);
 
   Vector<Scalar> dH;
@@ -61,7 +61,7 @@ Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0,
     Vector<Scalar> ζNew = ζ - dζ;
     Vector<Scalar> zNew = stereographicToEuclidean(ζNew);
 
-    double WNew;
+    Real WNew;
     std::tie(WNew, std::ignore) = WdW(J, zNew);
 
     if (WNew < W) { // If Newton's step lowered the objective, accept it!
@@ -90,20 +90,20 @@ Vector<Scalar> randomVector(unsigned N, Distribution d, Generator& r) {
   return z;
 }
 
-template <class Scalar, int p, class Distribution, class Generator>
-std::tuple<double, Vector<Scalar>> metropolis(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0,
-    std::function<double(const Tensor<Scalar, p>&, const Vector<Scalar>&)>& energy,
-    double T, double γ, unsigned N, Distribution d, Generator& r) {
+template <class Real, class Scalar, int p, class Distribution, class Generator>
+std::tuple<Real, Vector<Scalar>> metropolis(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0,
+    std::function<Real(const Tensor<Scalar, p>&, const Vector<Scalar>&)>& energy,
+    Real T, Real γ, unsigned N, Distribution d, Generator& r) {
   Vector<Scalar> z = z0;
 
-  double E = energy(J, z);
+  Real E = energy(J, z);
 
-  std::uniform_real_distribution<double> D(0, 1);
+  std::uniform_real_distribution<Real> D(0, 1);
 
   for (unsigned i = 0; i < N; i++) {
     Vector<Scalar> zNew = normalize(z + γ * randomVector<Scalar>(z.size(), d, r));
 
-    double ENew = energy(J, zNew);
+    Real ENew = energy(J, zNew);
 
     if (E - ENew > T * log(D(r))) {
       z = zNew;
diff --git a/langevin.cpp b/langevin.cpp
index 51cb2a0..8c191f3 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -14,14 +14,14 @@
 
 #define PSPIN_P 3
 const int p = PSPIN_P; // polynomial degree of Hamiltonian
-using Complex = std::complex<double>;
+using Complex = std::complex<Real>;
 using ComplexVector = Vector<Complex>;
 using ComplexMatrix = Matrix<Complex>;
 using ComplexTensor = Tensor<Complex, p>;
 
 using Rng = randutils::random_generator<pcg32>;
 
-std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_map<uint64_t, ComplexVector>& saddles, double δ) {
+std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_map<uint64_t, ComplexVector>& saddles, Real δ) {
   std::list<std::array<ComplexVector, 2>> pairs;
 
   for (auto it1 = saddles.begin(); it1 != saddles.end(); it1++) {
@@ -36,17 +36,17 @@ std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_m
 }
 
 template <class Generator>
-std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(const ComplexTensor& J0, const ComplexVector& z10, const ComplexVector& z20, double ε, double Δ, Generator& r) {
-  double σ = sqrt(factorial(p) / (2.0 * pow(z10.size(), p - 1)));
-  complex_normal_distribution<> dJ(0, σ, 0);
+std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(const ComplexTensor& J0, const ComplexVector& z10, const ComplexVector& z20, Real ε, Real Δ, Generator& r) {
+  Real σ = sqrt(factorial(p) / (2.0 * pow(z10.size(), p - 1)));
+  complex_normal_distribution<Real> dJ(0, σ, 0);
 
   ComplexTensor J = J0;
   Complex H1, H2;
   ComplexVector z1, z2;
   std::tie(H1, std::ignore, std::ignore) = hamGradHess(J, z10);
   std::tie(H2, std::ignore, std::ignore) = hamGradHess(J, z20);
-  double prevdist = abs(imag(H1-H2));
-  double γ = 0.1 * prevdist;
+  Real prevdist = abs(imag(H1-H2));
+  Real γ = 0.1 * prevdist;
 
   std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ =
     [&γ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) {
@@ -62,7 +62,7 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c
       z1 = findSaddle(Jp, z10, Δ);
       z2 = findSaddle(Jp, z20, Δ);
 
-      double dist = (z1 - z2).norm();
+      Real dist = (z1 - z2).norm();
 
       std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jp, z1);
       std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jp, z2);
@@ -87,16 +87,16 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c
 int main(int argc, char* argv[]) {
   // model parameters
   unsigned N = 10; // number of spins
-  double T = 1;    // temperature
-  double Rκ = 1;   // real part of distribution parameter
-  double Iκ = 0;   // imaginary part of distribution parameter
+  Real T = 1;    // temperature
+  Real Rκ = 1;   // real part of distribution parameter
+  Real Iκ = 0;   // imaginary part of distribution parameter
 
   // simulation parameters
-  double ε = 1e-12;
-  double εJ = 1e-5;
-  double δ = 1e-2;   // threshold for determining saddle
-  double Δ = 1e-3;
-  double γ = 1e-2;   // step size
+  Real ε = 1e-12;
+  Real εJ = 1e-5;
+  Real δ = 1e-2;   // threshold for determining saddle
+  Real Δ = 1e-3;
+  Real γ = 1e-2;   // step size
   unsigned t = 1000; // number of Langevin steps
   unsigned M = 100;
   unsigned n = 100;
@@ -140,18 +140,18 @@ int main(int argc, char* argv[]) {
   }
 
   Complex κ(Rκ, Iκ);
-  double σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1)));
+  Real σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1)));
 
   Rng r;
 
-  complex_normal_distribution<> d(0, 1, 0);
+  complex_normal_distribution<Real> d(0, 1, 0);
 
-  ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine());
+  ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine());
   ComplexVector z = normalize(randomVector<Complex>(N, d, r.engine()));
 
-  std::function<double(const ComplexTensor&, const ComplexVector&)> energyNormGrad = []
+  std::function<Real(const ComplexTensor&, const ComplexVector&)> energyNormGrad = []
     (const ComplexTensor& J, const ComplexVector& z) {
-      double W;
+      Real W;
       std::tie(W, std::ignore) = WdW(J, z);
       return W;
     };
@@ -185,18 +185,16 @@ int main(int argc, char* argv[]) {
   ComplexVector z1 = nearbySaddles.front()[0];
   ComplexVector z2 = nearbySaddles.front()[1];
 
-  for (unsigned j = 2; j < 8; j++) {
-    std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, pow(10, -2.0 * j), ε, r);
+  std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-14, ε, r);
 
-    Rope<Complex> stokes(10, z1, z2, J);
+  Rope stokes(10, z1, z2, J);
 
-    for (unsigned i = 0; i < 8; i++) {
-      stokes.relax(J, 1e5, 0.1);
+  for (unsigned i = 0; i < 9; i++) {
+    stokes.relaxDiscreteGradient(J, 1e6, 0.1, 0);
 
-      std::cout << stokes.n() << " " << stokes.cost(J) << " " << stokes.length() << " " << stokes.error(J) << std::endl;
+    std::cout << stokes.n() << " " << stokes.cost(J) << std::endl;
 
-      stokes = stokes.interpolate();
-    }
+    stokes = stokes.interpolate();
   }
 
   return 0;
diff --git a/p-spin.hpp b/p-spin.hpp
index bd3cacc..f1dc07f 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -2,12 +2,13 @@
 
 #include <eigen3/Eigen/Dense>
 
+#include "types.hpp"
 #include "tensor.hpp"
 #include "factorial.hpp"
 
 template <typename Derived>
 Vector<typename Derived::Scalar> normalize(const Eigen::MatrixBase<Derived>& z) {
-  return z * sqrt((double)z.size() / (typename Derived::Scalar)(z.transpose() * z));
+  return z * sqrt((Real)z.size() / (typename Derived::Scalar)(z.transpose() * z));
 }
 
 template <class Scalar, int p>
@@ -16,7 +17,7 @@ std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> hamGradHess(const Tensor<Scal
   Vector<Scalar> Jzz = Jz * z;
   Scalar Jzzz = Jzz.transpose() * z;
 
-  double pBang = factorial(p);
+  Real pBang = factorial(p);
 
   Matrix<Scalar> hessian = ((p - 1) * p / pBang) * Jz;
   Vector<Scalar> gradient = (p / pBang) * Jzz;
@@ -31,18 +32,18 @@ Vector<Scalar> zDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
 }
 
 template <class Scalar, int p>
-std::tuple<double, Vector<Scalar>> WdW(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
+std::tuple<Real, Vector<Scalar>> WdW(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
   Vector<Scalar> dH;
   Matrix<Scalar> ddH;
   std::tie(std::ignore, dH, ddH) = hamGradHess(J, z);
 
   Vector<Scalar> dzdt = zDot(z, dH);
 
-  double a = z.squaredNorm();
+  Real a = z.squaredNorm();
   Scalar A = (Scalar)(z.transpose() * dzdt) / a;
   Scalar B = dH.dot(z) / a;
 
-  double W = dzdt.squaredNorm();
+  Real W = dzdt.squaredNorm();
   Vector<Scalar> dW = ddH * (dzdt - A * z.conjugate())
     + 2 * (conj(A) * B * z).real()
     - conj(B) * dzdt - conj(A) * dH.conjugate();
diff --git a/stokes.hpp b/stokes.hpp
index 6cc5c4a..cbf5ebb 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,4 +1,6 @@
 #include "p-spin.hpp"
+#include "complex_normal.hpp"
+#include "dynamics.hpp"
 #include <iostream>
 
 class ropeRelaxationStallException: public std::exception {
@@ -9,13 +11,13 @@ class ropeRelaxationStallException: public std::exception {
 
 template <class Scalar>
 Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, const Vector<Scalar>& z´´, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH) {
-  double z² = z.squaredNorm();
-  double z´² = z´.squaredNorm();
+  Real z² = z.squaredNorm();
+  Real z´² = z´.squaredNorm();
 
   Vector<Scalar> ż = zDot(z, dH);
-  double ż² = ż.squaredNorm();
+  Real ż² = ż.squaredNorm();
 
-  double Reż·z´ = real(ż.dot(z´));
+  Real Reż·z´ = real(ż.dot(z´));
 
   Matrix<Scalar> dż = (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²;
   Matrix<Scalar> dżc = -ddH + (ddH * z.conjugate()) * z.transpose() / z²
@@ -33,7 +35,7 @@ Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, con
         ) * z.conjugate()
       ) / z²;
 
-  double dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´));
+  Real dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´));
 
   Vector<Scalar> ddtdLdz´ = - (
       (
@@ -72,7 +74,7 @@ class Rope {
       }
 
       for (unsigned i = 0; i < N + 2; i++) {
-        z[i] = normalize(z.front() + (z.back() - z.front()) * ((double)i / (N + 1.0)));
+        z[i] = normalize(z.front() + (z.back() - z.front()) * ((Real)i / (N + 1.0)));
       }
     }
 
@@ -80,8 +82,8 @@ class Rope {
       return z.size() - 2;
     }
 
-    double length() const {
-      double l = 0;
+    Real length() const {
+      Real l = 0;
 
       for (unsigned i = 0; i < z.size() - 1; i++) {
         l += (z[i + 1] - z[i]).norm();
@@ -91,14 +93,14 @@ class Rope {
     }
 
     template <int p>
-    double error(const Tensor<Scalar, p>& J) const {
+    Real error(const Tensor<Scalar, p>& J) const {
       Scalar H0, HN;
       std::tie(H0, std::ignore, std::ignore) = hamGradHess(J, z.front());
       std::tie(HN, std::ignore, std::ignore) = hamGradHess(J, z.back());
 
-      double ImH = imag((H0 + HN) / 2.0);
+      Real ImH = imag((H0 + HN) / 2.0);
 
-      double err = 0;
+      Real err = 0;
 
       for (unsigned i = 1; i < z.size() - 1; i++) {
         Scalar Hi;
@@ -119,7 +121,7 @@ class Rope {
     }
 
     template <int p>
-    std::vector<Vector<Scalar>> δz(const Tensor<Scalar, p>& J) const {
+    std::vector<Vector<Scalar>> generateGradientδz(const Tensor<Scalar, p>& J) const {
       std::vector<Vector<Scalar>> δz(z.size());
 
 #pragma omp parallel for
@@ -137,7 +139,7 @@ class Rope {
       }
 
       // We return a δz with average norm of one.
-      double mag = 0;
+      Real mag = 0;
       for (unsigned i = 1; i < z.size() - 1; i++) {
         mag += δz[i].norm();
       }
@@ -150,20 +152,87 @@ class Rope {
     }
 
     template <int p>
-    double perturb(const Tensor<Scalar, p>& J, double δ0) {
-      std::vector<Vector<Scalar>> Δz = δz(J);
+    std::vector<Vector<Scalar>> generateDiscreteGradientδz(const Tensor<Scalar, p>& J, Real γ) const {
+      std::vector<Vector<Scalar>> δz(z.size());
+
+      std::vector<Vector<Scalar>> ż(z.size());
+      std::vector<Vector<Scalar>> dH(z.size());
+      std::vector<Matrix<Scalar>> ddH(z.size());
+
+#pragma omp parallel for
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        std::tie(std::ignore, dH[i], ddH[i]) = hamGradHess(J, z[i]);
+        ż[i] = zDot(z[i], dH[i]);
+      }
+
+#pragma omp parallel for
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        Real z² = z[i].squaredNorm();
+        Matrix<Scalar> dż = (dH[i].conjugate() - (dH[i].dot(z[i]) / z²) * z[i].conjugate()) * z[i].adjoint() / z²;
+        Matrix<Scalar> dżc = -ddH[i] + (ddH[i] * z[i].conjugate()) * z[i].transpose() / z²
+          + (z[i].dot(dH[i]) / z²) * (
+              Matrix<Scalar>::Identity(ddH[i].rows(), ddH[i].cols()) - z[i].conjugate() * z[i].transpose() / z²
+            );
+
+        Vector<Scalar> dC = -0.5 * (dżc * dz(i) + dż * dz(i).conjugate() - (dżc * ż[i] + dż * ż[i].conjugate()) * real(ż[i].dot(dz(i))) / ż[i].squaredNorm()) / (ż[i].norm() * dz(i).norm());
+
+        if (i > 1) {
+          dC += -0.5 * (ż[i - 1].conjugate() - dz(i - 1).conjugate() * real(ż[i - 1].dot(dz(i - 1))) / dz(i - 1).squaredNorm()) / (ż[i - 1].norm() * dz(i - 1).norm());
+        }
+
+        if (i < z.size() - 2) {
+          dC +=  0.5 * (ż[i + 1].conjugate() - dz(i + 1).conjugate() * real(ż[i + 1].dot(dz(i + 1))) / dz(i + 1).squaredNorm()) / (ż[i + 1].norm() * dz(i + 1).norm());
+        }
+
+        dC += - γ * (z[i - 1] + z[i + 1]).conjugate();
 
+        δz[i] = dC;
+      }
+
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        δz[i] = δz[i].conjugate() - (δz[i].dot(z[i]) / z[i].squaredNorm()) * z[i].conjugate();
+      }
+
+      // We return a δz with average norm of one.
+      Real mag = 0;
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        mag += δz[i].norm();
+      }
+
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        δz[i] /= mag / n();
+      }
+
+      return δz;
+    }
+
+    template<class Gen>
+    std::vector<Vector<Scalar>> generateRandomδz(Gen& r) const {
+      std::vector<Vector<Scalar>> δz(z.size());
+
+      complex_normal_distribution<> d(0, 1, 0);
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        δz[i] = randomVector<Scalar>(z[0].size(), d, r);
+      }
+
+      return δz;
+    }
+
+    template<int p>
+    Real perturb(const Tensor<Scalar, p>& J, Real δ₀, const std::vector<Vector<Scalar>>& δz, Real γ = 0) {
       Rope rNew = *this;
-      double δ = δ0;
+      Real δ = δ₀;
       // We rescale the step size by the distance between points.
-      double Δl = length() / (n() + 1);
+      Real Δl = length() / (n() + 1);
 
       while (true) {
         for (unsigned i = 1; i < z.size() - 1; i++) {
-          rNew.z[i] = normalize(z[i] - (δ * Δl) * Δz[i]);
+          rNew.z[i] = normalize(z[i] - (δ * Δl) * δz[i]);
         }
 
-        if (rNew.cost(J) < cost(J)) {
+        rNew.spread();
+
+        if (rNew.cost(J, γ) < cost(J, γ)) {
           break;
         } else {
           δ /= 2;
@@ -180,15 +249,15 @@ class Rope {
     }
 
     void spread() {
-      double l = length();
+      Real l = length();
 
-      double a = 0;
+      Real a = 0;
       unsigned pos = 0;
 
       std::vector<Vector<Scalar>> zNew = z;
 
       for (unsigned i = 1; i < z.size() - 1; i++) {
-        double b = i * l / (z.size() - 1);
+        Real b = i * l / (z.size() - 1);
 
         while (b > a) {
           pos++;
@@ -204,27 +273,44 @@ class Rope {
     }
 
     template <int p>
-    double step(const Tensor<Scalar, p>& J, double δ0) {
-      double δ = perturb(J, δ0);
-      spread();
-
-      return δ;
+    void relaxGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0) {
+      Real δ = δ0;
+      try {
+        for (unsigned i = 0; i < N; i++) {
+          std::vector<Vector<Scalar>> δz = generateGradientδz(J);
+          δ = 1.1 * perturb(J, δ, δz);
+        }
+      } catch (std::exception& e) {
+      }
     }
 
     template <int p>
-    void relax(const Tensor<Scalar, p>& J, unsigned N, double δ0) {
-      double δ = δ0;
+    void relaxDiscreteGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Real γ) {
+      Real δ = δ0;
       try {
         for (unsigned i = 0; i < N; i++) {
-          δ = 1.1 * step(J, δ);
+          std::vector<Vector<Scalar>> δz = generateDiscreteGradientδz(J, γ);
+          δ = 1.1 * perturb(J, δ, δz, γ);
         }
       } catch (std::exception& e) {
       }
     }
 
+    template <int p, class Gen>
+    void relaxRandom(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Gen& r) {
+      Real δ = δ0;
+        for (unsigned i = 0; i < N; i++) {
+          try {
+            std::vector<Vector<Scalar>> δz = generateRandomδz(r);
+            δ = 1.1 * perturb(J, δ, δz);
+          } catch (std::exception& e) {
+        }
+      }
+    }
+
     template <int p>
-    double cost(const Tensor<Scalar, p>& J) const {
-      double c = 0;
+    Real cost(const Tensor<Scalar, p>& J, Real γ = 0) const {
+      Real c = 0;
 
       for (unsigned i = 1; i < z.size() - 1; i++) {
         Vector<Scalar> dH;
@@ -235,6 +321,10 @@ class Rope {
         c += 1.0 - real(zD.dot(dz(i))) / zD.norm() / dz(i).norm();
       }
 
+      for (unsigned i = 0; i < z.size() - 1; i++) {
+        c += γ * (z[i + 1] - z[i]).squaredNorm();
+      }
+
       return c;
     }
 
diff --git a/stokes_test.cpp b/stokes_test.cpp
index 7d96c8c..a9347a7 100644
--- a/stokes_test.cpp
+++ b/stokes_test.cpp
@@ -4,8 +4,8 @@ using namespace std::complex_literals;
 
 int main() {
 
-  Vector<std::complex<double>> z(2), dz(2), ddz(2), dH(2);
-  Matrix<std::complex<double>> ddH(2, 2);
+  Vector<std::complex<Real>> z(2), dz(2), ddz(2), dH(2);
+  Matrix<std::complex<Real>> ddH(2, 2);
 
   z << -0.75067 - 0.190132 * 1i, -0.625994 + 0.665987 * 1i;
   dz << -0.0636149 - 0.469166 * 1i, 0.820037 - 0.449064 * 1i;
diff --git a/types.hpp b/types.hpp
new file mode 100644
index 0000000..86ccb8f
--- /dev/null
+++ b/types.hpp
@@ -0,0 +1,3 @@
+#pragma once
+
+using Real = double;