diff --git a/CMakeLists.txt b/CMakeLists.txt
index e9bbc4e..f98ca48 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -73,7 +73,7 @@ add_library(${MY_TARGET} STATIC
 	include/qboot/context.hpp
 	include/qboot/hor_formula.hpp
 	include/qboot/hor_recursion.hpp
-	include/qboot/my_filesystem.hpp
+	include/qboot/my_concepts.hpp
 	include/qboot/polynomial_program.hpp
 	include/qboot/primary_op.hpp
 	include/qboot/qboot.hpp
@@ -93,6 +93,9 @@ target_link_libraries(
 )
 if(NOT MSVC)
 	target_link_libraries(${MY_TARGET} INTERFACE stdc++fs)
+	if(NOT CMAKE_COMPILER_IS_GNUCXX)
+		target_link_libraries(${MY_TARGET} INTERFACE c++abi)
+	endif()
 endif()
 
 target_include_directories(${MY_TARGET} PUBLIC
diff --git a/README.md b/README.md
index 0da18fc..789d657 100644
--- a/README.md
+++ b/README.md
@@ -28,7 +28,7 @@ Some codes are taken from [cboot](https://github.com/tohtsky/cboot.git).
 
 - [mpfr](http://mpfr.org/)
 
-- [gcc](http://gcc.gnu.org/) (`7.4.0+`) or [clang](http://clang.llvm.org/) (`8.0.0+`)
+- [gcc](http://gcc.gnu.org/) (`10.2.0+`) or [clang](http://clang.llvm.org/) (`11.0.1+`)
 
 ### Windows (MSVC)
 
@@ -36,7 +36,7 @@ Some codes are taken from [cboot](https://github.com/tohtsky/cboot.git).
 
 - [mpfr](https://github.com/BrianGladman/mpfr.git)
 
-- [Visual Studio](https://visualstudio.microsoft.com/) (`2017+`)
+- [Visual Studio](https://visualstudio.microsoft.com/) (`2019+`)
 
 Please build `mpir` and `mpfr` in both `Debug` and `Release` mode.
 
diff --git a/cmake/utils.cmake b/cmake/utils.cmake
index 4052e09..ca33883 100644
--- a/cmake/utils.cmake
+++ b/cmake/utils.cmake
@@ -21,10 +21,11 @@ macro(set_flags)
 		if(CMAKE_COMPILER_IS_GNUCXX)
 			set(cxx_flags ${cxx_flags} -Weffc++ -Wno-noexcept -Walloca -Wcast-align -Wcast-qual -Wconditionally-supported -Wconversion -Wctor-dtor-privacy -Wdate-time -Wdisabled-optimization -Wdouble-promotion -Wduplicated-cond -Wfloat-conversion -Wfloat-equal -Wlogical-op -Wmissing-declarations -Wmissing-include-dirs -Wmultichar -Wmultiple-inheritance -Wnull-dereference -Wold-style-cast -Woverlength-strings -Woverloaded-virtual -Wpacked -Wredundant-decls -Wrestrict -Wsign-conversion -Wsign-promo -Wstrict-null-sentinel -Wsuggest-override -Wswitch-default -Wswitch-enum -Wuninitialized -Wunused-macros -Wunused-parameter -Wuseless-cast -Wvariadic-macros -Wvector-operation-performance -Wvirtual-inheritance)
 		else()
+			set(cxx_flags ${cxx_flags} -stdlib=libc++)
 			if (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 9.0.0)
 				set(cxx_flags ${cxx_flags} -Wno-ctad-maybe-unsupported)
 			endif()
-			set(cxx_flags ${cxx_flags} -Weverything -Wno-c++98-compat -Wno-shadow-field-in-constructor -Wno-c++98-compat-pedantic -Wno-global-constructors -Wno-exit-time-destructors -Wno-covered-switch-default)
+			set(cxx_flags ${cxx_flags} -Weverything -Wno-zero-as-null-pointer-constant -Wno-c++98-compat -Wno-shadow-field-in-constructor -Wno-c++98-compat-pedantic -Wno-global-constructors -Wno-exit-time-destructors -Wno-covered-switch-default)
 		endif()
 	endif()
 endmacro()
@@ -43,7 +44,7 @@ function(set_default_property target dir)
 		DEBUG_POSTFIX "-debug"
 		COMPILE_OPTIONS "${cxx_flags}"
 		CXX_EXTENSIONS OFF
-		CXX_STANDARD 17
+		CXX_STANDARD 20
 		CXX_STANDARD_REQUIRED ON
 		RUNTIME_OUTPUT_DIRECTORY_DEBUG ${dir}/bin
 		RUNTIME_OUTPUT_DIRECTORY_RELEASE ${dir}/bin
diff --git a/include/qboot/algebra/complex_function.hpp b/include/qboot/algebra/complex_function.hpp
index 973cd77..0f1b97f 100644
--- a/include/qboot/algebra/complex_function.hpp
+++ b/include/qboot/algebra/complex_function.hpp
@@ -57,15 +57,15 @@ namespace qboot::algebra
 	//   \sum_{n = 0}^{lambda / 2} \sum_{m = 0}^{lambda - 2 n} this->at(n, m) (x - x0) ^ m (y - y0) ^ n + ...
 	// we take (x0, y0) = (1 / 2, 0)
 	// if a nontrivial symmetry even (resp. odd) is given, m runs over even (resp. odd) number only.
-	// Ring must be mp::real or Polynomial
-	template <class Ring>
+	// R must be mp::real or Polynomial
+	template <Ring R>
 	class ComplexFunction
 	{
-		template <class Ring2>
+		template <Ring R2>
 		friend class ComplexFunction;
 		FunctionSymmetry sym_ = FunctionSymmetry::Mixed;
 		uint32_t lambda_;
-		Vector<Ring> coeffs_;
+		Vector<R> coeffs_;
 		// indices are aligned by lexicographical order of (dy, dx)
 		// i.e., (dx, dy) is the index(dx, dy)-th element in
 		//   [(dx, dy) for dy in range(lambda // 2) for dx in range(lambda - 2 * dy) if sym is Mixed or dx % 2 == sym]
@@ -79,7 +79,7 @@ namespace qboot::algebra
 			default: return (lambda_ + 2 - dy) * dy + dx;
 			}
 		}
-		ComplexFunction(Vector<Ring>&& v, uint32_t l, FunctionSymmetry sym)
+		ComplexFunction(Vector<R>&& v, uint32_t l, FunctionSymmetry sym)
 		    : sym_(sym), lambda_(l), coeffs_(std::move(v))
 		{
 		}
@@ -112,9 +112,9 @@ namespace qboot::algebra
 		// get coefficient of the term (x - x0) ^ {dx} (y - y0) ^ {dy}
 		// 0 <= dx, dy and dx + 2 dy <= lambda
 		// if symmetry is even or odd, the parity of dx must equals symmetry
-		[[nodiscard]] Ring& at(uint32_t dx, uint32_t dy) & { return coeffs_[index(dx, dy)]; }
-		[[nodiscard]] const Ring& at(uint32_t dx, uint32_t dy) const& { return coeffs_[index(dx, dy)]; }
-		[[nodiscard]] Vector<Ring> flatten() && { return std::move(coeffs_); }
+		[[nodiscard]] R& at(uint32_t dx, uint32_t dy) & { return coeffs_[index(dx, dy)]; }
+		[[nodiscard]] const R& at(uint32_t dx, uint32_t dy) const& { return coeffs_[index(dx, dy)]; }
+		[[nodiscard]] Vector<R> flatten() && { return std::move(coeffs_); }
 
 		[[nodiscard]] auto abs() const { return coeffs_.abs(); }
 		[[nodiscard]] auto norm() const { return coeffs_.norm(); }
@@ -144,14 +144,14 @@ namespace qboot::algebra
 			coeffs_ -= f.coeffs_;
 			return *this;
 		}
-		template <class T>
-		ComplexFunction& operator*=(const T& r) &
+		template <class S>
+		ComplexFunction& operator*=(const S& r) &
 		{
 			coeffs_ *= r;
 			return *this;
 		}
-		template <class T>
-		ComplexFunction& operator/=(const T& r) &
+		template <class S>
+		ComplexFunction& operator/=(const S& r) &
 		{
 			coeffs_ /= r;
 			return *this;
@@ -204,23 +204,23 @@ namespace qboot::algebra
 				}
 			return z;
 		}
-		template <class R>
-		friend ComplexFunction mul_scalar(const R& r, const ComplexFunction& x)
+		template <class S>
+		friend ComplexFunction mul_scalar(const S& r, const ComplexFunction& x)
 		{
 			return ComplexFunction(mul_scalar(r, x.coeffs_), x.lambda_, x.sym_);
 		}
-		template <class R>
-		friend ComplexFunction mul_scalar(const R& r, ComplexFunction&& x)
+		template <class S>
+		friend ComplexFunction mul_scalar(const S& r, ComplexFunction&& x)
 		{
 			return std::move(x *= r);
 		}
-		template <class R>
-		friend ComplexFunction operator/(const ComplexFunction& x, const R& r)
+		template <class S>
+		friend ComplexFunction operator/(const ComplexFunction& x, const S& r)
 		{
 			return ComplexFunction(x.coeffs_ / r, x.lambda_, x.sym_);
 		}
-		template <class R>
-		friend ComplexFunction operator/(ComplexFunction&& x, const R& r)
+		template <class S>
+		friend ComplexFunction operator/(ComplexFunction&& x, const S& r)
 		{
 			return std::move(x /= r);
 		}
@@ -229,22 +229,22 @@ namespace qboot::algebra
 			return x.lambda_ == y.lambda_ && x.sym_ == y.sym_ && x.coeffs_ == y.coeffs_;
 		}
 		friend bool operator!=(const ComplexFunction& x, const ComplexFunction& y) { return !(x == y); }
-		[[nodiscard]] ComplexFunction<_evaluated_t<Ring>> eval(const mp::real& x) const
+		[[nodiscard]] ComplexFunction<_evaluated_t<R>> eval(const mp::real& x) const
 		{
-			return ComplexFunction<_evaluated_t<Ring>>(coeffs_.eval(x), lambda_, sym_);
+			return ComplexFunction<_evaluated_t<R>>(coeffs_.eval(x), lambda_, sym_);
 		}
 	};
 	template <>
 	ComplexFunction<mp::real> ComplexFunction<mp::real>::proj(FunctionSymmetry sym) &&;
 	template <>
 	ComplexFunction<Polynomial> ComplexFunction<Polynomial>::proj(FunctionSymmetry sym) &&;
-	template <class Ring>
-	ComplexFunction<_polynomialize_t<Ring>> to_pol(Vector<ComplexFunction<Ring>>* coeffs)
+	template <Ring R>
+	ComplexFunction<_polynomialize_t<R>> to_pol(Vector<ComplexFunction<R>>* coeffs)
 	{
 		uint32_t lambda = coeffs->at(0).lambda(), len = coeffs->size();
 		auto sym = coeffs->at(0).symmetry();
-		ComplexFunction<_polynomialize_t<Ring>> ans(lambda, sym);
-		Vector<Ring> v(len);
+		ComplexFunction<_polynomialize_t<R>> ans(lambda, sym);
+		Vector<R> v(len);
 		for (uint32_t dy = 0; dy <= lambda / 2; ++dy)
 			for (uint32_t dx = 0; dx + 2 * dy <= lambda; ++dx)
 				if (_matches(sym, dx))
@@ -254,7 +254,7 @@ namespace qboot::algebra
 				}
 		return ans;
 	}
-	template <class R>
+	template <Ring R>
 	std::ostream& operator<<(std::ostream& out, const ComplexFunction<R>& v)
 	{
 		out << "[";
diff --git a/include/qboot/algebra/matrix.hpp b/include/qboot/algebra/matrix.hpp
index 8a4e8e2..115fcaa 100644
--- a/include/qboot/algebra/matrix.hpp
+++ b/include/qboot/algebra/matrix.hpp
@@ -2,36 +2,36 @@
 #define QBOOT_ALGEBRA_MATRIX_HPP_
 
 #include <cassert>           // for assert
+#include <concepts>          // for same_as
 #include <cstdint>           // for uint32_t
 #include <initializer_list>  // for initializer_list
 #include <memory>            // for unique_ptr, make_unique
 #include <ostream>           // for ostream
-#include <type_traits>       // for true_type, false_type, is_same_v, enable_if, void_t
 #include <utility>           // for move, swap
 
-#include "qboot/mp/real.hpp"  // for real
+#include "qboot/mp/real.hpp"      // for real
+#include "qboot/my_concepts.hpp"  // for default_initializable
 
 namespace qboot::algebra
 {
-	template <class, template <class> class, class = std::void_t<>>
-	struct _detect : std::false_type
-	{
-	};
-	template <class T, template <class> class Check>
-	struct _detect<T, Check, std::void_t<Check<T>>> : std::true_type
+	template <class T>
+	concept _zero_checkable = requires(const T c)
 	{
+		// assert(T{}.iszero());
+		{
+			c.iszero()
+		}
+		->std::same_as<bool>;
 	};
-	template <class T>
-	using _has_iszero_checker = decltype(std::declval<const T&>().iszero());
-	template <class T>
-	inline constexpr bool _has_iszero = _detect<T, _has_iszero_checker>::value;
 	template <class R>
 	bool iszero(const R& v)
 	{
-		if constexpr (_has_iszero<R>)
-			return v.iszero();
-		else
-			return v == 0;
+		return v == 0;
+	}
+	template <_zero_checkable R>
+	bool iszero(const R& v)
+	{
+		return v.iszero();
 	}
 	template <class T>
 	struct _evaluated;
@@ -47,7 +47,7 @@ namespace qboot::algebra
 	{
 		using type = mp::real;
 	};
-	template <class Ring, class S>
+	template <class R, class S>
 	struct _substitute;
 	// substitute the most inner template argument R by S,
 	// i.e. _substitute_t<A<B<...<C<R>>...>>> = A<B<...<C<S>>...>>
@@ -64,25 +64,58 @@ namespace qboot::algebra
 		using type = Vec<_substitute_t<R, S>>;
 	};
 	inline mp::real eval(const mp::real& v, [[maybe_unused]] const mp::real& x) { return v; }
-	// interface Swappable<T> {
-	//	void swap(T&);
-	// };
-	// interface Cloneable<T> {
-	//	T clone() const;
-	// };
-	// interface ZeroCheckable<T> {
-	//	bool iszero() const;
-	// }
-	// assert(T{}.iszero());  // where T: DefaultConstructible, ZeroCheckable
-	// Ring: Swappable, Clonable, ZeroCheckable, DefaultConstructible
-	// Polynomial: Swappable, Clonable, ZeroCheckable, DefaultConstructible
-	// Vector<Ring>: Swappable, Clonable
-	// Matrix<Ring>: Swappable, Clonable
-	// definitions
-	template <class Ring>
+
+	template <class T>
+	concept _ring_base = requires(T x, T y, const T c)
+	{
+		{
+			x.swap(y)
+		}
+		->std::same_as<void>;
+		{
+			c.clone()
+		}
+		->std::same_as<T>;
+	}
+	&&_zero_checkable<T>&& default_initializable<T>;
+	template <class T>
+	concept _ring_ops = requires(const T& x, const T& y, T& r)
+	{
+		r += y;
+		r -= y;
+		r.negate();
+		{
+			x + y
+		}
+		->std::same_as<T>;
+		{
+			+x
+		}
+		->std::same_as<T>;
+		{
+			-x
+		}
+		->std::same_as<T>;
+		{
+			x == y
+		}
+		->std::same_as<bool>;
+		x.norm();
+	};
+	template <class T>
+	concept Ring = _ring_base<T>&& _ring_ops<T>;
+	template <class T>
+	concept Algebra = requires(const T x, const T y)
+	{
+		{
+			mul(x, y)
+		}
+		->std::same_as<T>;
+	};
+	template <Ring R>
 	class Vector
 	{
-		std::unique_ptr<Ring[]> arr_;
+		std::unique_ptr<R[]> arr_;
 		uint32_t sz_;
 
 	public:
@@ -92,8 +125,8 @@ namespace qboot::algebra
 			sz_ = 0;
 		}
 		Vector() : arr_{}, sz_(0) {}
-		explicit Vector(uint32_t len) : arr_(std::make_unique<Ring[]>(len)), sz_(len) {}
-		Vector(uint32_t len, const Ring& val) : Vector(len)
+		explicit Vector(uint32_t len) : arr_(std::make_unique<R[]>(len)), sz_(len) {}
+		Vector(uint32_t len, const R& val) : Vector(len)
 		{
 			for (uint32_t i = 0; i < len; ++i) arr_[i] = val.clone();
 		}
@@ -102,20 +135,20 @@ namespace qboot::algebra
 		~Vector() = default;
 		Vector(const Vector& v) = delete;
 		Vector& operator=(const Vector& v) = delete;
-		Vector(std::initializer_list<Ring> v) : Vector(uint32_t(v.size()))
+		Vector(std::initializer_list<R> v) : Vector(uint32_t(v.size()))
 		{
 			uint32_t i = 0;
 			for (auto& t : v) arr_[i++] = t.clone();
 		}
-		[[nodiscard]] Ring& at(uint32_t i) & { return arr_[i]; }
-		[[nodiscard]] const Ring& at(uint32_t i) const& { return arr_[i]; }
-		[[nodiscard]] Ring& operator[](uint32_t i) & { return at(i); }
-		[[nodiscard]] const Ring& operator[](uint32_t i) const& { return at(i); }
+		[[nodiscard]] R& at(uint32_t i) & { return arr_[i]; }
+		[[nodiscard]] const R& at(uint32_t i) const& { return arr_[i]; }
+		[[nodiscard]] R& operator[](uint32_t i) & { return at(i); }
+		[[nodiscard]] const R& operator[](uint32_t i) const& { return at(i); }
 		[[nodiscard]] const uint32_t& size() const noexcept { return sz_; }
-		[[nodiscard]] const Ring* begin() const& noexcept { return arr_.get(); }
-		[[nodiscard]] const Ring* end() const& noexcept { return arr_.get() + sz_; }
-		[[nodiscard]] Ring* begin() & noexcept { return arr_.get(); }
-		[[nodiscard]] Ring* end() & noexcept { return arr_.get() + sz_; }
+		[[nodiscard]] const R* begin() const& noexcept { return arr_.get(); }
+		[[nodiscard]] const R* end() const& noexcept { return arr_.get() + sz_; }
+		[[nodiscard]] R* begin() & noexcept { return arr_.get(); }
+		[[nodiscard]] R* end() & noexcept { return arr_.get() + sz_; }
 		[[nodiscard]] Vector clone() const
 		{
 			Vector v(sz_);
@@ -156,14 +189,14 @@ namespace qboot::algebra
 			for (uint32_t i = 0; i < sz_; ++i) arr_[i] -= v.arr_[i];
 			return *this;
 		}
-		template <class T>
-		Vector& operator*=(const T& v) &
+		template <class S>
+		Vector& operator*=(const S& v) &
 		{
 			for (uint32_t i = 0; i < sz_; ++i) arr_[i] *= v;
 			return *this;
 		}
-		template <class T>
-		Vector& operator/=(const T& v) &
+		template <class S>
+		Vector& operator/=(const S& v) &
 		{
 			for (uint32_t i = 0; i < sz_; ++i) arr_[i] /= v;
 			return *this;
@@ -203,27 +236,27 @@ namespace qboot::algebra
 			std::move(y)._reset();
 			return std::move(x);
 		}
-		template <class R>
-		friend Vector mul_scalar(const R& r, const Vector& x)
+		template <class S>
+		friend Vector mul_scalar(const S& r, const Vector& x)
 		{
 			Vector z(x.sz_);
 			for (uint32_t i = 0; i < x.sz_; ++i) z[i] = mul_scalar(r, x.arr_[i]);
 			return z;
 		}
-		template <class R>
-		friend Vector mul_scalar(const R& r, Vector&& x)
+		template <class S>
+		friend Vector mul_scalar(const S& r, Vector&& x)
 		{
 			return std::move(x *= r);
 		}
-		template <class R>
-		friend Vector operator/(const Vector& x, const R& r)
+		template <class S>
+		friend Vector operator/(const Vector& x, const S& r)
 		{
 			Vector z(x.sz_);
 			for (uint32_t i = 0; i < x.sz_; ++i) z.arr_[i] = x.arr_[i] / r;
 			return z;
 		}
-		template <class R>
-		friend Vector operator/(Vector&& x, const R& r)
+		template <class S>
+		friend Vector operator/(Vector&& x, const S& r)
 		{
 			return std::move(x /= r);
 		}
@@ -248,31 +281,30 @@ namespace qboot::algebra
 			negate();
 			return std::move(*this);
 		}
-		template <class Ring2 = Ring, class = std::enable_if<std::is_same_v<Ring, Ring2>>>
-		friend Ring dot(const Vector& x, const Vector<Ring2>& y)
+		friend R dot(const Vector& x, const Vector& y) requires Algebra<R>
 		{
 			assert(x.sz_ == y.sz_);
 			if (x.sz_ == 0) return {};
-			auto s = mul(x.arr_[0], y.arr_[0]);
+			R s = mul(x.arr_[0], y.arr_[0]);
 			for (uint32_t i = 1; i < x.sz_; ++i) s += mul(x.arr_[i], y.arr_[i]);
 			return s;
 		}
-		[[nodiscard]] Vector<_evaluated_t<Ring>> eval(const mp::real& x) const
+		[[nodiscard]] Vector<_evaluated_t<R>> eval(const mp::real& x) const
 		{
-			Vector<_evaluated_t<Ring>> ans(sz_);
+			Vector<_evaluated_t<R>> ans(sz_);
 			for (uint32_t i = 0; i < sz_; ++i) ans[i] = at(i).eval(x);
 			return ans;
 		}
 	};
 
-	template <class Ring>
+	template <Ring R>
 	class Matrix
 	{
-		template <class Ring2>
+		template <Ring R2>
 		friend class Matrix;
-		Vector<Ring> arr_;
+		Vector<R> arr_;
 		uint32_t row_, col_;
-		Matrix(Vector<Ring>&& v, uint32_t r, uint32_t c) : arr_(std::move(v)), row_(r), col_(c) {}
+		Matrix(Vector<R>&& v, uint32_t r, uint32_t c) : arr_(std::move(v)), row_(r), col_(c) {}
 
 	public:
 		void _reset() &&
@@ -287,8 +319,8 @@ namespace qboot::algebra
 		~Matrix() = default;
 		Matrix(const Matrix& v) = delete;
 		Matrix& operator=(const Matrix& v) = delete;
-		[[nodiscard]] Ring& at(uint32_t r, uint32_t c) & { return arr_[r * col_ + c]; }
-		[[nodiscard]] const Ring& at(uint32_t r, uint32_t c) const& { return arr_[r * col_ + c]; }
+		[[nodiscard]] R& at(uint32_t r, uint32_t c) & { return arr_[r * col_ + c]; }
+		[[nodiscard]] const R& at(uint32_t r, uint32_t c) const& { return arr_[r * col_ + c]; }
 		[[nodiscard]] const uint32_t& row() const noexcept { return row_; }
 		[[nodiscard]] const uint32_t& column() const noexcept { return col_; }
 		[[nodiscard]] bool is_square() const noexcept { return row_ == col_; }
@@ -315,26 +347,26 @@ namespace qboot::algebra
 			arr_ -= v.arr_;
 			return *this;
 		}
-		template <class T>
-		Matrix& operator*=(const T& v) &
+		template <class S>
+		Matrix& operator*=(const S& v) &
 		{
 			arr_ *= v;
 			return *this;
 		}
-		template <class T>
-		Matrix& operator/=(const T& v) &
+		template <class S>
+		Matrix& operator/=(const S& v) &
 		{
 			arr_ /= v;
 			return *this;
 		}
-		Vector<Ring> flatten() && { return std::move(arr_); }
+		Vector<R> flatten() && { return std::move(arr_); }
 		void transpose() &
 		{
 			assert(is_square());
 			for (uint32_t i = 0; i < row_; ++i)
 				for (uint32_t j = 0; j < i; ++j) at(i, j).swap(at(j, i));
 		}
-		static Matrix constant(const Ring& c, uint32_t n)
+		static Matrix constant(const R& c, uint32_t n)
 		{
 			Matrix m(n, n);
 			for (uint32_t i = 0; i < n; ++i) m.at(i, i) = c.clone();
@@ -349,8 +381,7 @@ namespace qboot::algebra
 			return std::move(*this);
 		}
 		// v^t M v
-		template <class = std::enable_if<std::is_same_v<Ring, mp::real>>>
-		[[nodiscard]] Ring inner_product(const Vector<Ring>& v) const
+		[[nodiscard]] R inner_product(const Vector<R>& v) const requires std::same_as<R, mp::real>
 		{
 			assert(is_square() && row_ == v.size());
 			mp::real s{};
@@ -389,23 +420,23 @@ namespace qboot::algebra
 			std::move(y)._reset();
 			return std::move(x);
 		}
-		template <class R>
-		friend Matrix mul_scalar(const R& r, const Matrix& x)
+		template <class S>
+		friend Matrix mul_scalar(const S& r, const Matrix& x)
 		{
 			return Matrix(mul_scalar(r, x.arr_), x.row_, x.col_);
 		}
-		template <class R>
-		friend Matrix mul_scalar(const R& r, Matrix&& x)
+		template <class S>
+		friend Matrix mul_scalar(const S& r, Matrix&& x)
 		{
 			return std::move(x *= r);
 		}
-		template <class R>
-		friend Matrix operator/(const Matrix& x, const R& r)
+		template <class S>
+		friend Matrix operator/(const Matrix& x, const S& r)
 		{
 			return Matrix(x.arr_ / r, x.row_, x.col_);
 		}
-		template <class R>
-		friend Matrix operator/(Matrix&& x, const R& r)
+		template <class S>
+		friend Matrix operator/(Matrix&& x, const S& r)
 		{
 			return std::move(x /= r);
 		}
@@ -428,11 +459,11 @@ namespace qboot::algebra
 			return z;
 		}
 		friend Matrix mul(const Matrix& x, const Matrix& y) { return dot(x, y); }
-		template <class R>
-		friend Vector<R> dot(const Matrix& x, const Vector<R>& y)
+		template <class S>
+		friend Vector<S> dot(const Matrix& x, const Vector<S>& y)
 		{
 			assert(x.col_ == y.size());
-			Vector<R> z(x.row_);
+			Vector<S> z(x.row_);
 			if (x.col_ > 0)
 				for (uint32_t i = 0; i < x.row_; ++i)
 				{
@@ -441,11 +472,11 @@ namespace qboot::algebra
 				}
 			return z;
 		}
-		template <class R>
-		friend Vector<R> dot(const Vector<R>& x, const Matrix& y)
+		template <class S>
+		friend Vector<S> dot(const Vector<S>& x, const Matrix& y)
 		{
 			assert(x.size() == y.row_);
-			Vector<R> z(y.col_);
+			Vector<S> z(y.col_);
 			if (y.row_ > 0)
 				for (uint32_t i = 0; i < y.col_; ++i)
 				{
@@ -454,9 +485,9 @@ namespace qboot::algebra
 				}
 			return z;
 		}
-		[[nodiscard]] Matrix<_evaluated_t<Ring>> eval(const mp::real& x) const
+		[[nodiscard]] Matrix<_evaluated_t<R>> eval(const mp::real& x) const
 		{
-			return Matrix<_evaluated_t<Ring>>(arr_.eval(x), row_, col_);
+			return Matrix<_evaluated_t<R>>(arr_.eval(x), row_, col_);
 		}
 	};
 
@@ -470,8 +501,8 @@ namespace qboot::algebra
 	// calculate the inverse matrix of lower triangular matrix
 	[[nodiscard]] Matrix<mp::real> lower_triangular_inverse(const Matrix<mp::real>& mat);
 
-	template <class Ring>
-	std::ostream& operator<<(std::ostream& out, const Vector<Ring>& v)
+	template <Ring R>
+	std::ostream& operator<<(std::ostream& out, const Vector<R>& v)
 	{
 		out << "[";
 		auto f = false;
@@ -484,8 +515,8 @@ namespace qboot::algebra
 		return out << "]";
 	}
 
-	template <class Ring>
-	std::ostream& operator<<(std::ostream& out, const Matrix<Ring>& v)
+	template <Ring R>
+	std::ostream& operator<<(std::ostream& out, const Matrix<R>& v)
 	{
 		out << "[";
 		auto f = false;
diff --git a/include/qboot/algebra/polynomial.hpp b/include/qboot/algebra/polynomial.hpp
index 0ac92ff..158bb3c 100644
--- a/include/qboot/algebra/polynomial.hpp
+++ b/include/qboot/algebra/polynomial.hpp
@@ -71,7 +71,7 @@ namespace qboot::algebra
 			auto s = coeff_[d];
 			for (uint32_t i = d - 1; i <= d; --i)
 			{
-				if constexpr (std::is_same_v<R, mp::real>)
+				if constexpr (std::same_as<R, mp::real>)
 					mp::fma(s, s, x, coeff_[i]);
 				else
 				{
@@ -163,14 +163,14 @@ namespace qboot::algebra
 	};
 	template <class T>
 	using _polynomialize_t = _substitute_t<T, Polynomial>;
-	// schematically, to_pol(Vector<Ring>{a, b, c, ...}) = a + b x + c x ^ 2 + ...
+	// schematically, to_pol(Vector<R>{a, b, c, ...}) = a + b x + c x ^ 2 + ...
 	inline auto to_pol(Vector<mp::real>* coeffs) { return Polynomial(coeffs->clone()); }
-	template <class Ring>
-	Matrix<_polynomialize_t<Ring>> to_pol(Vector<Matrix<Ring>>* coeffs)
+	template <Ring R>
+	Matrix<_polynomialize_t<R>> to_pol(Vector<Matrix<R>>* coeffs)
 	{
 		uint32_t row = coeffs->at(0).row(), column = coeffs->at(0).column(), len = coeffs->size();
-		Matrix<_polynomialize_t<Ring>> ans(row, column);
-		Vector<Ring> v(len);
+		Matrix<_polynomialize_t<R>> ans(row, column);
+		Vector<R> v(len);
 		for (uint32_t r = 0; r < row; ++r)
 			for (uint32_t c = 0; c < column; ++c)
 			{
@@ -179,12 +179,12 @@ namespace qboot::algebra
 			}
 		return ans;
 	}
-	template <class Ring>
-	Vector<_polynomialize_t<Ring>> to_pol(Vector<Vector<Ring>>* coeffs)
+	template <Ring R>
+	Vector<_polynomialize_t<R>> to_pol(Vector<Vector<R>>* coeffs)
 	{
 		uint32_t sz = coeffs->at(0).size(), len = coeffs->size();
-		Vector<_polynomialize_t<Ring>> ans(sz);
-		Vector<Ring> v(len);
+		Vector<_polynomialize_t<R>> ans(sz);
+		Vector<R> v(len);
 		for (uint32_t r = 0; r < sz; ++r)
 		{
 			for (uint32_t i = 0; i < len; ++i) v[i].swap(coeffs->at(i).at(r));
@@ -193,9 +193,8 @@ namespace qboot::algebra
 		return ans;
 	}
 	Matrix<mp::real> interpolation_matrix(const Vector<mp::real>& points);
-	template <class Ring>
-	_polynomialize_t<Ring> polynomial_interpolate(const Vector<Ring>& vals,
-	                                              const Matrix<mp::real>& interpolation_matrix)
+	template <Ring R>
+	_polynomialize_t<R> polynomial_interpolate(const Vector<R>& vals, const Matrix<mp::real>& interpolation_matrix)
 	{
 		assert(vals.size() == interpolation_matrix.row() && vals.size() > 0 && interpolation_matrix.is_square());
 		auto coeffs = dot(interpolation_matrix, vals);
@@ -204,16 +203,16 @@ namespace qboot::algebra
 	// calculate coefficients c of polynomial f(x) s.t. for each i, f(points[i]) = vals[i]
 	// vals[i] = c[0] + c[1] points[i] + c[2] points[i] ^ 2 + ... + c[deg] points[i] ^ {deg}
 	// evals(polynomial_interpolate(vals, points), points) == vals (up to rounding errors)
-	template <class Ring>
-	_polynomialize_t<Ring> polynomial_interpolate(const Vector<Ring>& vals, const Vector<mp::real>& points)
+	template <Ring R>
+	_polynomialize_t<R> polynomial_interpolate(const Vector<R>& vals, const Vector<mp::real>& points)
 	{
 		assert(vals.size() == points.size() && points.size() > 0);
 		return polynomial_interpolate(vals, interpolation_matrix(points));
 	}
-	template <class Ring>
-	Vector<_evaluated_t<Ring>> evals(const Ring& v, const Vector<mp::real>& xs)
+	template <Ring R>
+	Vector<_evaluated_t<R>> evals(const R& v, const Vector<mp::real>& xs)
 	{
-		Vector<_evaluated_t<Ring>> ans(xs.size());
+		Vector<_evaluated_t<R>> ans(xs.size());
 		for (uint32_t i = 0; i < xs.size(); ++i) ans[i] = v.eval(xs[i]);
 		return ans;
 	}
diff --git a/include/qboot/algebra/real_function.hpp b/include/qboot/algebra/real_function.hpp
index 2e8cd5d..665f75d 100644
--- a/include/qboot/algebra/real_function.hpp
+++ b/include/qboot/algebra/real_function.hpp
@@ -16,15 +16,15 @@ namespace qboot::algebra
 	// take derivatives (der x) ^ k upto k <= lambda
 	// namely, a function is represented as
 	//   \sum_{k = 0}^{lambda} this->at(k) x ^ k + O(x ^ {lambda + 1})
-	// Ring must be mp::real or Polynomial
-	template <class Ring>
+	// R must be mp::real or Polynomial
+	template <Ring R>
 	class RealFunction
 	{
-		template <class Ring2>
+		template <Ring R2>
 		friend class RealFunction;
 		uint32_t lambda_;
-		Vector<Ring> coeffs_;
-		explicit RealFunction(Vector<Ring>&& vec) : lambda_(vec.size() - 1), coeffs_(std::move(vec)) {}
+		Vector<R> coeffs_;
+		explicit RealFunction(Vector<R>&& vec) : lambda_(vec.size() - 1), coeffs_(std::move(vec)) {}
 
 	public:
 		friend class RealConverter;
@@ -48,8 +48,8 @@ namespace qboot::algebra
 		[[nodiscard]] RealFunction clone() const { return RealFunction(coeffs_.clone()); }
 		// get coefficient of the term x ^ k
 		// 0 <= k <= lambda
-		[[nodiscard]] Ring& at(uint32_t k) & { return coeffs_[k]; }
-		[[nodiscard]] const Ring& at(uint32_t k) const& { return coeffs_[k]; }
+		[[nodiscard]] R& at(uint32_t k) & { return coeffs_[k]; }
+		[[nodiscard]] const R& at(uint32_t k) const& { return coeffs_[k]; }
 
 		[[nodiscard]] auto abs() const { return coeffs_.abs(); }
 		[[nodiscard]] auto norm() const { return coeffs_.norm(); }
@@ -70,14 +70,14 @@ namespace qboot::algebra
 			coeffs_ -= v.coeffs_;
 			return *this;
 		}
-		template <class T>
-		RealFunction& operator*=(const T& r) &
+		template <class S>
+		RealFunction& operator*=(const S& r) &
 		{
 			coeffs_ *= r;
 			return *this;
 		}
-		template <class T>
-		RealFunction& operator/=(const T& r) &
+		template <class S>
+		RealFunction& operator/=(const S& r) &
 		{
 			coeffs_ /= r;
 			return *this;
@@ -126,23 +126,23 @@ namespace qboot::algebra
 				for (uint32_t k2 = 0; k1 + k2 <= x.lambda_; ++k2) z.at(k1 + k2) += mul(x.at(k1), y.at(k2));
 			return z;
 		}
-		template <class R>
-		friend RealFunction mul_scalar(const R& r, const RealFunction& x)
+		template <class S>
+		friend RealFunction mul_scalar(const S& r, const RealFunction& x)
 		{
 			return RealFunction(mul_scalar(r, x.coeffs_));
 		}
-		template <class R>
-		friend RealFunction mul_scalar(const R& r, RealFunction&& x)
+		template <class S>
+		friend RealFunction mul_scalar(const S& r, RealFunction&& x)
 		{
 			return std::move(x *= r);
 		}
-		template <class R>
-		friend RealFunction operator/(const RealFunction& x, const R& r)
+		template <class S>
+		friend RealFunction operator/(const RealFunction& x, const S& r)
 		{
 			return RealFunction(x.coeffs_ / r);
 		}
-		template <class R>
-		friend RealFunction operator/(RealFunction&& x, const R& r)
+		template <class S>
+		friend RealFunction operator/(RealFunction&& x, const S& r)
 		{
 			return std::move(x /= r);
 		}
@@ -151,17 +151,17 @@ namespace qboot::algebra
 			return x.lambda_ == y.lambda_ && x.coeffs_ == y.coeffs_;
 		}
 		friend bool operator!=(const RealFunction& x, const RealFunction& y) { return !(x == y); }
-		[[nodiscard]] RealFunction<_evaluated_t<Ring>> eval(const mp::real& x) const
+		[[nodiscard]] RealFunction<_evaluated_t<R>> eval(const mp::real& x) const
 		{
-			return RealFunction<_evaluated_t<Ring>>(coeffs_.eval(x));
+			return RealFunction<_evaluated_t<R>>(coeffs_.eval(x));
 		}
 	};
-	template <class Ring>
-	RealFunction<_polynomialize_t<Ring>> to_pol(Vector<RealFunction<Ring>>* coeffs)
+	template <Ring R>
+	RealFunction<_polynomialize_t<R>> to_pol(Vector<RealFunction<R>>* coeffs)
 	{
 		uint32_t lambda = coeffs->at(0).lambda(), len = coeffs->size();
-		RealFunction<_polynomialize_t<Ring>> ans(lambda);
-		Vector<Ring> v(len);
+		RealFunction<_polynomialize_t<R>> ans(lambda);
+		Vector<R> v(len);
 		for (uint32_t r = 0; r <= lambda; ++r)
 		{
 			for (uint32_t i = 0; i < len; ++i) v[i].swap(coeffs->at(i).at(r));
@@ -169,7 +169,7 @@ namespace qboot::algebra
 		}
 		return ans;
 	}
-	template <class R>
+	template <Ring R>
 	std::ostream& operator<<(std::ostream& out, const RealFunction<R>& v)
 	{
 		auto f = false;
@@ -205,7 +205,7 @@ namespace qboot::algebra
 		explicit RealConverter(const RealFunction<mp::real>& func);
 		[[nodiscard]] RealConverter inverse() const;
 		// convert a function f of x to a function of y where x = func(y)
-		template <class R>
+		template <Ring R>
 		[[nodiscard]] RealFunction<R> convert(const RealFunction<R>& f) const
 		{
 			assert(lambda_ == f.lambda());
@@ -251,7 +251,7 @@ namespace qboot::algebra
 			mp::real s{};
 			for (uint32_t i = f_.lambda(); i <= f_.lambda(); --i)
 			{
-				if constexpr (std::is_same_v<R, mp::real>)
+				if constexpr (std::same_as<R, mp::real>)
 					mp::fma(s, s, x, f_.at(i));
 				else
 				{
diff --git a/include/qboot/block.hpp b/include/qboot/block.hpp
index 92d46fd..caabc39 100644
--- a/include/qboot/block.hpp
+++ b/include/qboot/block.hpp
@@ -1,10 +1,10 @@
 #ifndef QBOOT_BLOCK_HPP_
 #define QBOOT_BLOCK_HPP_
 
-#include <cstdint>      // for uint32_t
-#include <string>       // for string
-#include <type_traits>  // enable_if, is_same_v
-#include <variant>      // for varinat
+#include <concepts>  // for same_as
+#include <cstdint>   // for uint32_t
+#include <string>    // for string
+#include <variant>   // for varinat
 
 #include "qboot/algebra/complex_function.hpp"  // for FunctionSymmetry
 #include "qboot/mp/real.hpp"                   // for real
@@ -22,8 +22,12 @@ namespace qboot
 	{
 		return _delta_P(d1 - d2, d3 - d4);
 	}
+	class PrimaryOperator;
+	class GeneralPrimaryOperator;
+	template <class T>
+	concept _operator = std::same_as<T, PrimaryOperator> || std::same_as<T, GeneralPrimaryOperator>;
 	// F_{\mp, op}^{d1 d2, d3 d4}
-	template <class Operator>
+	template <_operator Operator>
 	class ConformalBlock
 	{
 		Operator op_;
@@ -63,8 +67,7 @@ namespace qboot
 		    : ConformalBlock(op, op1.delta(), op2.delta(), op3.delta(), op4.delta(), sym)
 		{
 		}
-		template <class = std::enable_if<std::is_same_v<Operator, PrimaryOperator>>>
-		[[nodiscard]] const mp::real& delta() const
+		[[nodiscard]] const mp::real& delta() const requires std::same_as<Operator, PrimaryOperator>
 		{
 			return op_.delta();
 		}
@@ -76,13 +79,13 @@ namespace qboot
 		[[nodiscard]] bool include_odd() const { return d12_ != 0 && d34_ != 0; }
 		[[nodiscard]] algebra::FunctionSymmetry symmetry() const { return sym_; }
 		[[nodiscard]] const Operator& get_op() const { return op_; }
-		template <class = std::enable_if<!std::is_same_v<Operator, PrimaryOperator>>>
 		[[nodiscard]] PrimaryOperator get_op(const mp::real& delta) const
+		    requires std::same_as<Operator, GeneralPrimaryOperator>
 		{
 			return op_.fix_delta(delta);
 		}
-		template <class = std::enable_if<!std::is_same_v<Operator, PrimaryOperator>>>
 		[[nodiscard]] ConformalBlock<PrimaryOperator> fix_delta(const mp::real& delta) const
+		    requires std::same_as<Operator, GeneralPrimaryOperator>
 		{
 			return ConformalBlock<PrimaryOperator>(op_.fix_delta(delta), d12_, d34_, d23h_, sym_);
 		}
@@ -130,7 +133,7 @@ namespace qboot
 		[[nodiscard]] const mp::real& P() const { return P_; }
 		[[nodiscard]] bool include_odd() const { return d12_ != 0 && d34_ != 0; }
 		[[nodiscard]] algebra::FunctionSymmetry symmetry() const { return sym_; }
-		template <class Operator>
+		template <_operator Operator>
 		[[nodiscard]] ConformalBlock<Operator> fix_op(const Operator& op) const
 		{
 			return ConformalBlock<Operator>(op, d12_, d34_, d23h_, sym_);
diff --git a/include/qboot/bootstrap_equation.hpp b/include/qboot/bootstrap_equation.hpp
index f80a6d4..33fa2a5 100644
--- a/include/qboot/bootstrap_equation.hpp
+++ b/include/qboot/bootstrap_equation.hpp
@@ -4,6 +4,7 @@
 #include <algorithm>    // for min
 #include <array>        // for array
 #include <cstdint>      // for uint32_t
+#include <filesystem>   // for path
 #include <functional>   // for function
 #include <map>          // for map
 #include <memory>       // for unique_ptr
@@ -21,7 +22,6 @@
 #include "qboot/context.hpp"                   // for Context
 #include "qboot/mp/rational.hpp"               // for rational
 #include "qboot/mp/real.hpp"                   // for real
-#include "qboot/my_filesystem.hpp"             // for path
 #include "qboot/polynomial_program.hpp"        // for PolynomialProgram
 #include "qboot/primary_op.hpp"                // for GeneralPrimaryOperator, PrimaryOperator
 #include "qboot/task_queue.hpp"                // for _event_base
@@ -234,7 +234,7 @@ namespace qboot
 		return std::make_unique<_extremal_OPE>(maximize, target, norm);
 	}
 	// read raw functional from a result of sdpb
-	algebra::Vector<mp::real> read_raw_functional(const fs::path& y_txt);
+	algebra::Vector<mp::real> read_raw_functional(const std::filesystem::path& y_txt);
 	class Equation;
 	// one bootrap equation contains
 	// - constant terms (from the unit operator)
diff --git a/include/qboot/context.hpp b/include/qboot/context.hpp
index fd4d2fb..368c6a6 100644
--- a/include/qboot/context.hpp
+++ b/include/qboot/context.hpp
@@ -86,7 +86,7 @@ namespace qboot
 		}
 		// cutoff of the power series expansion of conformal blocks at rho = 0
 		[[nodiscard]] uint32_t n_Max() const { return n_Max_; }
-		// passed to the constructor of ComplexFunction<Ring>
+		// passed to the constructor of ComplexFunction<R>
 		[[nodiscard]] uint32_t lambda() const { return lambda_; }
 		// the dimension of the spacetime
 		[[nodiscard]] const mp::rational& dimension() const { return dim_; }
diff --git a/include/qboot/mp/integer.hpp b/include/qboot/mp/integer.hpp
index 2add642..9f1eda0 100644
--- a/include/qboot/mp/integer.hpp
+++ b/include/qboot/mp/integer.hpp
@@ -1,6 +1,8 @@
 #ifndef QBOOT_MP_INTEGER_HPP_
 #define QBOOT_MP_INTEGER_HPP_
 
+#include <compare>      // for strong_ordering
+#include <concepts>     // for same_as
 #include <istream>      // for basic_istream
 #include <limits>       // for numeric_limits
 #include <optional>     // for optional
@@ -8,11 +10,12 @@
 #include <stdexcept>    // for runtime_error
 #include <string>       // for to_string, string_literals
 #include <string_view>  // for string_view
-#include <type_traits>  // for enable_if_t, is_same_v, is_integral_v, is_signed_v, conjunction_v, integral_constant
 #include <utility>      // for move
 
-#include "gmpxx.h"
-#include "mpfr.h"
+#include <gmpxx.h>
+#include <mpfr.h>
+
+#include "qboot/my_concepts.hpp"  // for integral, signed_integral, unsigned_integral
 
 namespace qboot::mp
 {
@@ -25,26 +28,18 @@ namespace qboot::mp
 	class real;
 
 	template <class Tp>
-	inline constexpr bool _is_mp =
-	    std::is_same_v<Tp, integer> || std::is_same_v<Tp, rational> || std::is_same_v<Tp, real>;
+	concept _is_mp = std::same_as<Tp, integer> || std::same_as<Tp, rational> || std::same_as<Tp, real>;
 
 	// check all values in integral class I1 are included in integral class I2
 	// and I1, I2 has the same signed property
 	template <class I1, class I2>
-	inline constexpr bool _is_included_v = std::numeric_limits<I2>::min() <= std::numeric_limits<I1>::min() &&
-	                                       std::numeric_limits<I1>::max() <= std::numeric_limits<I2>::max();
-
-	template <class I1, class I2>
-	struct _is_included : std::integral_constant<bool, _is_included_v<I1, I2>>
-	{
-	};
+	concept _is_included = std::numeric_limits<I2>::min() <= std::numeric_limits<I1>::min() &&
+	                       std::numeric_limits<I1>::max() <= std::numeric_limits<I2>::max();
 
 	template <class Tp>
-	inline constexpr bool _ulong_convertible_v =
-	    std::conjunction_v<std::is_integral<Tp>, std::is_unsigned<Tp>, _is_included<Tp, _ulong>>;
+	concept _ulong_convertible = unsigned_integral<Tp>&& _is_included<Tp, _ulong>;
 	template <class Tp>
-	inline constexpr bool _long_convertible_v =
-	    std::conjunction_v<std::is_integral<Tp>, std::is_signed<Tp>, _is_included<Tp, _long>>;
+	concept _long_convertible = signed_integral<Tp>&& _is_included<Tp, _long>;
 
 	inline bool _is_even(mpz_srcptr p) { return mpz_even_p(p) != 0; }           // NOLINT
 	inline bool _is_odd(mpz_srcptr p) { return mpz_odd_p(p) != 0; }             // NOLINT
@@ -54,7 +49,10 @@ namespace qboot::mp
 	inline int _cmp_si(mpq_srcptr p, _long o) { return mpq_cmp_si(p, o, 1); }   // NOLINT
 
 	template <class Tp>
-	inline constexpr bool _mpz_is_other_operands = _long_convertible_v<Tp> || _ulong_convertible_v<Tp>;
+	concept _mpz_is_other_operands = _long_convertible<Tp> || _ulong_convertible<Tp>;
+
+	template <class Tp>
+	concept _mpz_is_other_operands_d = _mpz_is_other_operands<Tp> || std::same_as<Tp, double>;
 
 	template <class Tp>
 	struct _mp_ops;
@@ -401,37 +399,10 @@ namespace qboot::mp
 		inline static int cmp(mpfr_srcptr op1, double op2) { return mpfr_cmp_d(op1, op2); }
 	};
 
-	// generate relational operators from _cmp (using ADL)
-
-	template <class Tp1, class Tp2, class = std::enable_if_t<_is_mp<Tp1> || _is_mp<Tp2>>>
+	template <_is_mp Tp1, class Tp2>
 	inline bool operator==(const Tp1& r1, const Tp2& r2)
 	{
-		return _cmp(r1, r2) == 0;
-	}
-	template <class Tp1, class Tp2, class = std::enable_if_t<_is_mp<Tp1> || _is_mp<Tp2>>>
-	inline bool operator!=(const Tp1& r1, const Tp2& r2)
-	{
-		return _cmp(r1, r2) != 0;
-	}
-	template <class Tp1, class Tp2, class = std::enable_if_t<_is_mp<Tp1> || _is_mp<Tp2>>>
-	inline bool operator<(const Tp1& r1, const Tp2& r2)
-	{
-		return _cmp(r1, r2) < 0;
-	}
-	template <class Tp1, class Tp2, class = std::enable_if_t<_is_mp<Tp1> || _is_mp<Tp2>>>
-	inline bool operator>(const Tp1& r1, const Tp2& r2)
-	{
-		return _cmp(r1, r2) > 0;
-	}
-	template <class Tp1, class Tp2, class = std::enable_if_t<_is_mp<Tp1> || _is_mp<Tp2>>>
-	inline bool operator<=(const Tp1& r1, const Tp2& r2)
-	{
-		return _cmp(r1, r2) <= 0;
-	}
-	template <class Tp1, class Tp2, class = std::enable_if_t<_is_mp<Tp1> || _is_mp<Tp2>>>
-	inline bool operator>=(const Tp1& r1, const Tp2& r2)
-	{
-		return _cmp(r1, r2) >= 0;
+		return (r1 <=> r2) == 0;
 	}
 
 	inline std::optional<rational> parse(std::string_view str);
@@ -487,6 +458,13 @@ namespace qboot::mp
 		[[nodiscard]] bool divisible_by(const integer& o) const { return mpz_divisible_p(_x, o._x) != 0; }
 		[[nodiscard]] bool divisible_by(_ulong o) const { return mpz_divisible_ui_p(_x, o) != 0; }
 
+		[[nodiscard]] integer norm() const& { return *this * *this; }
+		[[nodiscard]] integer norm() &&
+		{
+			*this *= *this;
+			return std::move(*this);
+		}
+
 		template <class T>
 		[[nodiscard]] integer eval([[maybe_unused]] const T& x) const
 		{
@@ -508,7 +486,7 @@ namespace qboot::mp
 			return {std::move(x)};
 		}
 
-		template <class T, class = std::enable_if_t<_mpz_is_other_operands<T> || std::is_same_v<T, double>>>
+		template <_mpz_is_other_operands_d T>
 		explicit integer(T o)
 		{
 			_mp_ops<T>::init_set(_x, o);
@@ -523,7 +501,7 @@ namespace qboot::mp
 				throw std::runtime_error("in qboot::mp::integer(string_view):\n  invalid input format "s += o);
 			}
 		}
-		template <class T, class = std::enable_if_t<_mpz_is_other_operands<T> || std::is_same_v<T, double>>>
+		template <_mpz_is_other_operands_d T>
 		integer& operator=(T o) &
 		{
 			_mp_ops<T>::set(_x, o);
@@ -540,30 +518,26 @@ namespace qboot::mp
 				throw std::runtime_error("in qboot::mp::integer(string_view):\n  invalid input format "s += o);
 			return *this;
 		}
-		template <class T, class = std::enable_if_t<_mpz_is_other_operands<T> || std::is_same_v<T, double>>>
+		template <_mpz_is_other_operands_d T>
 		explicit operator T() const
 		{
 			return _mp_ops<T>::get(_x);
 		}
 
-		// _cmp(a, b) returns the sign of a - b
-
-		friend int _cmp(const integer& r1, const integer& r2) { return mpz_cmp(r1._x, r2._x); }
-		template <class T, class = std::enable_if_t<_mpz_is_other_operands<T> || std::is_same_v<T, double>>>
-		friend int _cmp(const integer& r1, T r2)
+		friend std::strong_ordering operator<=>(const integer& r1, const integer& r2)
 		{
-			return _mp_ops<T>::cmp(r1._x, r2);
+			return mpz_cmp(r1._x, r2._x) <=> 0;
 		}
-		template <class T, class = std::enable_if_t<_mpz_is_other_operands<T> || std::is_same_v<T, double>>>
-		friend int _cmp(T r1, const integer& r2)
+		template <_mpz_is_other_operands_d T>
+		friend std::strong_ordering operator<=>(const integer& r1, T r2)
 		{
-			return -_cmp(r2, r1);
+			return _mp_ops<T>::cmp(r1._x, r2) <=> 0;
 		}
 
 		template <class Tp>
 		integer& operator+=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, integer>)
+			if constexpr (std::same_as<Tp, integer>)
 				mpz_add(_x, _x, o._x);
 			else
 				_mp_ops<Tp>::add(_x, _x, o);
@@ -573,7 +547,7 @@ namespace qboot::mp
 		template <class Tp>
 		integer& operator-=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, integer>)
+			if constexpr (std::same_as<Tp, integer>)
 				mpz_sub(_x, _x, o._x);
 			else
 				_mp_ops<Tp>::sub_a(_x, _x, o);
@@ -583,7 +557,7 @@ namespace qboot::mp
 		template <class Tp>
 		integer& operator*=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, integer>)
+			if constexpr (std::same_as<Tp, integer>)
 				mpz_mul(_x, _x, o._x);
 			else
 				_mp_ops<Tp>::mul(_x, _x, o);
@@ -719,73 +693,73 @@ namespace qboot::mp
 		}
 		friend integer operator<<(integer&& a, mp_bitcnt_t o) { return std::move(a <<= o); }
 
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator+(const integer& r1, const Tp& r2)
 		{
 			integer temp;
 			_mp_ops<Tp>::add(temp._x, r1._x, r2);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator+(integer&& r1, const Tp& r2)
 		{
 			return std::move(r1 += r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator+(const Tp& r1, const integer& r2)
 		{
 			return r2 + r1;
 		}
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator+(const Tp& r1, integer&& r2)
 		{
 			return std::move(r2 += r1);
 		}
 
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator-(const integer& r1, const Tp& r2)
 		{
 			integer temp;
 			_mp_ops<Tp>::sub_a(temp._x, r1._x, r2);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator-(integer&& r1, const Tp& r2)
 		{
 			return std::move(r1 -= r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator-(const Tp& r1, const integer& r2)
 		{
 			integer temp;
 			_mp_ops<Tp>::sub_b(temp._x, r1, r2._x);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator-(const Tp& r1, integer&& r2)
 		{
 			_mp_ops<Tp>::sub_b(r2._x, r1, r2._x);
 			return std::move(r2);
 		}
 
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator*(const integer& r1, const Tp& r2)
 		{
 			integer temp;
 			_mp_ops<Tp>::mul(temp._x, r1._x, r2);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator*(integer&& r1, const Tp& r2)
 		{
 			return std::move(r1 *= r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator*(const Tp& r1, const integer& r2)
 		{
 			return r2 * r1;
 		}
-		template <class Tp, class = std::enable_if_t<_mpz_is_other_operands<Tp>>>
+		template <_mpz_is_other_operands Tp>
 		friend integer operator*(const Tp& r1, integer&& r2)
 		{
 			return std::move(r2 *= r1);
diff --git a/include/qboot/mp/rational.hpp b/include/qboot/mp/rational.hpp
index 079c3b0..7c9d8ff 100644
--- a/include/qboot/mp/rational.hpp
+++ b/include/qboot/mp/rational.hpp
@@ -3,25 +3,28 @@
 
 #include <array>        // for array
 #include <cassert>      // for assert
+#include <compare>      // for strong_ordering
 #include <istream>      // for basic_istream
 #include <optional>     // for optional
 #include <ostream>      // for basic_ostream
 #include <stdexcept>    // for runtime_error
 #include <string>       // for to_string, string_literals
 #include <string_view>  // for string_view
-#include <type_traits>  // for enable_if_t, is_same_v, is_integral_v, is_signed_v, is_unsigned_v
 #include <utility>      // for move
 #include <vector>       // for vector
 
-#include "gmpxx.h"
-#include "mpfr.h"
+#include <gmpxx.h>
+#include <mpfr.h>
 
 #include "qboot/mp/integer.hpp"
 
 namespace qboot::mp
 {
 	template <class Tp>
-	inline constexpr bool _mpq_is_other_operands = _mpz_is_other_operands<Tp> || std::is_same_v<Tp, integer>;
+	concept _mpq_is_other_operands = _mpz_is_other_operands<Tp> || std::same_as<Tp, integer>;
+
+	template <class Tp>
+	concept _mpq_is_other_operands_d = _mpq_is_other_operands<Tp> || std::same_as<Tp, double>;
 
 	inline integer ceil(const rational& q);
 	inline integer floor(const rational& q);
@@ -83,6 +86,13 @@ namespace qboot::mp
 		}
 		[[nodiscard]] bool isinteger() const { return _cmp_ui(mpq_denref(_x), 1) == 0; }
 
+		[[nodiscard]] rational norm() const& { return *this * *this; }
+		[[nodiscard]] rational norm() &&
+		{
+			*this *= *this;
+			return std::move(*this);
+		}
+
 		template <class T>
 		[[nodiscard]] rational eval([[maybe_unused]] const T& x) const
 		{
@@ -124,16 +134,15 @@ namespace qboot::mp
 			return {std::move(x)};
 		}
 
-		template <class T, class = std::enable_if_t<_mpq_is_other_operands<T> || std::is_same_v<T, double>>>
+		template <_mpq_is_other_operands_d T>
 		explicit rational(T o) : rational()
 		{
 			_mp_ops<T>::set(_x, o);
 		}
-		template <class T1, class T2,
-		          class = std::enable_if_t<std::is_integral_v<T1> && std::is_integral_v<T2> && std::is_unsigned_v<T2>>>
+		template <integral T1, unsigned_integral T2>
 		rational(T1 num, T2 den) : rational()
 		{
-			if constexpr (std::is_signed_v<T1>)
+			if constexpr (signed_integral<T1>)
 				mpq_set_si(_x, num, den);
 			else
 				mpq_set_ui(_x, num, den);
@@ -158,7 +167,7 @@ namespace qboot::mp
 			mpq_canonicalize(_x);
 		}
 
-		template <class T, class = std::enable_if_t<_mpq_is_other_operands<T> || std::is_same_v<T, double>>>
+		template <_mpq_is_other_operands_d T>
 		rational& operator=(T o) &
 		{
 			_mp_ops<T>::set(_x, o);
@@ -177,7 +186,7 @@ namespace qboot::mp
 			return *this;
 		}
 
-		template <class T, class = std::enable_if_t<_mpz_is_other_operands<T> || std::is_same_v<T, double>>>
+		template <_mpz_is_other_operands_d T>
 		explicit operator T() const
 		{
 			return _mp_ops<T>::get(_x);
@@ -185,24 +194,20 @@ namespace qboot::mp
 
 		explicit operator integer() const { return floor(*this); }
 
-		// _cmp(a, b) returns the sign of a - b
-
-		friend int _cmp(const rational& r1, const rational& r2) { return mpq_cmp(r1._x, r2._x); }
-		template <class T, class = std::enable_if_t<_mpq_is_other_operands<T>>>
-		friend int _cmp(const rational& r1, T r2)
+		friend std::strong_ordering operator<=>(const rational& r1, const rational& r2)
 		{
-			return _mp_ops<T>::cmp(r1._x, r2);
+			return mpq_cmp(r1._x, r2._x) <=> 0;
 		}
-		template <class T, class = std::enable_if_t<_mpq_is_other_operands<T>>>
-		friend int _cmp(T r1, const rational& r2)
+		template <_mpq_is_other_operands T>
+		friend std::strong_ordering operator<=>(const rational& r1, T r2)
 		{
-			return -_cmp(r2, r1);
+			return _mp_ops<T>::cmp(r1._x, r2) <=> 0;
 		}
 
 		template <class Tp>
 		rational& operator+=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, rational>)
+			if constexpr (std::same_as<Tp, rational>)
 				mpq_add(_x, _x, o._x);
 			else
 				_mp_ops<Tp>::add(_x, _x, o);
@@ -212,7 +217,7 @@ namespace qboot::mp
 		template <class Tp>
 		rational& operator-=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, rational>)
+			if constexpr (std::same_as<Tp, rational>)
 				mpq_sub(_x, _x, o._x);
 			else
 				_mp_ops<Tp>::sub_a(_x, _x, o);
@@ -222,7 +227,7 @@ namespace qboot::mp
 		template <class Tp>
 		rational& operator*=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, rational>)
+			if constexpr (std::same_as<Tp, rational>)
 				mpq_mul(_x, _x, o._x);
 			else
 			{
@@ -235,7 +240,7 @@ namespace qboot::mp
 		template <class Tp>
 		rational& operator/=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, rational>)
+			if constexpr (std::same_as<Tp, rational>)
 				mpq_div(_x, _x, o._x);
 			else
 			{
@@ -322,56 +327,56 @@ namespace qboot::mp
 			return std::move(a);
 		}
 
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator+(const rational& r1, const Tp& r2)
 		{
 			rational temp;
 			_mp_ops<Tp>::add(temp._x, r1._x, r2);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator+(rational&& r1, const Tp& r2)
 		{
 			return std::move(r1 += r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator+(const Tp& r1, const rational& r2)
 		{
 			return r2 + r1;
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator+(const Tp& r1, rational&& r2)
 		{
 			return std::move(r2 += r1);
 		}
 
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator-(const rational& r1, const Tp& r2)
 		{
 			rational temp;
 			_mp_ops<Tp>::sub_a(temp._x, r1._x, r2);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator-(rational&& r1, const Tp& r2)
 		{
 			return std::move(r1 -= r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator-(const Tp& r1, const rational& r2)
 		{
 			rational temp;
 			_mp_ops<Tp>::sub_b(temp._x, r1, r2._x);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator-(const Tp& r1, rational&& r2)
 		{
 			_mp_ops<Tp>::sub_b(r2._x, r1, r2._x);
 			return std::move(r2);
 		}
 
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator*(const rational& r1, const Tp& r2)
 		{
 			rational temp;
@@ -379,23 +384,23 @@ namespace qboot::mp
 			mpq_canonicalize(temp._x);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator*(rational&& r1, const Tp& r2)
 		{
 			return std::move(r1 *= r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator*(const Tp& r1, const rational& r2) noexcept
 		{
 			return std::move(r2 * r1);
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator*(const Tp& r1, rational&& r2)
 		{
 			return std::move(r2 *= r1);
 		}
 
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator/(const rational& r1, const Tp& r2)
 		{
 			rational temp;
@@ -403,12 +408,12 @@ namespace qboot::mp
 			mpq_canonicalize(temp._x);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator/(rational&& r1, const Tp& r2)
 		{
 			return std::move(r1 /= r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator/(const Tp& r1, const rational& r2)
 		{
 			rational temp;
@@ -416,7 +421,7 @@ namespace qboot::mp
 			mpq_canonicalize(temp._x);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpq_is_other_operands<Tp>>>
+		template <_mpq_is_other_operands Tp>
 		friend rational operator/(const Tp& r1, rational&& r2)
 		{
 			_mp_ops<Tp>::div_b(r2._x, r1, r2._x);
diff --git a/include/qboot/mp/real.hpp b/include/qboot/mp/real.hpp
index e9699af..0548d22 100644
--- a/include/qboot/mp/real.hpp
+++ b/include/qboot/mp/real.hpp
@@ -33,6 +33,8 @@
 #ifndef QBOOT_MP_REAL_HPP_
 #define QBOOT_MP_REAL_HPP_
 
+#include <compare>      // for strong_ordering
+#include <concepts>     // for same_as
 #include <cstddef>      // for size_t
 #include <cstdint>      // for intmax_t
 #include <cstring>      // for strlen
@@ -46,11 +48,10 @@
 #include <stdexcept>    // for runtime_error
 #include <string>       // for string, to_string, basic_string, string_literals
 #include <string_view>  // for string_view
-#include <type_traits>  // for enable_if_t, is_same_v, is_integral_v, is_signed_v
 #include <utility>      // for move
 
-#include "gmpxx.h"
-#include "mpfr.h"
+#include <gmpxx.h>
+#include <mpfr.h>
 
 #include "qboot/mp/integer.hpp"
 #include "qboot/mp/rational.hpp"
@@ -75,8 +76,11 @@ namespace qboot::mp
 	inline int _signbit(mpfr_srcptr x) { return mpfr_signbit(x); }             // NOLINT
 
 	template <class Tp>
-	inline constexpr bool _mpfr_is_other_operands = _mpz_is_other_operands<Tp> || std::is_same_v<Tp, integer> ||
-	                                                std::is_same_v<Tp, rational> || std::is_same_v<Tp, double>;
+	concept _mpfr_is_other_operands = _mpz_is_other_operands<Tp> || std::same_as<Tp, integer> ||
+	                                  std::same_as<Tp, rational> || std::same_as<Tp, double>;
+
+	template <class Tp>
+	concept _mpfr_is_other_operands_int = _mpfr_is_other_operands<Tp> || integral<Tp>;
 
 	inline void _reset(mpfr_ptr x)
 	{
@@ -253,32 +257,32 @@ namespace qboot::mp
 			}
 		}
 
-		template <class T, class = std::enable_if_t<std::is_integral_v<T> || _mpfr_is_other_operands<T>>>
+		template <_mpfr_is_other_operands_int T>
 		explicit real(const T& o)
 		{
 			mpfr_init2(_x, global_prec);
 			if constexpr (_mpfr_is_other_operands<T>)
 				_mp_ops<T>::set(_x, o, global_rnd);
-			else if constexpr (std::is_signed_v<T>)
+			else if constexpr (signed_integral<T>)
 				mpfr_set_sj(_x, o, global_rnd);
 			else
 				mpfr_set_uj(_x, o, global_rnd);
 		}
 
-		template <class T, class = std::enable_if_t<std::is_integral_v<T>>>
+		template <integral T>
 		real(T op, mpfr_exp_t e)
 		{
 			mpfr_init2(_x, global_prec);
-			if constexpr (std::is_signed_v<T>)
+			if constexpr (signed_integral<T>)
 			{
-				if constexpr (_is_included_v<T, _long>)
+				if constexpr (_is_included<T, _long>)
 					mpfr_set_si_2exp(_x, op, e, global_rnd);
 				else
 					mpfr_set_sj_2exp(_x, op, e, global_rnd);
 			}
 			else
 			{
-				if constexpr (_is_included_v<T, _ulong>)
+				if constexpr (_is_included<T, _ulong>)
 					mpfr_set_ui_2exp(_x, op, e, global_rnd);
 				else
 					mpfr_set_uj_2exp(_x, op, e, global_rnd);
@@ -287,12 +291,12 @@ namespace qboot::mp
 
 		// converting assignment operators
 
-		template <class T, class = std::enable_if_t<std::is_integral_v<T> || _mpfr_is_other_operands<T>>>
+		template <_mpfr_is_other_operands_int T>
 		real& operator=(const T& o) &
 		{
 			if constexpr (_mpfr_is_other_operands<T>)
 				_mp_ops<T>::set(_x, o, global_rnd);
-			else if constexpr (std::is_signed_v<T>)
+			else if constexpr (signed_integral<T>)
 				mpfr_set_sj(_x, o, global_rnd);
 			else
 				mpfr_set_uj(_x, o, global_rnd);
@@ -306,24 +310,17 @@ namespace qboot::mp
 		// generic operators
 		/////////////////////////////////////////////////////////////////
 
-		// _cmp(a, b) returns the sign of a - b
-
-		friend int _cmp(const real& r1, const real& r2) { return mpfr_cmp(r1._x, r2._x); }
-		template <class T, class = std::enable_if_t<_mpfr_is_other_operands<T>>>
-		friend int _cmp(const real& r1, T r2)
-		{
-			return _mp_ops<T>::cmp(r1._x, r2);
-		}
-		template <class T, class = std::enable_if_t<_mpfr_is_other_operands<T>>>
-		friend int _cmp(T r1, const real& r2)
+		friend std::strong_ordering operator<=>(const real& r1, const real& r2) { return mpfr_cmp(r1._x, r2._x) <=> 0; }
+		template <_mpfr_is_other_operands T>
+		friend std::strong_ordering operator<=>(const real& r1, T r2)
 		{
-			return -_cmp(r2, r1);
+			return _mp_ops<T>::cmp(r1._x, r2) <=> 0;
 		}
 
 		template <class Tp>
 		real& operator+=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, real>)
+			if constexpr (std::same_as<Tp, real>)
 				mpfr_add(_x, _x, o._x, global_rnd);
 			else
 				_mp_ops<Tp>::add(_x, _x, o, global_rnd);
@@ -333,7 +330,7 @@ namespace qboot::mp
 		template <class Tp>
 		real& operator-=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, real>)
+			if constexpr (std::same_as<Tp, real>)
 				mpfr_sub(_x, _x, o._x, global_rnd);
 			else
 				_mp_ops<Tp>::sub_a(_x, _x, o, global_rnd);
@@ -343,7 +340,7 @@ namespace qboot::mp
 		template <class Tp>
 		real& operator*=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, real>)
+			if constexpr (std::same_as<Tp, real>)
 				mpfr_mul(_x, _x, o._x, global_rnd);
 			else
 				_mp_ops<Tp>::mul(_x, _x, o, global_rnd);
@@ -353,7 +350,7 @@ namespace qboot::mp
 		template <class Tp>
 		real& operator/=(const Tp& o) &
 		{
-			if constexpr (std::is_same_v<Tp, real>)
+			if constexpr (std::same_as<Tp, real>)
 				mpfr_div(_x, _x, o._x, global_rnd);
 			else
 				_mp_ops<Tp>::div_a(_x, _x, o, global_rnd);
@@ -421,98 +418,98 @@ namespace qboot::mp
 		}
 		friend real operator/(real&& r1, real&& r2) { return std::move(r1 /= r2); }
 
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator+(const real& r1, const Tp& r2)
 		{
 			real temp;
 			_mp_ops<Tp>::add(temp._x, r1._x, r2, global_rnd);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator+(real&& r1, const Tp& r2)
 		{
 			return std::move(r1 += r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator+(const Tp& r1, const real& r2)
 		{
 			return r2 + r1;
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator+(const Tp& r1, real&& r2)
 		{
 			return std::move(r2 += r1);
 		}
 
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator-(const real& r1, const Tp& r2)
 		{
 			real temp;
 			_mp_ops<Tp>::sub_a(temp._x, r1._x, r2, global_rnd);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator-(real&& r1, const Tp& r2)
 		{
 			return std::move(r1 -= r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator-(const Tp& r1, const real& r2)
 		{
 			real temp;
 			_mp_ops<Tp>::sub_b(temp._x, r1, r2._x, global_rnd);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator-(const Tp& r1, real&& r2)
 		{
 			_mp_ops<Tp>::sub_b(r2._x, r1, r2._x, global_rnd);
 			return std::move(r2);
 		}
 
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator*(const real& r1, const Tp& r2)
 		{
 			real temp;
 			_mp_ops<Tp>::mul(temp._x, r1._x, r2, global_rnd);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator*(real&& r1, const Tp& r2)
 		{
 			return std::move(r1 *= r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator*(const Tp& r1, const real& r2) noexcept
 		{
 			return std::move(r2 * r1);
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator*(const Tp& r1, real&& r2)
 		{
 			return std::move(r2 *= r1);
 		}
 
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator/(const real& r1, const Tp& r2)
 		{
 			real temp;
 			_mp_ops<Tp>::div_a(temp._x, r1._x, r2, global_rnd);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator/(real&& r1, const Tp& r2)
 		{
 			return std::move(r1 /= r2);
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator/(const Tp& r1, const real& r2)
 		{
 			real temp;
 			_mp_ops<Tp>::div_b(temp._x, r1, r2._x, global_rnd);
 			return temp;
 		}
-		template <class Tp, class = std::enable_if_t<_mpfr_is_other_operands<Tp>>>
+		template <_mpfr_is_other_operands Tp>
 		friend real operator/(const Tp& r1, real&& r2)
 		{
 			_mp_ops<Tp>::div_b(r2._x, r1, r2._x, global_rnd);
@@ -523,12 +520,12 @@ namespace qboot::mp
 		// conversion operators
 		/////////////////////////////////////////////////////////////////
 
-		template <class T, class = std::enable_if_t<std::is_integral_v<T> || _mpfr_is_other_operands<T>>>
+		template <_mpfr_is_other_operands_int T>
 		explicit operator T() const
 		{
 			if constexpr (_mpfr_is_other_operands<T>)
 				return _mp_ops<T>::get(_x, global_rnd);
-			else if constexpr (std::is_signed_v<T>)
+			else if constexpr (signed_integral<T>)
 				return mpfr_get_sj(_x, global_rnd);
 			else
 				return mpfr_get_uj(_x, global_rnd);
diff --git a/include/qboot/my_concepts.hpp b/include/qboot/my_concepts.hpp
new file mode 100644
index 0000000..98eaff4
--- /dev/null
+++ b/include/qboot/my_concepts.hpp
@@ -0,0 +1,28 @@
+#ifndef QBOOT_MY_CONCEPTS_HPP_
+#define QBOOT_MY_CONCEPTS_HPP_
+
+#include <type_traits>  // for is_constructible_v, is_integral_v, is_nothrow_destructible_v, is_signed_v
+
+namespace qboot
+{
+	// define concepts not defined by clang
+	template <class T>
+	concept integral = std::is_integral_v<T>;
+	template <class T>
+	concept signed_integral = integral<T>&& std::is_signed_v<T>;
+	template <class T>
+	concept unsigned_integral = integral<T> && !signed_integral<T>;
+
+	template <class T>
+	concept _destructible = std::is_nothrow_destructible_v<T>;
+	template <class T, class... Args>
+	concept _constructible_from = _destructible<T>&& std::is_constructible_v<T, Args...>;
+	template <class T>
+	concept default_initializable = _constructible_from<T>&& requires
+	{
+		T{};
+	}
+	&&requires { ::new (static_cast<void*>(nullptr)) T; };
+}  // namespace qboot
+
+#endif  // QBOOT_MY_CONCEPTS_HPP_
diff --git a/include/qboot/my_filesystem.hpp b/include/qboot/my_filesystem.hpp
deleted file mode 100644
index 16f3289..0000000
--- a/include/qboot/my_filesystem.hpp
+++ /dev/null
@@ -1,18 +0,0 @@
-#ifndef QBOOT_MY_FILESYSTEM_HPP_
-#define QBOOT_MY_FILESYSTEM_HPP_
-
-#if __has_include(<filesystem>)
-#include <filesystem>  // for path, create_directory
-namespace qboot
-{
-	namespace fs = std::filesystem;
-}  // namespace qboot
-#else
-#include <experimental/filesystem>
-namespace qboot
-{
-	namespace fs = std::experimental::filesystem;
-}  // namespace qboot
-#endif
-
-#endif  // QBOOT_MY_FILESYSTEM_HPP_
diff --git a/include/qboot/primary_op.hpp b/include/qboot/primary_op.hpp
index 72b2db4..157a5e6 100644
--- a/include/qboot/primary_op.hpp
+++ b/include/qboot/primary_op.hpp
@@ -1,13 +1,12 @@
 #ifndef QBOOT_PRIMARY_OP_HPP_
 #define QBOOT_PRIMARY_OP_HPP_
 
-#include <cassert>      // for assert
-#include <cstdint>      // for uint32_t
-#include <optional>     // for optional
-#include <sstream>      // for ostringstream
-#include <string>       // for string
-#include <type_traits>  // for true_type
-#include <variant>      // for variant
+#include <cassert>   // for assert
+#include <cstdint>   // for uint32_t
+#include <optional>  // for optional
+#include <sstream>   // for ostringstream
+#include <string>    // for string
+#include <variant>   // for variant
 
 #include "qboot/mp/rational.hpp"  // for rational
 #include "qboot/mp/real.hpp"      // for real, isinteger
diff --git a/include/qboot/qboot.hpp b/include/qboot/qboot.hpp
index 600169c..3e2c52c 100644
--- a/include/qboot/qboot.hpp
+++ b/include/qboot/qboot.hpp
@@ -16,7 +16,7 @@
 #include "qboot/context.hpp"
 #include "qboot/hor_formula.hpp"
 #include "qboot/hor_recursion.hpp"
-#include "qboot/my_filesystem.hpp"
+#include "qboot/my_concepts.hpp"
 #include "qboot/polynomial_program.hpp"
 #include "qboot/primary_op.hpp"
 #include "qboot/scale_factor.hpp"
diff --git a/include/qboot/sdpb_input.hpp b/include/qboot/sdpb_input.hpp
index 2ec024e..cb26e51 100644
--- a/include/qboot/sdpb_input.hpp
+++ b/include/qboot/sdpb_input.hpp
@@ -1,18 +1,18 @@
 #ifndef QBOOT_SDPB_INPUT_HPP_
 #define QBOOT_SDPB_INPUT_HPP_
 
-#include <array>     // for array
-#include <cassert>   // for assert
-#include <cstdint>   // for uint32_t
-#include <fstream>   // for ofstream
-#include <iomanip>   // for setprecision
-#include <ios>       // for fixed
-#include <memory>    // for unique_ptr
-#include <optional>  // for optional
+#include <array>       // for array
+#include <cassert>     // for assert
+#include <cstdint>     // for uint32_t
+#include <filesystem>  // for path, create_directory
+#include <fstream>     // for ofstream
+#include <iomanip>     // for setprecision
+#include <ios>         // for fixed
+#include <memory>      // for unique_ptr
+#include <optional>    // for optional
 
 #include "qboot/algebra/matrix.hpp"  // for Matrix, Vector
 #include "qboot/mp/real.hpp"         // for real
-#include "qboot/my_filesystem.hpp"   // for path, create_directory
 #include "qboot/task_queue.hpp"      // for _event_base
 
 namespace qboot
@@ -63,15 +63,15 @@ namespace qboot
 		algebra::Vector<mp::real> objectives_;
 		std::unique_ptr<std::optional<DualConstraint>[]> constraints_;
 		uint32_t num_constraints_;
-		void write_objectives(const fs::path& root) const;
+		void write_objectives(const std::filesystem::path& root) const;
 		void write_objectives(std::ostream& out) const;
-		void write_bilinear_bases(const fs::path& root) const;
+		void write_bilinear_bases(const std::filesystem::path& root) const;
 		void write_bilinear_bases(std::ostream& out) const;
-		void write_free_var_matrix(const fs::path& root, uint32_t i) const;
+		void write_free_var_matrix(const std::filesystem::path& root, uint32_t i) const;
 		void write_free_var_matrix(std::ostream& out, const DualConstraint& cons) const;
-		void write_primal_objective_c(const fs::path& root, uint32_t i) const;
+		void write_primal_objective_c(const std::filesystem::path& root, uint32_t i) const;
 		void write_primal_objective_c(std::ostream& out, const DualConstraint& cons) const;
-		void write_blocks(const fs::path& root) const;
+		void write_blocks(const std::filesystem::path& root) const;
 		void write_blocks(std::ostream& out) const;
 
 	public:
@@ -93,8 +93,10 @@ namespace qboot
 		}
 		// call this for index = 0, ..., num_constraints() - 1 before call of write
 		void register_constraint(uint32_t index, DualConstraint&& c) &;
-		void write(const fs::path& root_, uint32_t parallel = 1, const std::unique_ptr<_event_base>& event = {}) const&;
-		void write(const fs::path& root_, uint32_t parallel = 1, const std::unique_ptr<_event_base>& event = {}) &&;
+		void write(const std::filesystem::path& root_, uint32_t parallel = 1,
+		           const std::unique_ptr<_event_base>& event = {}) const&;
+		void write(const std::filesystem::path& root_, uint32_t parallel = 1,
+		           const std::unique_ptr<_event_base>& event = {}) &&;
 	};
 }  // namespace qboot
 
diff --git a/include/qboot/task_queue.hpp b/include/qboot/task_queue.hpp
index 431dcbc..f22e9ca 100644
--- a/include/qboot/task_queue.hpp
+++ b/include/qboot/task_queue.hpp
@@ -2,6 +2,7 @@
 #define QBOOT_TASK_QUEUE_HPP_
 
 #include <atomic>              // for atomic
+#include <concepts>            // for same_as
 #include <condition_variable>  // for condition_variable
 #include <functional>          // for function
 #include <future>              // for future, promise
@@ -11,13 +12,14 @@
 #include <string>              // for string
 #include <string_view>         // for string_view
 #include <thread>              // for thread
-#include <type_traits>         // for is_default_constructible_v
 #include <utility>             // for declval, move
 #include <vector>              // for vector
 
+#include "qboot/my_concepts.hpp"  // for default_initializable
+
 namespace qboot
 {
-	template <class T>
+	template <default_initializable T>
 	std::vector<T> _seq_eval(const std::vector<std::function<T()>>& fs)
 	{
 		std::vector<T> ans;
@@ -30,11 +32,10 @@ namespace qboot
 	}
 	// evaluate fs in parallel and returns its value
 	// if p <= 1, evaluate sequentially
-	template <class T>
+	template <default_initializable T>
 	std::vector<T> _parallel_evaluate(const std::vector<std::function<T()>>& fs,
 	                                  uint32_t p = std::thread::hardware_concurrency())
 	{
-		static_assert(std::is_default_constructible_v<T>);
 		if (p <= 1) return _seq_eval(fs);
 		auto N = fs.size();
 		std::vector<T> ans(N);
@@ -100,7 +101,7 @@ namespace qboot
 		~_task() override = default;
 		void run() override
 		{
-			if constexpr (std::is_same_v<T, void>)
+			if constexpr (std::same_as<T, void>)
 			{
 				f_();
 				p_.set_value();
diff --git a/include/qboot/xml_input.hpp b/include/qboot/xml_input.hpp
index 978c7eb..09aba13 100644
--- a/include/qboot/xml_input.hpp
+++ b/include/qboot/xml_input.hpp
@@ -1,13 +1,13 @@
 #ifndef QBOOT_XML_INPUT_HPP_
 #define QBOOT_XML_INPUT_HPP_
 
-#include <memory>    // for unique_ptr
-#include <optional>  // for optional
+#include <filesystem>  // for path
+#include <memory>      // for unique_ptr
+#include <optional>    // for optional
 
 #include "qboot/algebra/matrix.hpp"      // for Vector, Matrix
 #include "qboot/algebra/polynomial.hpp"  // for Polynomial
 #include "qboot/mp/real.hpp"             // for real
-#include "qboot/my_filesystem.hpp"       // for path
 
 namespace qboot
 {
@@ -54,7 +54,7 @@ namespace qboot
 		XMLInput(mp::real&& constant, algebra::Vector<mp::real>&& obj, uint32_t num_constraints);
 		[[nodiscard]] uint32_t num_constraints() const { return num_constraints_; }
 		void register_constraint(uint32_t index, PVM&& c) &;
-		void write(const fs::path& path) const;
+		void write(const std::filesystem::path& path) const;
 	};
 }  // namespace qboot
 
diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt
index b13c953..e7b5237 100644
--- a/sample/CMakeLists.txt
+++ b/sample/CMakeLists.txt
@@ -19,7 +19,7 @@ set_target_properties(sample PROPERTIES
 	DEBUG_POSTFIX "-debug"
 	COMPILE_OPTIONS "${cxx_flags}"
 	CXX_EXTENSIONS OFF
-	CXX_STANDARD 17
+	CXX_STANDARD 20
 	CXX_STANDARD_REQUIRED ON
 	RUNTIME_OUTPUT_DIRECTORY bin
 	RUNTIME_OUTPUT_DIRECTORY_DEBUG bin
diff --git a/src/algebra/complex_function.cpp b/src/algebra/complex_function.cpp
index 4bcd051..e8328d9 100644
--- a/src/algebra/complex_function.cpp
+++ b/src/algebra/complex_function.cpp
@@ -2,17 +2,17 @@
 
 #include "qboot/mp/real.hpp"  // for real, pow
 
-using qboot::algebra::ComplexFunction, qboot::algebra::FunctionSymmetry;
+using qboot::algebra::ComplexFunction, qboot::algebra::FunctionSymmetry, qboot::algebra::Ring;
 using qboot::mp::real;
 using std::move;
 
 namespace
 {
-	template <class Ring>
-	ComplexFunction<Ring> _proj(ComplexFunction<Ring>&& f, FunctionSymmetry sym)
+	template <Ring R>
+	ComplexFunction<R> _proj(ComplexFunction<R>&& f, FunctionSymmetry sym)
 	{
 		if (f.symmetry() == sym) return move(f);
-		ComplexFunction<Ring> z(f.lambda(), sym);
+		ComplexFunction<R> z(f.lambda(), sym);
 		if (f.symmetry() == FunctionSymmetry::Mixed)
 		{
 			for (uint32_t dy = 0; dy <= f.lambda() / 2; ++dy)
diff --git a/src/algebra/matrix.cpp b/src/algebra/matrix.cpp
index f90d763..581c802 100644
--- a/src/algebra/matrix.cpp
+++ b/src/algebra/matrix.cpp
@@ -9,20 +9,27 @@ using std::move;
 
 namespace qboot::algebra
 {
-	template class Vector<real>;
-	template class Matrix<real>;
-	template class Vector<Vector<real>>;
-	template class Vector<Matrix<real>>;
-	template class Vector<Polynomial>;
-	template class Matrix<Polynomial>;
-	template class Vector<RealFunction<real>>;
-	template class Matrix<RealFunction<real>>;
-	template class Vector<ComplexFunction<real>>;
-	template class Matrix<ComplexFunction<real>>;
-	template class Vector<RealFunction<Polynomial>>;
-	template class Vector<ComplexFunction<Polynomial>>;
-	template class Matrix<RealFunction<Polynomial>>;
-	template class Matrix<ComplexFunction<Polynomial>>;
+	// these instantiation causes compilation error on MSVC
+	// template class Vector<Vector<real>>;
+	// template class Vector<Matrix<real>>;
+	[[maybe_unused]] static void _instantiate()
+	{
+		Vector<real> _;
+		Matrix<real> _0;
+		// but you can use them... why?
+		Vector<Vector<real>> _1;
+		Vector<Matrix<real>> _2;
+		Vector<Polynomial> _3;
+		Matrix<Polynomial> _4;
+		Vector<RealFunction<real>> _5;
+		Matrix<RealFunction<real>> _6;
+		Vector<ComplexFunction<real>> _7;
+		Matrix<ComplexFunction<real>> _8;
+		Vector<RealFunction<Polynomial>> _9;
+		Vector<ComplexFunction<Polynomial>> _10;
+		Matrix<RealFunction<Polynomial>> _11;
+		Matrix<ComplexFunction<Polynomial>> _12;
+	}
 
 	static void _add_row(Matrix<real>* mat, uint32_t f, uint32_t t, const real& x, uint32_t c0 = 0)
 	{
diff --git a/src/block.cpp b/src/block.cpp
index 2ac3f62..9c2dde8 100644
--- a/src/block.cpp
+++ b/src/block.cpp
@@ -7,10 +7,11 @@ using std::string, std::ostringstream;
 
 namespace qboot
 {
-	template class ConformalBlock<PrimaryOperator>;
-	template class ConformalBlock<GeneralPrimaryOperator>;
+	[[maybe_unused]] static void _instantiate(ConformalBlock<PrimaryOperator>, ConformalBlock<GeneralPrimaryOperator>)
+	{
+	}
 
-	template <class Operator>
+	template <_operator Operator>
 	[[nodiscard]] string ConformalBlock<Operator>::str() const
 	{
 		ostringstream os;
diff --git a/src/bootstrap_equation.cpp b/src/bootstrap_equation.cpp
index 844b07c..eae9fcb 100644
--- a/src/bootstrap_equation.cpp
+++ b/src/bootstrap_equation.cpp
@@ -354,7 +354,7 @@ namespace qboot
 		return prg;
 	}
 
-	Vector<real> read_raw_functional(const fs::path& y_txt)
+	Vector<real> read_raw_functional(const std::filesystem::path& y_txt)
 	{
 		std::ifstream is(y_txt);
 		uint32_t sz, M;
@@ -369,7 +369,7 @@ namespace qboot
 	Matrix<real> BootstrapEquation::apply_functional(const Vector<real>& func, string_view disc_sector) const
 	{
 		auto id = get_id(disc_sector);
-		auto sec = sector(id);
+		const auto& sec = sector(id);
 		assert(sec.type() == SectorType::Discrete);
 		auto mats = make_disc_mat(id);
 		auto sz = sec.size();
@@ -473,7 +473,7 @@ namespace qboot
 			if (sector(id).type_ != SectorType::Continuous) continue;
 			auto sz = sector(id).size();
 			for (const auto& op : sector(id).ops_)
-				spectrums.emplace_back([this, id = id, sec = sec, op = op, sz, &recovered_func, &inv, &event]() {
+				spectrums.emplace_back([this, id = id, sec = sec, &op, sz, &recovered_func, &inv, &event]() {
 					vector<PrimaryOperator> spectrum;
 					auto tag = op.str();
 					tag += " in ";
@@ -576,7 +576,7 @@ namespace qboot
 					});
 			else
 				for (const auto& op : sector(id).ops_)
-					ineqs.emplace_back([this, id = id, sz, op = op, sec = sec, &event] {
+					ineqs.emplace_back([this, id = id, sz, &op, sec = sec, &event] {
 						auto tag = op.str();
 						tag += " in ";
 						tag += sec;
diff --git a/src/polynomial_program.cpp b/src/polynomial_program.cpp
index 4a6fcae..940038f 100644
--- a/src/polynomial_program.cpp
+++ b/src/polynomial_program.cpp
@@ -173,18 +173,18 @@ namespace qboot
 		for (uint32_t j = 0; j < inequality_.size(); ++j)
 			tasks.emplace_back([this, &sdpb, j, M, eq_sz, &event] {
 				_scoped_event scope(std::to_string(j), event);
-				auto ineq = std::move(inequality_[j]);
+				auto ineq = move(inequality_[j].value());
 				// convert ineq to DualConstraint
-				uint32_t sz = ineq->size(), deg = ineq->max_degree(), schur_sz = (deg + 1) * sz * (sz + 1) / 2,
+				uint32_t sz = ineq.size(), deg = ineq.max_degree(), schur_sz = (deg + 1) * sz * (sz + 1) / 2,
 				         d0 = deg / 2, d1 = deg == 0 ? 0 : (deg - 1) / 2;
 				Vector<real> d_c(schur_sz);
 				Matrix<real> d_B(schur_sz, M);
 				Matrix<real> q0(d0 + 1, deg + 1);
 				Matrix<real> q1(d1 + 1, deg + 1);
 				{
-					const auto& xs = ineq->sample_points();
-					const auto& scs = ineq->sample_scalings();
-					auto q = ineq->bilinear_bases();
+					const auto& xs = ineq.sample_points();
+					const auto& scs = ineq.sample_scalings();
+					auto q = ineq.bilinear_bases();
 					for (uint32_t m = 0; m <= d0; ++m)
 						for (uint32_t k = 0; k <= deg; ++k) q0.at(m, k) = q[m].eval(xs[k]) * mp::sqrt(scs[k]);
 					for (uint32_t m = 0; m <= d1; ++m)
@@ -193,13 +193,13 @@ namespace qboot
 				{
 					Vector<Matrix<real>> e_c(deg + 1);
 					Vector<Vector<Matrix<real>>> e_B(N_);
-					for (uint32_t k = 0; k <= deg; ++k) e_c[k] = -ineq->target_eval_with_scale(k);
+					for (uint32_t k = 0; k <= deg; ++k) e_c[k] = -ineq.target_eval_with_scale(k);
 					for (uint32_t n = 0; n < N_; ++n)
 					{
 						e_B[n] = Vector<Matrix<real>>{deg + 1};
-						for (uint32_t k = 0; k <= deg; ++k) e_B[n][k] = -ineq->matrix_eval_with_scale(n, k);
+						for (uint32_t k = 0; k <= deg; ++k) e_B[n][k] = -ineq.matrix_eval_with_scale(n, k);
 					}
-					ineq.reset();
+					move(ineq)._reset();
 					// Tr(A_p Y) + (e_B y)_p = (e_c)_p
 					// convert to Tr(A_p Z) + (d_B z)_p = (d_c)_p
 					uint32_t p = 0;
@@ -247,17 +247,17 @@ namespace qboot
 		for (uint32_t j = 0; j < inequality_.size(); ++j)
 			tasks.emplace_back([this, &sdpb, j, M, eq_sz, &event] {
 				_scoped_event scope(std::to_string(j), event);
-				auto ineq = std::move(inequality_[j]);
-				uint32_t sz = ineq->size();
+				auto ineq = move(inequality_[j].value());
+				uint32_t sz = ineq.size();
 				Matrix<Vector<Polynomial>> mat(sz, sz);
 				for (uint32_t r = 0; r < sz; ++r)
 					for (uint32_t c = 0; c < sz; ++c) mat.at(r, c) = Vector<Polynomial>(M + 1);
-				auto target = -ineq->target_polynomial();
+				auto target = -ineq.target_polynomial();
 				Vector<Matrix<Polynomial>> new_mat(M);
-				for (uint32_t m = 0; m < M; ++m) new_mat[m] = ineq->matrix_polynomial(free_indices_[m]);
+				for (uint32_t m = 0; m < M; ++m) new_mat[m] = ineq.matrix_polynomial(free_indices_[m]);
 				for (uint32_t e = 0; e < eq_sz; ++e)
 				{
-					auto t = ineq->matrix_polynomial(leading_indices_[e]);
+					auto t = ineq.matrix_polynomial(leading_indices_[e]);
 					for (uint32_t m = 0; m < M; ++m) new_mat[m] -= mul_scalar(equation_[e][free_indices_[m]], t);
 					target += mul_scalar(equation_targets_[e], move(t));
 				}
@@ -268,8 +268,8 @@ namespace qboot
 						for (uint32_t m = 0; m < M; ++m) mat.at(r, c).at(m + 1) = move(new_mat[m].at(r, c));
 					}
 				sdpb.register_constraint(
-				    j, PVM(move(mat), ineq->sample_points(), ineq->sample_scalings(), ineq->bilinear_bases()));
-				ineq.reset();
+				    j, PVM(move(mat), ineq.sample_points(), ineq.sample_scalings(), ineq.bilinear_bases()));
+				move(ineq)._reset();
 			});
 		_parallel_evaluate(tasks, parallel);
 		equation_.clear();
diff --git a/src/sdpb_input.cpp b/src/sdpb_input.cpp
index 4e44f75..a94f1f2 100644
--- a/src/sdpb_input.cpp
+++ b/src/sdpb_input.cpp
@@ -5,7 +5,7 @@
 
 #include "qboot/task_queue.hpp"  // for _parallel_evaluate, _event_base
 
-namespace fs = qboot::fs;
+namespace fs = std::filesystem;
 
 using qboot::algebra::Vector, qboot::algebra::Matrix;
 using qboot::mp::real;
diff --git a/src/xml_input.cpp b/src/xml_input.cpp
index 9e854e5..752ecb7 100644
--- a/src/xml_input.cpp
+++ b/src/xml_input.cpp
@@ -7,7 +7,7 @@
 #include <string_view>  // for string_view
 #include <utility>      // for move
 
-namespace fs = qboot::fs;
+namespace fs = std::filesystem;
 
 using qboot::algebra::Vector, qboot::algebra::Matrix, qboot::algebra::Polynomial;
 using qboot::mp::real;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index d7ee717..441f611 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,9 +1,15 @@
 set(MY_TARGET test)
 add_executable(${MY_TARGET} src/main.cpp)
 target_link_libraries(${MY_TARGET} PRIVATE qboot)
+if(NOT MSVC AND NOT CMAKE_COMPILER_IS_GNUCXX)
+	target_link_libraries(${MY_TARGET} PRIVATE c++abi)
+endif()
 set_default_property(${MY_TARGET} ..)
 
 set(MY_TARGET comp)
 add_executable(${MY_TARGET} src/comp.cpp)
 target_link_libraries(${MY_TARGET} PRIVATE qboot)
+if(NOT MSVC AND NOT CMAKE_COMPILER_IS_GNUCXX)
+	target_link_libraries(${MY_TARGET} PRIVATE c++abi)
+endif()
 set_default_property(${MY_TARGET} ..)
diff --git a/test/src/comp.cpp b/test/src/comp.cpp
index b00404b..d057c85 100644
--- a/test/src/comp.cpp
+++ b/test/src/comp.cpp
@@ -1,5 +1,6 @@
 #include <array>
 #include <cassert>
+#include <filesystem>
 #include <fstream>
 #include <future>
 #include <iostream>
@@ -11,11 +12,10 @@
 #include <vector>
 
 #include "qboot/mp/real.hpp"
-#include "qboot/my_filesystem.hpp"
 #include "qboot/task_queue.hpp"
 
 namespace mp = qboot::mp;
-namespace fs = qboot::fs;
+namespace fs = std::filesystem;
 using fs::path;
 using mp::real, mp::rational, mp::integer;
 using std::ostringstream, std::ifstream, std::endl, std::cout;
diff --git a/test/src/main.cpp b/test/src/main.cpp
index 077b1a4..2277606 100644
--- a/test/src/main.cpp
+++ b/test/src/main.cpp
@@ -17,7 +17,6 @@
 #include "qboot/hor_formula.hpp"
 #include "qboot/hor_recursion.hpp"
 #include "qboot/mp/real.hpp"
-#include "qboot/my_filesystem.hpp"
 #include "qboot/polynomial_program.hpp"
 #include "qboot/primary_op.hpp"
 #include "qboot/sdpb_input.hpp"
@@ -29,7 +28,7 @@ using qboot::algebra::Vector, qboot::algebra::Matrix, qboot::algebra::Polynomial
 using qboot::mp::real, qboot::mp::rational, qboot::mp::parse;
 using std::array, std::unique_ptr, std::cout, std::endl, std::map, std::optional, std::make_unique, std::move,
     std::pair, std::vector, std::string;
-namespace fs = qboot::fs;
+namespace fs = std::filesystem;
 
 template <class T>
 using dict = map<string, T, std::less<>>;