diff --git a/Tests/FermionMatrix/test_FermionMatrix_hermitian.cpp b/Tests/FermionMatrix/test_FermionMatrix_hermitian.cpp
index d1ab974c..39fd8895 100644
--- a/Tests/FermionMatrix/test_FermionMatrix_hermitian.cpp
+++ b/Tests/FermionMatrix/test_FermionMatrix_hermitian.cpp
@@ -51,7 +51,7 @@ COMPLEX_NSL_TEST_CASE( "FermionMatrix: MdaggerM hermitian", "[FermionMatrix, Mda
 //Test if MMdagger is hermitian
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_FermionMatrix_MMdagger_hermitian(const NSL::size_t nt, LatticeType & Lattice, const std::string latticeName, const Type & beta) {
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
+    typedef NSL::complex<NSL::RealTypeOf<Type>> ComplexType;
     NSL::size_t nx = Lattice.sites();
 	NSL::Tensor<Type> phi(nt, nx);
     phi.rand();
@@ -90,7 +90,7 @@ void test_FermionMatrix_MMdagger_hermitian(const NSL::size_t nt, LatticeType & L
 //Test if MdaggerM is hermitian
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_FermionMatrix_MdaggerM_hermitian(const NSL::size_t nt, LatticeType & Lattice, const std::string latticeName, const Type & beta) {
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
+    typedef NSL::complex<NSL::RealTypeOf<Type>> ComplexType;
     NSL::size_t nx = Lattice.sites();
 	NSL::Tensor<Type> phi(nt, nx);
     phi.rand();
diff --git a/Tests/FermionMatrix/test_fermionMatrixHubbardDiag.cpp b/Tests/FermionMatrix/test_fermionMatrixHubbardDiag.cpp
index dd0626e6..34a60f03 100644
--- a/Tests/FermionMatrix/test_fermionMatrixHubbardDiag.cpp
+++ b/Tests/FermionMatrix/test_fermionMatrixHubbardDiag.cpp
@@ -128,7 +128,6 @@ COMPLEX_NSL_TEST_CASE( "fermionMatrixHubbardDiag: logDetM_noninteracting", "[fer
 
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardDiag_M_dense(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::Tensor<Type> phi(nt, nx);
@@ -164,7 +163,6 @@ void test_fermionMatrixHubbardDiag_M_dense(const NSL::size_t nt, LatticeType & L
 //Test for the function Mdagger(psi)
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardDiag_Mdagger(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::Tensor<Type> phi(nt, nx);
@@ -206,7 +204,6 @@ void test_fermionMatrixHubbardDiag_Mdagger(const NSL::size_t nt, LatticeType & L
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardDiag_MMdagger(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
 
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::Tensor<Type> phi(nt, nx);
@@ -216,7 +213,7 @@ void test_fermionMatrixHubbardDiag_MMdagger(const NSL::size_t nt, LatticeType &
 
     NSL::FermionMatrix::HubbardDiag M(Lattice,nt,beta);
     M.populate(phi);
-    ComplexType I={0,1};
+    Type I={0,1};
  
     auto direct = M.MMdagger(psi);
     auto indirect = M.M(M.Mdagger(psi));
@@ -233,7 +230,6 @@ void test_fermionMatrixHubbardDiag_MMdagger(const NSL::size_t nt, LatticeType &
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardDiag_MdaggerM(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
 
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     //hardcoding the calculation done in the method MdaggerM of fermionMatrixHubbardDiag class
@@ -244,7 +240,7 @@ void test_fermionMatrixHubbardDiag_MdaggerM(const NSL::size_t nt, LatticeType &
 
     NSL::FermionMatrix::HubbardDiag M(Lattice,nt,beta);
     M.populate(phi,NSL::Hubbard::Particle);
-    ComplexType I={0,1};
+    Type I={0,1};
  
     auto direct = M.MdaggerM(psi);
     auto indirect = M.Mdagger(M.M(psi));
@@ -265,7 +261,6 @@ void test_logDetM_time_shift_invariance(const NSL::size_t nt, LatticeType & Latt
     // We should find that shifting phi in time doesn't change the determinant.
     int slices_to_shift_by=4;
 
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     size_t nx = Lattice.sites();
     NSL::Tensor<Type> phi(nt, nx), phiShift(nt, nx);
     phi.rand();
@@ -301,8 +296,6 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c
 
     // We should find that by shifting any element of phi by 2π the real part of the determinant doesn't change.
 
-    typedef typename NSL::RT_extractor<Type>::value_type RealType;
-    typedef NSL::complex<RealType> ComplexType;
     NSL::size_t nx = Lattice.sites();
     Type delta = beta/nt;
 
@@ -321,7 +314,7 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c
     random.rand();
     NSL::Tensor<Type> orbits = static_cast<NSL::Tensor<int>>(10*random);
 
-    RealType two_pi = 2*std::numbers::pi_v<RealType>;
+    NSL::RealTypeOf<Type> two_pi = 2*std::numbers::pi_v<NSL::RealTypeOf<Type>>;
     phiShift = phi + two_pi * orbits;
 
     NSL::FermionMatrix::HubbardDiag<Type,LatticeType> M     (Lattice,nt     ,beta);
@@ -331,9 +324,9 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c
 
     Type result = M.logDetM();
     Type result_shift = Mshift.logDetM();
-    RealType diff_imag_mod_two_pi = std::remainder(
-            static_cast<RealType>(NSL::imag(result - result_shift)),
-            static_cast<RealType>(two_pi)
+    NSL::RealTypeOf<Type> diff_imag_mod_two_pi = std::remainder(
+            static_cast<NSL::RealTypeOf<Type>>(NSL::imag(result - result_shift)),
+            static_cast<NSL::RealTypeOf<Type>>(two_pi)
             );
 
 
@@ -345,7 +338,7 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c
 
     //comparing only the real parts
     REQUIRE(almost_equal(result_shift.real(),result.real(),std::numeric_limits<Type>::digits10-1));
-    REQUIRE(almost_equal(static_cast<RealType>(0), 
+    REQUIRE(almost_equal(static_cast<NSL::RealTypeOf<Type>>(0), 
         diff_imag_mod_two_pi,
         std::numeric_limits<Type>::digits10-3));
 
@@ -360,7 +353,6 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_logDetM_noninteracting(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
 
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::Tensor<Type> phi(nt, nx), sausage = NSL::Matrix::Identity<Type>(nx);
diff --git a/Tests/FermionMatrix/test_fermionMatrixHubbardExp.cpp b/Tests/FermionMatrix/test_fermionMatrixHubbardExp.cpp
index bcb3a9ac..21863cae 100644
--- a/Tests/FermionMatrix/test_fermionMatrixHubbardExp.cpp
+++ b/Tests/FermionMatrix/test_fermionMatrixHubbardExp.cpp
@@ -142,7 +142,6 @@ COMPLEX_NSL_TEST_CASE( "fermionMatrixHubbardExp: logDetM_uniform_timeslices", "[
 //Test for the function M(psi)
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardExp_M(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
     //hardcoding the calculation done in the method M of fermionMatrixHubbardExp class
     NSL::Tensor<Type> phi(nt, nx);
@@ -153,7 +152,7 @@ void test_fermionMatrixHubbardExp_M(const NSL::size_t nt, LatticeType & Lattice,
     Type delta = beta/nt;
     NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
     M.populate(phi);
-    ComplexType I ={0,1};
+    Type I ={0,1};
 
     // apply kronecker delta
     //NSL::Tensor<Type> psiShift = NSL::LinAlg::shift(psi,1);
@@ -183,7 +182,6 @@ void test_fermionMatrixHubbardExp_M(const NSL::size_t nt, LatticeType & Lattice,
 
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardExp_M_dense(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::Tensor<Type> phi(nt, nx);
@@ -219,19 +217,18 @@ void test_fermionMatrixHubbardExp_M_dense(const NSL::size_t nt, LatticeType & La
 //Test for the function Mdagger(psi)
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardExp_Mdagger(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
     INFO("nt: "+NSL::to_string(nt)+" nx: "+NSL::to_string(nx));
 
-    NSL::Tensor<ComplexType> phi(nt, nx);
-    NSL::Tensor<ComplexType> psi(nt, nx);
+    NSL::Tensor<Type> phi(nt, nx);
+    NSL::Tensor<Type> psi(nt, nx);
     phi.rand();
     psi.rand();
 
     // To simplify this test one can force the field
-    //phi = phi.real() + ComplexType(0,0);  // phi to be real
-    //psi = psi.real() + ComplexType(0,0);  // psi to be real
-    //psi = psi.imag() * ComplexType(0,1);  // psi to be imaginary
+    //phi = phi.real() + Type(0,0);  // phi to be real
+    //psi = psi.real() + Type(0,0);  // psi to be real
+    //psi = psi.imag() * Type(0,1);  // psi to be imaginary
     NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
     M.populate(phi);
 
@@ -239,22 +236,22 @@ void test_fermionMatrixHubbardExp_Mdagger(const NSL::size_t nt, LatticeType & La
     // First let's check Mdagger against the dagger of the dense implementation of M.
     
     // Construct a dense representation of M† from a dense representation of M
-    NSL::Tensor<ComplexType> Mdense_dagger = M.M_dense(nt).transpose(0,2).transpose(1,3).conj();
+    NSL::Tensor<Type> Mdense_dagger = M.M_dense(nt).transpose(0,2).transpose(1,3).conj();
 
     //  We can also apply Mdagger to the identity matrix in order to get a dense Mdagger.
     //  Follow the M_dense implementation:
 
-    NSL::Tensor<ComplexType> dense(nt, nx, nt, nx);
+    NSL::Tensor<Type> dense(nt, nx, nt, nx);
 
     // Construct the identity matrix.
-    NSL::Tensor<ComplexType> identity(nt, nx, nt, nx);
+    NSL::Tensor<Type> identity(nt, nx, nt, nx);
     for(int t = 0; t < nt; t++){
-        identity(t,NSL::Slice(), t, NSL::Slice()) = NSL::Matrix::Identity<ComplexType>(nx);
+        identity(t,NSL::Slice(), t, NSL::Slice()) = NSL::Matrix::Identity<Type>(nx);
     }
 
     // Ensure it's really the identity in the mat-vec sense.
     // Apply the identity to psi via obvious mat-vec
-    NSL::Tensor<ComplexType> Ipsi(nt, nx);
+    NSL::Tensor<Type> Ipsi(nt, nx);
     for(int t=0; t < nt; t++){
         for(int x=0; x < nx; x++){
             for(int i=0; i< nt; i++){
@@ -275,7 +272,7 @@ void test_fermionMatrixHubbardExp_Mdagger(const NSL::size_t nt, LatticeType & La
     }
     
     // So, we can compare M†.I to (M.I)†
-    REQUIRE( almost_equal(Mdense_dagger-dense, ComplexType(0,0)).all() );
+    REQUIRE( almost_equal(Mdense_dagger-dense, Type(0,0)).all() );
     // This REQUIREment looks funny; why not just check that the two tensors are almost_equal directly, as in
     //      REQUIRE( almost_equal(Mdense_dagger, dense).all() );
     // TODO: almost_equal of +0.0000... and -0.0000... evaluates to False and that's extremely misleading
@@ -284,7 +281,7 @@ void test_fermionMatrixHubbardExp_Mdagger(const NSL::size_t nt, LatticeType & La
 
 
     // Finally, compare two ways of computing M†ψ
-    NSL::Tensor<ComplexType>Mdense_dagger_psi(nt, nx), M_dagger_psi(nt, nx);
+    NSL::Tensor<Type>Mdense_dagger_psi(nt, nx), M_dagger_psi(nt, nx);
     // by doing the obvious mat-vec,
     for(int t=0; t < nt; t++){
         for(int x=0; x < nx; x++){
@@ -308,7 +305,6 @@ void test_fermionMatrixHubbardExp_Mdagger(const NSL::size_t nt, LatticeType & La
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardExp_MdaggerM(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
 
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     //hardcoding the calculation done in the method MdaggerM of fermionMatrixHubbardExp class
@@ -319,7 +315,7 @@ void test_fermionMatrixHubbardExp_MdaggerM(const NSL::size_t nt, LatticeType & L
 
     NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
     M.populate(phi);
-    ComplexType I={0,1};
+    Type I={0,1};
  
     auto direct = M.MdaggerM(psi);
     auto indirect = M.Mdagger(M.M(psi));
@@ -335,7 +331,6 @@ void test_fermionMatrixHubbardExp_MdaggerM(const NSL::size_t nt, LatticeType & L
 
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardExp_MMdagger(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::Tensor<Type> phi(nt, nx);
@@ -345,7 +340,7 @@ void test_fermionMatrixHubbardExp_MMdagger(const NSL::size_t nt, LatticeType & L
 
     NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
     M.populate(phi);
-    ComplexType I={0,1};
+    Type I={0,1};
  
     auto direct = M.MMdagger(psi);
     auto indirect = M.M(M.Mdagger(psi));
@@ -361,7 +356,6 @@ void test_fermionMatrixHubbardExp_MMdagger(const NSL::size_t nt, LatticeType & L
 
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardExp_M_batched(const NSL::size_t nt, LatticeType & Lattice, const Type & beta){
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::size_t Nbatch = 10;
@@ -373,7 +367,7 @@ void test_fermionMatrixHubbardExp_M_batched(const NSL::size_t nt, LatticeType &
 
     NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
     M.populate(phi);
-    ComplexType I={0,1};
+    Type I={0,1};
  
     auto direct = M.M(psi);
 
@@ -389,7 +383,6 @@ void test_fermionMatrixHubbardExp_M_batched(const NSL::size_t nt, LatticeType &
 
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardExp_Mdagger_batched(const NSL::size_t nt, LatticeType & Lattice, const Type & beta){
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::size_t Nbatch = 10;
@@ -401,7 +394,7 @@ void test_fermionMatrixHubbardExp_Mdagger_batched(const NSL::size_t nt, LatticeT
 
     NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
     M.populate(phi);
-    ComplexType I={0,1};
+    Type I={0,1};
  
     auto direct = M.Mdagger(psi);
 
@@ -417,7 +410,6 @@ void test_fermionMatrixHubbardExp_Mdagger_batched(const NSL::size_t nt, LatticeT
 
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardExp_MMdagger_batched(const NSL::size_t nt, LatticeType & Lattice, const Type & beta){
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::size_t Nbatch = 10;
@@ -429,7 +421,7 @@ void test_fermionMatrixHubbardExp_MMdagger_batched(const NSL::size_t nt, Lattice
 
     NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
     M.populate(phi);
-    ComplexType I={0,1};
+    Type I={0,1};
  
     auto direct = M.MMdagger(psi);
 
@@ -445,7 +437,6 @@ void test_fermionMatrixHubbardExp_MMdagger_batched(const NSL::size_t nt, Lattice
 
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_fermionMatrixHubbardExp_MdaggerM_batched(const NSL::size_t nt, LatticeType & Lattice, const Type & beta){
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::size_t Nbatch = 10;
@@ -457,7 +448,7 @@ void test_fermionMatrixHubbardExp_MdaggerM_batched(const NSL::size_t nt, Lattice
 
     NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
     M.populate(phi);
-    ComplexType I={0,1};
+    Type I={0,1};
  
     auto direct = M.MdaggerM(psi);
 
@@ -477,7 +468,6 @@ void test_logDetM_time_shift_invariance(const NSL::size_t nt, LatticeType & Latt
     // We should find that shifting phi in time doesn't change the determinant.
     int slices_to_shift_by=4;
 
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     size_t nx = Lattice.sites();
     NSL::Tensor<Type> phi(nt, nx), phiShift(nt, nx);
     phi.rand();
@@ -512,7 +502,6 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c
 
     // We should find that by shifting any element of phi by 2π the determinant doesn't change.
 
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
     Type delta = beta/nt;
 
@@ -552,7 +541,6 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_logDetM_noninteracting(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
 
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::Tensor<Type> phi(nt, nx);
@@ -583,12 +571,11 @@ void test_logDetM_noninteracting(const NSL::size_t nt, LatticeType & Lattice, co
 template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
 void test_logDetM_uniform_timeslices(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
 
-    typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
     NSL::size_t nx = Lattice.sites();
 
     NSL::Tensor<Type> phi(nt, nx);
     Type delta = beta/nt;
-    ComplexType I ={0,1};
+    Type I ={0,1};
 
     // When phi on a given timeslice is the same on every spatial site
     NSL::Tensor<Type> tmp (nt); tmp.rand();
@@ -598,9 +585,9 @@ void test_logDetM_uniform_timeslices(const NSL::size_t nt, LatticeType & Lattice
     // exp(i phi(t)) matrix is proportional to the identity matrix and can be
     // treated like a scalar.
     // When EVERY timeslice is like that we gather all the scalars together
-    NSL::Tensor<ComplexType> sum(1);
+    NSL::Tensor<Type> sum(1);
     sum(0) = I*phi( NSL::Slice(), 0).sum();
-    ComplexType expsum = NSL::LinAlg::exp(sum)(0);
+    Type expsum = NSL::LinAlg::exp(sum)(0);
 
     // to get
     //      logdet M = logdet( 1 + exp(sum(phi(t))) exp(kappa_tilde * Nt))
diff --git a/Tests/LinAlg/Solver/test_BiCGStab.cpp b/Tests/LinAlg/Solver/test_BiCGStab.cpp
index e0706f88..19d5edea 100644
--- a/Tests/LinAlg/Solver/test_BiCGStab.cpp
+++ b/Tests/LinAlg/Solver/test_BiCGStab.cpp
@@ -15,14 +15,14 @@
  * */
 
 template<typename Type>
-void test_BiCGStab_randomMatrix(const typename NSL::RT_extractor<Type>::type eps, NSL::size_t V);
+void test_BiCGStab_randomMatrix(const NSL::RealTypeOf<Type> eps, NSL::size_t V);
 
 // Notice it is a good idea to check that the `FermionMatrix`
 // is a viable Fermion Matrix, i.e. if it dervives from 
 // `NSL::FermionMatrix::FermionMatrix`. 
 // However, this is done savely in the construction of the BiCGStab regardless.
 template<typename Type, class FermionMatrix>
-void test_BiCGStab_fermionMatrix(FermionMatrix & M,const typename NSL::RT_extractor<Type>::type eps, NSL::size_t Nt, NSL::size_t Nx);
+void test_BiCGStab_fermionMatrix(FermionMatrix & M,const NSL::RealTypeOf<Type> eps, NSL::size_t Nt, NSL::size_t Nx);
 
 // =======================================================================
 // Test Cases
@@ -31,13 +31,13 @@ void test_BiCGStab_fermionMatrix(FermionMatrix & M,const typename NSL::RT_extrac
 FLOAT_NSL_TEST_CASE("BiCGStab - Random Matrix", "[BiCGStab,Random Matrix]"){
 
     // for double types we can demand higer precisions and volumes
-    if constexpr( std::is_same_v<double, typename NSL::RT_extractor<TestType>::type>) {
+    if constexpr( std::is_same_v<double, NSL::RealTypeOf<TestType>>) {
         NSL::size_t V_d = GENERATE(1,2,4,8,10,16,20,30,32,100,1000);
-        typename NSL::RT_extractor<TestType>::type eps_d = GENERATE(1e-1,1e-2,1e-3,1e-4,1e-6,1e-8,1e-10);
+        NSL::RealTypeOf<TestType> eps_d = GENERATE(1e-1,1e-2,1e-3,1e-4,1e-6,1e-8,1e-10);
         test_BiCGStab_randomMatrix<TestType>(eps_d, V_d);
     }  else {
         NSL::size_t V = GENERATE(1,2,4,8,10,16,20,30,32);
-        typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-1,1e-2,1e-3,1e-4);
+        NSL::RealTypeOf<TestType> eps = GENERATE(1e-1,1e-2,1e-3,1e-4);
         test_BiCGStab_randomMatrix<TestType>(eps, V);
     }
 }
@@ -47,7 +47,7 @@ COMPLEX_NSL_TEST_CASE("BiCGStab - Hubbard Exp - Complete", "[BiCGStab,Hubbard Ex
     NSL::size_t Nx = GENERATE(2,8);
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -64,7 +64,7 @@ COMPLEX_NSL_TEST_CASE("BiCGStab - Hubbard Exp - Ring", "[BiCGStab,Hubbard Exp,Ri
     NSL::size_t Nx = GENERATE(2,8);
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -80,7 +80,7 @@ COMPLEX_NSL_TEST_CASE("BiCGStab - Hubbard Exp - Square 1D", "[BiCGStab,Hubbard E
     NSL::size_t Nx = GENERATE(2,8);
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -98,7 +98,7 @@ COMPLEX_NSL_TEST_CASE("BiCGStab - Hubbard Exp - Square 2D", "[BiCGStab,Hubbard E
     NSL::size_t Nx = Nx1*Nx2;
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -117,7 +117,7 @@ COMPLEX_NSL_TEST_CASE("BiCGStab - Hubbard Exp - Square 3D", "[BiCGStab,Hubbard E
     NSL::size_t Nx = Nx1*Nx2*Nx3;
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -137,7 +137,7 @@ COMPLEX_NSL_TEST_CASE("BiCGStab - Hubbard Exp - Square 4D", "[BiCGStab,Hubbard E
     NSL::size_t Nx = Nx1*Nx2*Nx3*Nx4;
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -154,7 +154,7 @@ COMPLEX_NSL_TEST_CASE("BiCGStab - Hubbard Exp - Square 4D", "[BiCGStab,Hubbard E
 // ======================================================================
 
 template<typename Type>
-void test_BiCGStab_randomMatrix(const typename NSL::RT_extractor<Type>::type eps, NSL::size_t V){
+void test_BiCGStab_randomMatrix(const NSL::RealTypeOf<Type> eps, NSL::size_t V){
     INFO( std::string("V = ") + std::to_string(V) );
     INFO( std::string("eps = ") + std::to_string(eps) );
     INFO( std::string("Matching Digits = ") + std::to_string(getMatchingDigits(eps)) );
@@ -187,7 +187,7 @@ void test_BiCGStab_randomMatrix(const typename NSL::RT_extractor<Type>::type eps
 // Implementation details: test_BiCGStab_fermionMatrix
 // ======================================================================
 template<typename Type, class FermionMatrix>
-void test_BiCGStab_fermionMatrix(FermionMatrix & M,const typename NSL::RT_extractor<Type>::type eps, NSL::size_t Nt, NSL::size_t Nx){
+void test_BiCGStab_fermionMatrix(FermionMatrix & M,const NSL::RealTypeOf<Type> eps, NSL::size_t Nt, NSL::size_t Nx){
     NSL::Tensor<Type> b(Nt,Nx);b.rand();
 
     // MdaggerM =========================================================
diff --git a/Tests/LinAlg/Solver/test_CG.cpp b/Tests/LinAlg/Solver/test_CG.cpp
index 3a530fec..59f6e718 100644
--- a/Tests/LinAlg/Solver/test_CG.cpp
+++ b/Tests/LinAlg/Solver/test_CG.cpp
@@ -15,14 +15,14 @@
  * */
 
 template<typename Type>
-void test_CG_randomMatrix(const typename NSL::RT_extractor<Type>::type eps, NSL::size_t V);
+void test_CG_randomMatrix(const NSL::RealTypeOf<Type> eps, NSL::size_t V);
 
 // Notice it is a good idea to check that the `FermionMatrix`
 // is a viable Fermion Matrix, i.e. if it dervives from 
 // `NSL::FermionMatrix::FermionMatrix`. 
 // However, this is done savely in the construction of the CG regardless.
 template<typename Type, class FermionMatrix>
-void test_CG_fermionMatrix(FermionMatrix & M,const typename NSL::RT_extractor<Type>::type eps, NSL::size_t Nt, NSL::size_t Nx);
+void test_CG_fermionMatrix(FermionMatrix & M,const NSL::RealTypeOf<Type> eps, NSL::size_t Nt, NSL::size_t Nx);
 
 // =======================================================================
 // Test Cases
@@ -31,13 +31,13 @@ void test_CG_fermionMatrix(FermionMatrix & M,const typename NSL::RT_extractor<Ty
 FLOAT_NSL_TEST_CASE("CG - Random Matrix", "[CG,Random Matrix]"){
 
     // for double types we can demand higer precisions and volumes
-    if constexpr( std::is_same_v<double, typename NSL::RT_extractor<TestType>::type>) {
+    if constexpr( std::is_same_v<double, NSL::RealTypeOf<TestType>>) {
         NSL::size_t V_d = GENERATE(1,2,4,8,10,16,20,30,32,100,1000);
-        typename NSL::RT_extractor<TestType>::type eps_d = GENERATE(1e-1,1e-2,1e-3,1e-4,1e-6,1e-8,1e-10);
+        NSL::RealTypeOf<TestType> eps_d = GENERATE(1e-1,1e-2,1e-3,1e-4,1e-6,1e-8,1e-10);
         test_CG_randomMatrix<TestType>(eps_d, V_d);
     }  else {
         NSL::size_t V = GENERATE(1,2,4,8,10,16,20,30,32);
-        typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-1,1e-2,1e-3,1e-4);
+        NSL::RealTypeOf<TestType> eps = GENERATE(1e-1,1e-2,1e-3,1e-4);
         test_CG_randomMatrix<TestType>(eps, V);
     }
 }
@@ -47,7 +47,7 @@ COMPLEX_NSL_TEST_CASE("CG - Hubbard Exp - Complete", "[CG,Hubbard Exp,Complete]"
     NSL::size_t Nx = GENERATE(2,8);
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -64,7 +64,7 @@ COMPLEX_NSL_TEST_CASE("CG - Hubbard Exp - Ring", "[CG,Hubbard Exp,Ring]"){
     NSL::size_t Nx = GENERATE(2,8);
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -80,7 +80,7 @@ COMPLEX_NSL_TEST_CASE("CG - Hubbard Exp - Square 1D", "[CG,Hubbard Exp,Square,1D
     NSL::size_t Nx = GENERATE(2,8);
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -98,7 +98,7 @@ COMPLEX_NSL_TEST_CASE("CG - Hubbard Exp - Square 2D", "[CG,Hubbard Exp,Square,2D
     NSL::size_t Nx = Nx1*Nx2;
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -117,7 +117,7 @@ COMPLEX_NSL_TEST_CASE("CG - Hubbard Exp - Square 3D", "[CG,Hubbard Exp,Square,3D
     NSL::size_t Nx = Nx1*Nx2*Nx3;
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -137,7 +137,7 @@ COMPLEX_NSL_TEST_CASE("CG - Hubbard Exp - Square 4D", "[CG,Hubbard Exp,Square,4D
     NSL::size_t Nx = Nx1*Nx2*Nx3*Nx4;
     TestType beta = GENERATE(1,1.5);
 
-    typename NSL::RT_extractor<TestType>::type eps = GENERATE(1e-3,1e-5);
+    NSL::RealTypeOf<TestType> eps = GENERATE(1e-3,1e-5);
 
     INFO(std::string("beta= ") + std::to_string(NSL::real(beta)));
 
@@ -154,7 +154,7 @@ COMPLEX_NSL_TEST_CASE("CG - Hubbard Exp - Square 4D", "[CG,Hubbard Exp,Square,4D
 // ======================================================================
 
 template<typename Type>
-void test_CG_randomMatrix(const typename NSL::RT_extractor<Type>::type eps, NSL::size_t V){
+void test_CG_randomMatrix(const NSL::RealTypeOf<Type> eps, NSL::size_t V){
     INFO( std::string("V = ") + std::to_string(V) );
     INFO( std::string("eps = ") + std::to_string(eps) );
     INFO( std::string("Matching Digits = ") + std::to_string(getMatchingDigits(eps)) );
@@ -187,7 +187,7 @@ void test_CG_randomMatrix(const typename NSL::RT_extractor<Type>::type eps, NSL:
 // Implementation details: test_CG_fermionMatrix
 // ======================================================================
 template<typename Type, class FermionMatrix>
-void test_CG_fermionMatrix(FermionMatrix & M,const typename NSL::RT_extractor<Type>::type eps, NSL::size_t Nt, NSL::size_t Nx){
+void test_CG_fermionMatrix(FermionMatrix & M,const NSL::RealTypeOf<Type> eps, NSL::size_t Nt, NSL::size_t Nx){
     NSL::Tensor<Type> b(Nt,Nx);b.rand();
 
     // MdaggerM =========================================================
diff --git a/Tests/LinAlg/test_LinAlg_mat_exp.cpp b/Tests/LinAlg/test_LinAlg_mat_exp.cpp
index a1e10392..eb93cc78 100644
--- a/Tests/LinAlg/test_LinAlg_mat_exp.cpp
+++ b/Tests/LinAlg/test_LinAlg_mat_exp.cpp
@@ -60,7 +60,7 @@ void test_exponential_of_diagonal(const size_type & size){
     }
 }
 
-template<typename T, typename RT = typename NSL::RT_extractor<T>::value_type>
+template<typename T, typename RT = NSL::RealTypeOf<T>>
 void test_exponential_of_hermitian(const size_type & size){
     INFO("Type = " << typeid(T).name());
     INFO("size = " << size);
diff --git a/Tests/Tensor/test_Tensor_trig.cpp b/Tests/Tensor/test_Tensor_trig.cpp
index bccd3149..66d7e1b6 100644
--- a/Tests/Tensor/test_Tensor_trig.cpp
+++ b/Tests/Tensor/test_Tensor_trig.cpp
@@ -153,21 +153,21 @@ void test_half_periods(SizeTypes ... Ns){
     // sine
     {
         NSL::Tensor<Type> unshifted(A,true);
-        NSL::Tensor<Type> shifted(A+std::numbers::pi_v<typename NSL::RT_extractor<Type>::type>,true);
+        NSL::Tensor<Type> shifted(A+std::numbers::pi_v<NSL::RealTypeOf<Type>>,true);
         REQUIRE( almost_equal(unshifted.sin() + shifted.sin(), Type(0) ).all() );
     }
 
     // cosine
     {
         NSL::Tensor<Type> unshifted(A,true);
-        NSL::Tensor<Type> shifted(A+std::numbers::pi_v<typename NSL::RT_extractor<Type>::type>,true);
+        NSL::Tensor<Type> shifted(A+std::numbers::pi_v<NSL::RealTypeOf<Type>>,true);
         REQUIRE( almost_equal(unshifted.cos() + shifted.cos(), Type(0) ).all() );
     }
 
     // tan
     {
         NSL::Tensor<Type> unshifted(A,true);
-        NSL::Tensor<Type> shifted(A+std::numbers::pi_v<typename NSL::RT_extractor<Type>::type>,true);
+        NSL::Tensor<Type> shifted(A+std::numbers::pi_v<NSL::RealTypeOf<Type>>,true);
         REQUIRE( almost_equal(unshifted.tan(), shifted.tan()).all() );
     }
 
diff --git a/Tests/test.hpp b/Tests/test.hpp
index 783fbe28..cea66c07 100644
--- a/Tests/test.hpp
+++ b/Tests/test.hpp
@@ -63,8 +63,8 @@
  *
  * */
 template<NSL::Concept::isFloatingPoint T>
-bool compare_floating_point(T a, T b, typename NSL::RT_extractor<T>::type factor = 10){
-    return std::abs(static_cast<T>(1) - a/b) <= factor*std::numeric_limits<typename NSL::RT_extractor<T>::type>::epsilon();
+bool compare_floating_point(T a, T b, NSL::RealTypeOf<T> factor = 10){
+    return std::abs(static_cast<T>(1) - a/b) <= factor*std::numeric_limits<NSL::RealTypeOf<T>>::epsilon();
 }
 
 //! compare two integer numbers a,b
@@ -132,7 +132,7 @@ template<typename Type>
 NSL::Tensor<bool> almost_equal(NSL::Tensor<Type> x, NSL::Tensor<Type> y, int matchingDigits = std::numeric_limits<Type>::digits10){
     assertm( y.shape() == x.shape(), "To be almost equal two tensors must be the same shape.");
 
-    NSL::Tensor<bool> result(static_cast<NSL::Tensor<typename NSL::RT_extractor<Type>::value_type>>(x));
+    NSL::Tensor<bool> result(static_cast<NSL::Tensor<NSL::RealTypeOf<Type>>>(x));
     result = false;
     NSL::size_t elements = x.numel();
     for(NSL::size_t i = 0; i < elements; i++){
@@ -144,7 +144,7 @@ NSL::Tensor<bool> almost_equal(NSL::Tensor<Type> x, NSL::Tensor<Type> y, int mat
 
 template<typename Type>
 NSL::Tensor<bool> almost_equal(NSL::Tensor<Type> x, Type y, int matchingDigits = std::numeric_limits<Type>::digits10){
-    NSL::Tensor<bool> result(static_cast<NSL::Tensor<typename NSL::RT_extractor<Type>::value_type>>(x));
+    NSL::Tensor<bool> result(static_cast<NSL::Tensor<NSL::RealTypeOf<Type>>>(x));
     result = false;
     NSL::size_t elements = x.numel();
     for(NSL::size_t i = 0; i < elements; i++){
diff --git a/src/NSL/LinAlg/Solver/Impl/BiCGStab.hpp b/src/NSL/LinAlg/Solver/Impl/BiCGStab.hpp
index b9ebdcec..b057f401 100644
--- a/src/NSL/LinAlg/Solver/Impl/BiCGStab.hpp
+++ b/src/NSL/LinAlg/Solver/Impl/BiCGStab.hpp
@@ -29,7 +29,7 @@ class BiCGStab: public NSL::LinAlg::Solver<Type> {
          * This Solver implementation uses the conjugate gradient (BiCGStab) algorithm.
          * */
         BiCGStab(std::function<NSL::Tensor<Type>(const NSL::Tensor<Type> &)> M,
-               const typename NSL::RT_extractor<Type>::type eps = 1e-6, const NSL::size_t maxIter = 10000) : 
+               const NSL::RealTypeOf<Type> eps = 1e-6, const NSL::size_t maxIter = 10000) : 
             NSL::LinAlg::Solver<Type>(M),
             errSq_(eps*eps),
             maxIter_(maxIter),
@@ -85,7 +85,7 @@ class BiCGStab: public NSL::LinAlg::Solver<Type> {
             // to ensure that the required interface is given.
             requires( NSL::Concept::isDerived<FermionMatrix<Type,LatticeType>,NSL::FermionMatrix::FermionMatrix<Type,LatticeType>> )
         BiCGStab(FermionMatrix<Type,LatticeType> & M,
-               const typename NSL::RT_extractor<Type>::type eps = 1e-6, const NSL::size_t maxIter = 10000) : 
+               const NSL::RealTypeOf<Type> eps = 1e-6, const NSL::size_t maxIter = 10000) : 
             NSL::LinAlg::Solver<Type>(M, NSL::FermionMatrix::M),
             errSq_(eps*eps),
             maxIter_(maxIter),
@@ -150,7 +150,7 @@ class BiCGStab: public NSL::LinAlg::Solver<Type> {
             requires( NSL::Concept::isDerived<FermionMatrix<Type,LatticeType>,NSL::FermionMatrix::FermionMatrix<Type,LatticeType>> )
         BiCGStab(FermionMatrix<Type,LatticeType> & M, 
                NSL::FermionMatrix::MatrixCombination matrixCombination,
-               const typename NSL::RT_extractor<Type>::type eps = 1e-6, const NSL::size_t maxIter = 10000) : 
+               const NSL::RealTypeOf<Type> eps = 1e-6, const NSL::size_t maxIter = 10000) : 
             NSL::LinAlg::Solver<Type>(M,matrixCombination),
             errSq_(eps*eps),
             maxIter_(maxIter),
@@ -184,7 +184,7 @@ class BiCGStab: public NSL::LinAlg::Solver<Type> {
     private:
 
         // precision at which the algorithm is stopped
-        const typename NSL::RT_extractor<Type>::type errSq_;
+        const NSL::RealTypeOf<Type> errSq_;
         // maximum of iterations as fall back in case we don't converge
         const NSL::size_t maxIter_;
 
diff --git a/src/NSL/LinAlg/Solver/Impl/BiCGStab.tpp b/src/NSL/LinAlg/Solver/Impl/BiCGStab.tpp
index 93131e78..f33c7ff6 100644
--- a/src/NSL/LinAlg/Solver/Impl/BiCGStab.tpp
+++ b/src/NSL/LinAlg/Solver/Impl/BiCGStab.tpp
@@ -35,8 +35,8 @@ NSL::Tensor<Type> BiCGStab<Type>::operator()(const NSL::Tensor<Type> & b ){
     // error (this is a simple efficiency optimization)
     // inner_product returns a number of type `Type` from which the real 
     // part is extracted, the imaginary part is 0 by construction
-    typename NSL::RT_extractor<Type>::type rsqr_curr = NSL::real( NSL::LinAlg::inner_product(r_,r_) );
-    typename NSL::RT_extractor<Type>::type rsqr_prev = rsqr_curr;
+    NSL::RealTypeOf<Type> rsqr_curr = NSL::real( NSL::LinAlg::inner_product(r_,r_) );
+    NSL::RealTypeOf<Type> rsqr_prev = rsqr_curr;
     
     // if the guess is already good enough return
     if (rsqr_curr <= errSq_) {
@@ -118,7 +118,7 @@ NSL::Tensor<Type> BiCGStab<Type>::operator()(const NSL::Tensor<Type> & b ){
 
         // compute the momentum update scale
         // beta{i} = (r{i+1},r{i+1})/(r{i},r{i}
-        typename NSL::RT_extractor<Type>::type beta = rsqr_curr / rsqr_prev;
+        NSL::RealTypeOf<Type> beta = rsqr_curr / rsqr_prev;
 
         // now prepare the previous residual square for the next iteration
         rsqr_prev = rsqr_curr;
diff --git a/src/NSL/LinAlg/Solver/Impl/CG.hpp b/src/NSL/LinAlg/Solver/Impl/CG.hpp
index b4f9039d..f8439bf0 100644
--- a/src/NSL/LinAlg/Solver/Impl/CG.hpp
+++ b/src/NSL/LinAlg/Solver/Impl/CG.hpp
@@ -29,7 +29,7 @@ class CG: public NSL::LinAlg::Solver<Type> {
          * This Solver implementation uses the conjugate gradient (CG) algorithm.
          * */
         CG(std::function<NSL::Tensor<Type>(const NSL::Tensor<Type> &)> M,
-               const typename NSL::RT_extractor<Type>::type eps = 1e-12, const NSL::size_t maxIter = 10000) : 
+               const NSL::RealTypeOf<Type> eps = 1e-12, const NSL::size_t maxIter = 10000) : 
             NSL::LinAlg::Solver<Type>(M),
             errSq_(eps*eps),
             maxIter_(maxIter),
@@ -81,7 +81,7 @@ class CG: public NSL::LinAlg::Solver<Type> {
             // to ensure that the required interface is given.
             requires( NSL::Concept::isDerived<FermionMatrix<Type,LatticeType>,NSL::FermionMatrix::FermionMatrix<Type,LatticeType>> )
         CG(FermionMatrix<Type,LatticeType> & M,
-               const typename NSL::RT_extractor<Type>::type eps = 1e-12, const NSL::size_t maxIter = 10000) : 
+               const NSL::RealTypeOf<Type> eps = 1e-12, const NSL::size_t maxIter = 10000) : 
             NSL::LinAlg::Solver<Type>(M, NSL::FermionMatrix::M),
             errSq_(eps*eps),
             maxIter_(maxIter),
@@ -141,7 +141,7 @@ class CG: public NSL::LinAlg::Solver<Type> {
             requires( NSL::Concept::isDerived<FermionMatrix<Type,LatticeType>,NSL::FermionMatrix::FermionMatrix<Type,LatticeType>> )
         CG(FermionMatrix<Type,LatticeType> & M, 
                NSL::FermionMatrix::MatrixCombination matrixCombination,
-               const typename NSL::RT_extractor<Type>::type eps = 1e-12, const NSL::size_t maxIter = 10000) : 
+               const NSL::RealTypeOf<Type> eps = 1e-12, const NSL::size_t maxIter = 10000) : 
             NSL::LinAlg::Solver<Type>(M,matrixCombination),
             errSq_(eps*eps),
             maxIter_(maxIter),
@@ -170,7 +170,7 @@ class CG: public NSL::LinAlg::Solver<Type> {
     private:
 
         // precision at which the algorithm is stopped
-        const typename NSL::RT_extractor<Type>::type errSq_;
+        const NSL::RealTypeOf<Type> errSq_;
         // maximum of iterations as fall back in case we don't converge
         const NSL::size_t maxIter_;
 
diff --git a/src/NSL/LinAlg/Solver/Impl/CG.tpp b/src/NSL/LinAlg/Solver/Impl/CG.tpp
index 3cd1d6f1..2da07476 100644
--- a/src/NSL/LinAlg/Solver/Impl/CG.tpp
+++ b/src/NSL/LinAlg/Solver/Impl/CG.tpp
@@ -35,8 +35,8 @@ NSL::Tensor<Type> CG<Type>::operator()(const NSL::Tensor<Type> & b ){
     // error (this is a simple efficiency optimization)
     // inner_product returns a number of type `Type` from which the real 
     // part is extracted, the imaginary part is 0 by construction
-    typename NSL::RT_extractor<Type>::type rsqr_curr = NSL::real( NSL::LinAlg::inner_product(r_,r_) );
-    typename NSL::RT_extractor<Type>::type rsqr_prev = rsqr_curr;
+    NSL::RealTypeOf<Type> rsqr_curr = NSL::real( NSL::LinAlg::inner_product(r_,r_) );
+    NSL::RealTypeOf<Type> rsqr_prev = rsqr_curr;
     
     // if the guess is already good enough return
     if (rsqr_curr <= errSq_) {
@@ -79,7 +79,7 @@ NSL::Tensor<Type> CG<Type>::operator()(const NSL::Tensor<Type> & b ){
 
         // compute the momentum update scale
         // beta{i} = (r{i+1},r{i+1})/(r{i},r{i}
-        typename NSL::RT_extractor<Type>::type beta = rsqr_curr / rsqr_prev;
+        NSL::RealTypeOf<Type> beta = rsqr_curr / rsqr_prev;
 
         // update the momentum
         // p{i+1} = r{i+1} + beta{i} * p{i}
diff --git a/src/NSL/LinAlg/abs.tpp b/src/NSL/LinAlg/abs.tpp
index fe342bea..73d66811 100644
--- a/src/NSL/LinAlg/abs.tpp
+++ b/src/NSL/LinAlg/abs.tpp
@@ -16,7 +16,7 @@ namespace NSL::LinAlg {
 
     //! Returns the real-type absolute value, regardless of whether the passed value is real or `complex<>`.
 template<NSL::Concept::isNumber Type>
-typename NSL::RT_extractor<Type>::type abs(const Type &value){
+NSL::RealTypeOf<Type> abs(const Type &value){
     if constexpr(is_complex<Type>()) {
         // See NOTE above for std::explanation.
         return std::abs(value);
@@ -27,7 +27,7 @@ typename NSL::RT_extractor<Type>::type abs(const Type &value){
 }
 
 template<typename Type>
-inline NSL::Tensor<typename RT_extractor<Type>::value_type> abs(const NSL::Tensor<Type> &T){
+inline NSL::Tensor<NSL::RealTypeOf<Type>> abs(const NSL::Tensor<Type> &T){
     // preform a deep copy of the tensor;
     NSL::Tensor<Type> Tcopy(T,true);
     return Tcopy.abs();
diff --git a/src/NSL/LinAlg/complex.tpp b/src/NSL/LinAlg/complex.tpp
index 1994a401..3d35fbe4 100644
--- a/src/NSL/LinAlg/complex.tpp
+++ b/src/NSL/LinAlg/complex.tpp
@@ -13,18 +13,18 @@ namespace NSL::LinAlg {
 
 //! Returns the complex conjugate, maintaining type (`complex<>` if `complex<>`, not if not).
 template<NSL::Concept::isNumber Type>
-inline typename NSL::RT_extractor<Type>::type arg(const Type &value){
+inline NSL::RealTypeOf<Type> arg(const Type &value){
     if constexpr(is_complex<Type>()) {
         // See NOTE above for std::explanation.
         return std::arg(value);
     }
     else {
-        if(value > 0) return static_cast<typename NSL::RT_extractor<Type>::type>(0);
+        if(value > 0) return static_cast<NSL::RealTypeOf<Type>>(0);
         //! todo We should be very careful about the branch-cut of arg.
         // If we want arg to be single-valued, we should pick a finite interval,
         // say (-π,+π], which is 2π periodic.  However, the negative real axis
         // is right on the boundary.
-        return static_cast<typename NSL::RT_extractor<Type>::type>(+std::numbers::pi);
+        return static_cast<NSL::RealTypeOf<Type>>(+std::numbers::pi);
         // I picked + because in Mathematica
         //      Arg[-1] == +π
         // and I trust Wolfram to have these conventions sorted out.
diff --git a/src/NSL/Tensor/Impl/abs.tpp b/src/NSL/Tensor/Impl/abs.tpp
index 5ff22d48..bb8ec156 100644
--- a/src/NSL/Tensor/Impl/abs.tpp
+++ b/src/NSL/Tensor/Impl/abs.tpp
@@ -16,7 +16,7 @@ class TensorAbs:
      * Whether a complex type or a real type, the absolute value is real.
      *
      * */
-    NSL::Tensor<typename RT_extractor<Type>::type> abs(){
+    NSL::Tensor<NSL::RealTypeOf<Type>> abs(){
         return torch::abs(this->data_);
     }
 
diff --git a/src/NSL/Tensor/Impl/realImag.tpp b/src/NSL/Tensor/Impl/realImag.tpp
index bea9f578..57cefa49 100644
--- a/src/NSL/Tensor/Impl/realImag.tpp
+++ b/src/NSL/Tensor/Impl/realImag.tpp
@@ -17,9 +17,9 @@ class TensorReal:
      * Else `Type` refers to a RealType expression (e.g. `float`,`double`, ...)
      * and is simply returned.
      * */
-    NSL::Tensor<typename RT_extractor<Type>::type> real(){
+    NSL::Tensor<NSL::RealTypeOf<Type>> real(){
         if constexpr(NSL::is_complex<Type>()){
-            return NSL::Tensor<typename RT_extractor<Type>::type>(torch::real(this->data_));
+            return NSL::Tensor<NSL::RealTypeOf<Type>>(torch::real(this->data_));
         } else {
             return NSL::Tensor<Type>(this);
         }
@@ -37,11 +37,11 @@ class TensorImag:
      * Else `Type` refers to a RealType expression and does not have an imaginary
      * part, a Tensor with zeros is returned.
      * */
-        Tensor<typename RT_extractor<Type>::type> imag(){
+        Tensor<NSL::RealTypeOf<Type>> imag(){
             if constexpr(NSL::is_complex<Type>()){
-                return Tensor<typename RT_extractor<Type>::type>(torch::imag(this->data_));
+                return Tensor<NSL::RealTypeOf<Type>>(torch::imag(this->data_));
             } else {
-                return Tensor<typename RT_extractor<Type>::type>(torch::zeros_like(this->data_));
+                return Tensor<NSL::RealTypeOf<Type>>(torch::zeros_like(this->data_));
             }
         }
 };
diff --git a/src/NSL/Tensor/tensor.hpp b/src/NSL/Tensor/tensor.hpp
index d5461025..75f824b9 100644
--- a/src/NSL/Tensor/tensor.hpp
+++ b/src/NSL/Tensor/tensor.hpp
@@ -3,7 +3,7 @@
 
 //! \file tensor.hpp
 
-#include "../complex.hpp" // get NSL::RT_extractor
+#include "../complex.hpp" // get NSL::RealTypeOf
 #include "../concepts.hpp" // get NSL::Concept::
 
 
diff --git a/src/NSL/realImag.tpp b/src/NSL/realImag.tpp
index 9cc9ddde..5ee4115a 100644
--- a/src/NSL/realImag.tpp
+++ b/src/NSL/realImag.tpp
@@ -15,7 +15,7 @@ namespace NSL{
 
 //! If `complex<>`, returns the real part; otherwise returns the passed value.
 template<typename Type>
-typename NSL::RT_extractor<Type>::type real(const Type &value){
+NSL::RealTypeOf<Type> real(const Type &value){
     if constexpr(is_complex<Type>()) {
         // See NOTE above for std::explanation.
         return std::real(value);
@@ -27,7 +27,7 @@ typename NSL::RT_extractor<Type>::type real(const Type &value){
 
 //! If `complex<>`, returns the imaginary part; otherwise returns 0.
 template<typename Type>
-typename NSL::RT_extractor<Type>::type imag(const Type &value){
+NSL::RealTypeOf<Type> imag(const Type &value){
     if constexpr(is_complex<Type>()) {
         // See NOTE above for std::explanation.
         return std::imag(value);