diff --git a/.github/workflows/multiprecision.yml b/.github/workflows/multiprecision.yml
index 8e2e7726f..b316d9278 100644
--- a/.github/workflows/multiprecision.yml
+++ b/.github/workflows/multiprecision.yml
@@ -435,49 +435,6 @@ jobs:
       - name: b2 config
         run: cat bin.v2/config.log
         working-directory: ../boost-root/
-  windows_msvc_14_0:
-    runs-on: windows-2019
-    defaults:
-      run:
-        shell: cmd
-    env:
-      ARGS: toolset=${{ matrix.toolset }} address-model=64 cxxstd=${{ matrix.standard }}
-    strategy:
-      fail-fast: false
-      matrix:
-        toolset: [ msvc-14.0 ]
-        standard: [ 14, latest ]
-        suite: [ github_ci_block_1, github_ci_block_2 ]
-    steps:
-      - uses: actions/checkout@v4
-        with:
-          fetch-depth: '0'
-      - name: Checkout main boost
-        run: git clone -b develop --depth 1 https://github.com/boostorg/boost.git ../boost-root
-      - name: Update tools/boostdep
-        run: git submodule update --init tools/boostdep
-        working-directory: ../boost-root
-      - name: Copy files
-        run: xcopy /s /e /q %GITHUB_WORKSPACE% libs\multiprecision
-        working-directory: ../boost-root
-      - name: Install deps
-        run: python tools/boostdep/depinst/depinst.py multiprecision
-        working-directory: ../boost-root
-      - name: Bootstrap
-        run: bootstrap
-        working-directory: ../boost-root
-      - name: Generate headers
-        run: b2 headers
-        working-directory: ../boost-root
-      - name: Config Info
-        run: ..\..\..\b2 print_config_info print_math_info %ARGS%
-        working-directory: ../boost-root/libs/config/test
-      - name: Test
-        run: ..\..\..\b2 -j2 --hash %ARGS% define=CI_SUPPRESS_KNOWN_ISSUES ${{ matrix.suite }}
-        working-directory: ../boost-root/libs/multiprecision/test
-      - name: b2 config
-        run: cat bin.v2/config.log
-        working-directory: ../boost-root/
   windows_msvc_14_2:
     runs-on: windows-2019
     defaults:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d6cdd2861..d7d7dee25 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -38,6 +38,12 @@ else()
 
 endif()
 
+if(BUILD_TESTING AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/test/CMakeLists.txt")
+
+  add_subdirectory(test)
+
+endif()
+
 # Only enable tests when we're the root project
 if(CMAKE_SOURCE_DIR STREQUAL CMAKE_CURRENT_SOURCE_DIR)
 
diff --git a/include/boost/multiprecision/complex.hpp b/include/boost/multiprecision/complex.hpp
new file mode 100644
index 000000000..cadf53d9b
--- /dev/null
+++ b/include/boost/multiprecision/complex.hpp
@@ -0,0 +1,408 @@
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2024 Matt Borland. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MP_COMPLEX_HPP
+#define BOOST_MP_COMPLEX_HPP
+
+#include <boost/multiprecision/detail/assert.hpp>
+#include <boost/multiprecision/traits/is_backend.hpp>
+#include <boost/multiprecision/cpp_complex.hpp>
+#include <boost/multiprecision/fwd.hpp>
+#include <utility>
+#include <ostream>
+#include <istream>
+#include <complex>
+#include <cmath>
+
+namespace std {
+
+template <typename T, boost::multiprecision::expression_template_option ET>
+class complex<boost::multiprecision::number<T, ET>>
+{
+private:
+    using float_type = boost::multiprecision::number<T, ET>;
+    using self_type = complex<boost::multiprecision::number<T, ET>>;
+    using complex_type = decltype(boost::multiprecision::polar(std::declval<boost::multiprecision::number<T, ET>>(), std::declval<boost::multiprecision::number<T, ET>>()));
+
+    complex_type m_data_;
+
+public:
+    complex() = default;
+    explicit complex(float_type real) : m_data_ {real} {}
+    complex(float_type real, float_type imag) : m_data_ {real, imag} {}
+    complex(const complex_type& data) : m_data_ {data} {}
+    complex(complex_type&& data) : m_data_ {data} {}
+
+    complex_type& data() { return m_data_; }
+    const complex_type& data() const { return m_data_; }
+
+    float_type real() const { return m_data_.real(); }
+    float_type imag() const { return m_data_.imag(); }
+
+    self_type operator+() const { return *this; }
+    self_type operator-() const { return {-m_data_.real(), -m_data_.imag()}; }
+
+    friend self_type operator+(const self_type& lhs, const self_type& rhs)
+    {
+        return {lhs.m_data_ + rhs.m_data_};
+    }
+
+    friend self_type operator+(const self_type& lhs, const float_type& rhs)
+    {
+        return {lhs.m_data_ + rhs};
+    }
+
+    friend self_type operator+(const float_type& lhs, const self_type& rhs)
+    {
+        return {lhs + rhs.m_data_};
+    }
+
+    friend self_type operator-(const self_type& lhs, const self_type& rhs)
+    {
+        return {lhs.m_data_ - rhs.m_data_};
+    }
+
+    friend self_type operator-(const self_type& lhs, const float_type& rhs)
+    {
+        return {lhs.m_data_ - rhs};
+    }
+
+    friend self_type operator-(const float_type& lhs, const self_type& rhs)
+    {
+        return {lhs - rhs.m_data_};
+    }
+
+    friend self_type operator*(const self_type& lhs, const self_type& rhs)
+    {
+        return {lhs.m_data_ * rhs.m_data_};
+    }
+
+    friend self_type operator*(const self_type& lhs, const float_type& rhs)
+    {
+        return {lhs.m_data_ * rhs};
+    }
+
+    friend self_type operator*(const float_type& lhs, const self_type& rhs)
+    {
+        return {lhs * rhs.m_data_};
+    }
+
+    friend self_type operator/(const self_type& lhs, const self_type& rhs)
+    {
+        return {lhs.m_data_ / rhs.m_data_};
+    }
+
+    friend self_type operator/(const self_type& lhs, const float_type& rhs)
+    {
+        return {lhs.m_data_ / rhs};
+    }
+
+    friend self_type operator/(const float_type& lhs, const self_type& rhs)
+    {
+        return {lhs / rhs.m_data_};
+    }
+
+    self_type& operator+=(const self_type& rhs)
+    {
+        *this = *this + rhs;
+        return *this;
+    }
+
+    self_type& operator+=(const float_type& rhs)
+    {
+        *this = *this + rhs;
+        return *this;
+    }
+
+    self_type& operator-=(const self_type& rhs)
+    {
+        *this = *this - rhs;
+        return *this;
+    }
+
+    self_type& operator-=(const float_type& rhs)
+    {
+        *this = *this - rhs;
+        return *this;
+    }
+
+    self_type& operator*=(const self_type& rhs)
+    {
+        *this = *this * rhs;
+        return *this;
+    }
+
+    self_type& operator*=(const float_type& rhs)
+    {
+        *this = *this * rhs;
+        return *this;
+    }
+
+    self_type& operator/=(const self_type& rhs)
+    {
+        *this = *this / rhs;
+        return *this;
+    }
+
+    self_type& operator/=(const float_type& rhs)
+    {
+        *this = *this / rhs;
+        return *this;
+    }
+
+    bool operator==(const self_type& rhs) const
+    {
+        return this->m_data_ == rhs.m_data_;
+    }
+
+    bool operator!=(const self_type& rhs) const
+    {
+        return !(*this == rhs);
+    }
+
+    bool operator==(const float_type& rhs) const
+    {
+        return this->real() == rhs && this->imag() == float_type{0};
+    }
+
+    bool operator!=(const float_type& rhs) const
+    {
+        return !(*this == rhs);
+    }
+
+    template <typename CharT, typename Traits>
+    friend std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os, const self_type& z)
+    {
+        std::basic_ostringstream<CharT, Traits> s;
+        s.flags(os.flags());
+        s.imbue(os.getloc());
+        s.precision(os.precision());
+        s << '(' << z.real() << ',' << z.imag() << ')';
+
+        return os << s.str();
+    }
+
+    // Supported formats:
+    // 1) real
+    // 2) (real)
+    // 3) (real, imag)
+    template <typename CharT, typename Traits>
+    friend std::basic_istream<CharT, Traits>& operator>>(std::basic_istream<CharT, Traits>& is, self_type& z)
+    {
+        CharT ch {};
+        float_type real {};
+        float_type imag {};
+
+        is >> std::ws;
+        is.get(ch);
+
+        if (ch == '(')
+        {
+            // Expecting a format like 2 or 3
+            is >> std::ws >> real;
+            is.get(ch);
+            if (ch == ',')
+            {
+                // A comma indicates it's 3
+                is >> std::ws >> imag;
+                is.get(ch); // Should be ')'
+            }
+            else if (ch != ')')
+            {
+                // Syntax error: unexpected character
+                is.setstate(std::ios_base::failbit);
+                return is;
+            }
+        }
+        else
+        {
+            // No parentheses, just a real number from format 1
+            is.putback(ch);
+            is >> real;
+        }
+
+        if (!is.fail())
+        {
+            z = self_type{real, imag};
+        }
+
+        return is;
+    }
+};
+
+template <typename T, boost::multiprecision::expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> polar(const boost::multiprecision::number<T, ET>& a, const boost::multiprecision::number<T, ET>& b)
+{
+   return { boost::multiprecision::polar(a, b) };
+}
+
+} // namespace std
+
+namespace boost {
+namespace multiprecision {
+
+template <typename T, expression_template_option ET>
+inline boost::multiprecision::number<T, ET> real(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return z.real();
+}
+
+template <typename T, expression_template_option ET>
+inline boost::multiprecision::number<T, ET> imag(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return z.imag();
+}
+
+template <typename T, expression_template_option ET>
+inline boost::multiprecision::number<T, ET> abs(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return abs(z.data());
+}
+
+template <typename T, expression_template_option ET>
+inline boost::multiprecision::number<T, ET> arg(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return arg(z.data());
+}
+
+template <typename T, expression_template_option ET>
+inline boost::multiprecision::number<T, ET> norm(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return norm(z.data());
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> conj(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {conj(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> proj(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {proj(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> exp(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {exp(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> log(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {log(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> log10(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {log10(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline  std::complex<boost::multiprecision::number<T, ET>> pow(const std::complex<boost::multiprecision::number<T, ET>>& x,
+                                                               const std::complex<boost::multiprecision::number<T, ET>>& y)
+{
+    return {pow(x.data(), y.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> pow(const boost::multiprecision::number<T, ET>& x,
+                                                              const std::complex<boost::multiprecision::number<T, ET>>& y)
+{
+    return {pow(x, y.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> pow(const std::complex<boost::multiprecision::number<T, ET>>& x,
+                                                              const boost::multiprecision::number<T, ET>& y)
+{
+    return {pow(x.data(), y)};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> sqrt(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {sqrt(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> sin(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {sin(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> cos(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {cos(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> tan(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {tan(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> asin(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {asin(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> acos(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {acos(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> atan(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {atan(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> sinh(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {sinh(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> cosh(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {cosh(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> tanh(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {tanh(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> asinh(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {asinh(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> acosh(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {acosh(z.data())};
+}
+
+template <typename T, expression_template_option ET>
+inline std::complex<boost::multiprecision::number<T, ET>> atanh(const std::complex<boost::multiprecision::number<T, ET>>& z)
+{
+    return {atanh(z.data())};
+}
+
+} // namespace MP
+} // namespace boost
+
+#endif //BOOST_MP_COMPLEX_HPP
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
new file mode 100644
index 000000000..a9ff2f412
--- /dev/null
+++ b/test/CMakeLists.txt
@@ -0,0 +1,9 @@
+# Copyright 2018, 2019 Peter Dimov
+# Copyright 2024 Matt Borland
+# Distributed under the Boost Software License, Version 1.0.
+# See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt
+
+file(GLOB SOURCES "*.cpp")
+add_library(boost_mp_tests STATIC ${SOURCES})
+target_compile_features(boost_mp_tests PRIVATE cxx_std_17)
+target_link_libraries(boost_mp_tests PUBLIC Boost::multiprecision Boost::assert Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::math Boost::random)
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 3e490a13c..e4d5b0893 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -223,6 +223,8 @@ test-suite arithmetic_tests :
    [ run test_complex_signed_zero.cpp : : : <define>TEST_CPP_BIN_FLOAT [ requires cxx17_if_constexpr ] : test_complex_signed_zero_cpp_bin_float ]
    [ run test_complex_signed_zero.cpp mpc mpfr gmp : : : <define>TEST_MPC [ requires cxx17_if_constexpr ] [ check-target-builds ../config//has_mpc : : <build>no ] : test_complex_signed_zero_mpc ]
 
+   [ run test_complex_class.cpp ]
+
    [ run test_preserve_source_precision.cpp gmp : : : [ requires cxx17_if_constexpr ] [ check-target-builds ../config//has_gmp : : <build>no ] <define>TEST_MPF : test_preserve_source_precision_gmp ]
    [ run test_preserve_source_precision.cpp gmp mpfr : : : [ requires cxx17_if_constexpr ] [ check-target-builds ../config//has_mpfr : : <build>no ] <define>TEST_MPFR  : test_preserve_source_precision_mpfr ]
    [ run test_preserve_source_precision.cpp gmp mpfr mpc : : : [ requires cxx17_if_constexpr ] [ check-target-builds ../config//has_mpc : : <build>no ] <define>TEST_MPC : test_preserve_source_precision_mpc ]
@@ -1264,6 +1266,7 @@ test-suite misc :
       # Eigen interoperability:
       #
       [ run test_eigen_interop_cpp_int.cpp : : : release [ check-target-builds ../config//has_eigen : : <build>no ] ]
+      [ run boost_math_issue_1131.cpp : : : release [ check-target-builds ../config//has_eigen : : <build>no ] ]
       [ run test_eigen_interop_cpp_dec_float.cpp : : : release [ check-target-builds ../config//has_eigen : : <build>no ] ]
       [ run test_eigen_interop_cpp_dec_float_2.cpp : : : release [ check-target-builds ../config//has_eigen : : <build>no ] ]
       [ run test_eigen_interop_cpp_dec_float_3.cpp : : : release [ check-target-builds ../config//has_eigen : : <build>no ] ]
diff --git a/test/boost_math_issue_1131.cpp b/test/boost_math_issue_1131.cpp
new file mode 100644
index 000000000..57d0eb9d0
--- /dev/null
+++ b/test/boost_math_issue_1131.cpp
@@ -0,0 +1,24 @@
+// Copyright 2024 Nick Thompson
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+
+#include <Eigen/Dense>
+#include <Eigen/Eigenvalues>
+#include <boost/multiprecision/eigen.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <iostream>
+
+using boost::multiprecision::cpp_bin_float_100;
+
+int main() {
+    Eigen::Matrix<cpp_bin_float_100, Eigen::Dynamic, Eigen::Dynamic> A = Eigen::Matrix<cpp_bin_float_100, Eigen::Dynamic, Eigen::Dynamic>::Identity(3,3);
+    Eigen::EigenSolver<decltype(A)> es;
+    es.compute(A, /*computeEigenvectors=*/ false);
+
+    auto eigs = es.eigenvalues();
+    for (std::size_t i = 0; i < eigs.size(); ++i) {
+        std::cout << eigs[i] << "\n";
+    }
+
+    return 0;
+}
diff --git a/test/test_complex_class.cpp b/test/test_complex_class.cpp
new file mode 100644
index 000000000..cdfe5bcfe
--- /dev/null
+++ b/test/test_complex_class.cpp
@@ -0,0 +1,625 @@
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2024 Matt Borland. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/multiprecision/complex.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/math/constants/constants.hpp>
+#include <boost/core/lightweight_test.hpp>
+#include <limits>
+#include <complex>
+#include <iostream>
+#include <cmath>
+
+using namespace boost::multiprecision;
+
+template <typename T>
+bool test_equal(T lhs, T rhs, int tol = 10) noexcept
+{
+    using std::fabs;
+    const bool res = fabs(lhs - rhs) < static_cast<T>(tol) * std::numeric_limits<T>::epsilon();
+
+    if (!res)
+    {
+        // LCOV_EXCL_START
+        std::cerr << "LHS: " << lhs
+                  << "\nRHS: " << rhs
+                  << "\nDist: " << fabs(lhs - rhs) / std::numeric_limits<T>::epsilon() << std::endl;
+        // LCOV_EXCL_STOP
+    }
+
+    return res;
+}
+
+template <typename T, std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
+void test_construction()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+    std::cerr << typeid(complex_scalar).name() << std::endl;
+
+    complex_scalar v {};
+    BOOST_TEST_EQ(v.real(), T{0});
+    BOOST_TEST_EQ(v.imag(), T{0});
+
+    complex_scalar v1 {T{1}};
+    BOOST_TEST_EQ(v1.real(), T{1});
+    BOOST_TEST_EQ(v1.imag(), T{0});
+
+    complex_scalar v2 {T{2}, T{2}};
+    BOOST_TEST_EQ(v2.real(), T{2});
+    BOOST_TEST_EQ(v2.imag(), T{2});
+}
+
+template <typename T, std::enable_if_t<!std::is_floating_point<T>::value, bool> = true>
+void test_construction()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+    std::cerr << typeid(complex_scalar).name() << std::endl;
+
+    complex_scalar v {};
+    BOOST_TEST_EQ(v.real(), T{0});
+    BOOST_TEST_EQ(v.imag(), T{0});
+
+    complex_scalar v1 {T{1}};
+    BOOST_TEST_EQ(v1.real(), T{1});
+    BOOST_TEST_EQ(v1.imag(), T{0});
+
+    complex_scalar v2 {T{2}, T{2}};
+    BOOST_TEST_EQ(v2.real(), T{2});
+    BOOST_TEST_EQ(v2.imag(), T{2});
+
+   complex_scalar v3 = std::polar(T(1), T(2));
+   complex_scalar v4 = boost::multiprecision::polar(T(1), T(2));
+   BOOST_TEST_EQ(v3.real(), v4.real());
+   BOOST_TEST_EQ(v3.imag(), v4.imag());
+}
+
+template <typename T>
+void test_unary_operators()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+
+    const complex_scalar val {T{2}, T{-2}};
+    complex_scalar pos_val = +val;
+    BOOST_TEST_EQ(val.real(), pos_val.real());
+    BOOST_TEST_EQ(val.imag(), pos_val.imag());
+
+    complex_scalar neg_val = -val;
+    BOOST_TEST_EQ(neg_val.real(), T{-2});
+    BOOST_TEST_EQ(neg_val.imag(), T{2});
+}
+
+template <typename T>
+void test_addition()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs_1 {T{1}, T{1}};
+    complex_scalar rhs_1 {T{2}, T{2}};
+    complex_scalar res_1 = lhs_1 + rhs_1;
+
+    BOOST_TEST_EQ(res_1.real(), T{3});
+    BOOST_TEST_EQ(res_1.imag(), T{3});
+}
+
+template <typename T>
+void test_subtraction()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs_1 {T{1}, T{1}};
+    complex_scalar rhs_1 {T{2}, T{2}};
+    complex_scalar res_1 = lhs_1 - rhs_1;
+
+    BOOST_TEST_EQ(res_1.real(), T{-1});
+    BOOST_TEST_EQ(res_1.imag(), T{-1});
+}
+
+template <typename T>
+void test_multiplication()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs_1 {T{-3}, T{3}};
+    complex_scalar rhs_1 {T{2}, T{2}};
+    complex_scalar res_1 = lhs_1 * rhs_1;
+
+    BOOST_TEST_EQ(res_1.real(), T{-12});
+    BOOST_TEST_EQ(res_1.imag(), T{0});
+}
+
+template <typename T>
+void test_division()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs_1 {T{6}, T{2}};
+    complex_scalar rhs_1 {T{2}, T{2}};
+    complex_scalar res_1 = lhs_1 / rhs_1;
+
+    BOOST_TEST_EQ(res_1.real(), T{2});
+    BOOST_TEST_EQ(res_1.imag(), T{-1});
+}
+
+template <typename T>
+void test_equality()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{2}, T{-1}};
+    complex_scalar rhs {T{2}, T{1}};
+
+    BOOST_TEST(lhs != rhs);
+    BOOST_TEST(!(lhs == rhs));
+
+    T scalar_rhs {T{2}};
+
+    BOOST_TEST(lhs != scalar_rhs);
+    BOOST_TEST(!(lhs == scalar_rhs));
+
+    lhs = complex_scalar{T{2}, T{0}};
+    BOOST_TEST(lhs == scalar_rhs);
+    BOOST_TEST(!(lhs != scalar_rhs));
+}
+
+template <typename T>
+void test_non_member_real_imag()
+{
+    using std::complex;
+    using std::real;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{2}, T{-1}};
+
+    BOOST_TEST_EQ(real(lhs), lhs.real());
+    BOOST_TEST_EQ(imag(lhs), lhs.imag());
+}
+
+template <typename T>
+void test_abs()
+{
+    using std::complex;
+    using std::abs;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{1}, T{1}};
+
+    BOOST_TEST(test_equal(static_cast<T>(abs(lhs)), static_cast<T>(sqrt(T{2}))));
+}
+
+template <typename T>
+void test_arg()
+{
+    using std::complex;
+    using std::arg;
+    using boost::math::constants::pi;
+    using boost::math::constants::half_pi;
+    using complex_scalar = std::complex<T>;
+
+    BOOST_TEST(test_equal(arg(complex_scalar{T{1}, T{0}}), T{0}));
+    BOOST_TEST(test_equal(arg(complex_scalar{T{0}, T{0}}), T{0}));
+    BOOST_TEST(test_equal(arg(complex_scalar{T{0}, T{1}}), half_pi<T>()));
+    BOOST_TEST(test_equal(arg(complex_scalar{T{-1}, T{0}}), pi<T>()));
+}
+
+template <typename T>
+void test_norm()
+{
+    using std::complex;
+    using std::norm;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{3}, T{4}};
+
+    BOOST_TEST(test_equal(norm(lhs), T{25}));
+}
+
+template <typename T>
+void test_conj()
+{
+    using std::complex;
+    using std::conj;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{1}, T{1}};
+    complex_scalar rhs {T{1}, T{-1}};
+
+    BOOST_TEST_EQ(conj(lhs), rhs);
+}
+
+template <typename T>
+void test_proj()
+{
+    using std::complex;
+    using std::proj;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{1}, T{1}};
+    BOOST_TEST_EQ(lhs, proj(lhs));
+
+    lhs = complex_scalar{T{std::numeric_limits<T>::infinity()}, T{1}};
+    complex_scalar rhs = complex_scalar{T{std::numeric_limits<T>::infinity()}, T{0}};
+    BOOST_TEST_EQ(proj(lhs), rhs);
+
+    lhs = complex_scalar{T{std::numeric_limits<T>::infinity()}, T{-1}};
+    rhs = complex_scalar{T{std::numeric_limits<T>::infinity()}, T{-0}};
+    BOOST_TEST_EQ(proj(lhs), rhs);
+}
+
+template <typename T>
+void test_exp()
+{
+    using std::complex;
+    using std::exp;
+    using boost::math::constants::pi;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{0}, pi<T>()};
+    lhs = exp(lhs);
+    complex_scalar rhs {T{-1}, T{0}};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+}
+
+template <typename T>
+void test_log()
+{
+    using std::complex;
+    using std::log;
+    using boost::math::constants::half_pi;
+    using boost::math::constants::pi;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{0}, T{1}};
+    lhs = log(lhs);
+    complex_scalar rhs {T{0}, half_pi<T>()};
+
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+
+    lhs = {T{-1}, T{0}};
+    lhs = log(lhs);
+    rhs = {T{0}, pi<T>()};
+
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+
+    BOOST_IF_CONSTEXPR (!std::is_same<T, cpp_dec_float_50>::value)
+    {
+        // Other side of the cut line
+        lhs = {T {-1}, -T {0}};
+        lhs = log(lhs);
+        rhs = {T {0}, -pi<T>()};
+
+        BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+        BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+    }
+}
+
+template <typename T>
+void test_scalar_addition()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs_1 {T{1}, T{1}};
+    T rhs_1 {T{2}};
+    complex_scalar res_1 = lhs_1 + rhs_1;
+
+    BOOST_TEST_EQ(res_1.real(), T{3});
+    BOOST_TEST_EQ(res_1.imag(), T{1});
+}
+
+template <typename T>
+void test_scalar_subtraction()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs_1 {T{1}, T{1}};
+    T rhs_1 {T{2}};
+    complex_scalar res_1 = lhs_1 - rhs_1;
+
+    BOOST_TEST_EQ(res_1.real(), T{-1});
+    BOOST_TEST_EQ(res_1.imag(), T{1});
+}
+
+template <typename T>
+void test_scalar_multiplication()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs_1 {T{3}, T{2}};
+    T rhs_1 {T{2}};
+    complex_scalar res_1 = lhs_1 * rhs_1;
+
+    BOOST_TEST_EQ(res_1.real(), T{6});
+    BOOST_TEST_EQ(res_1.imag(), T{4});
+}
+
+template <typename T>
+void test_scalar_division()
+{
+    using std::complex;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs_1 {T{4}, T{2}};
+    T rhs_1 {T{2}};
+    complex_scalar res_1 = lhs_1 / rhs_1;
+
+    BOOST_TEST_EQ(res_1.real(), T{2});
+    BOOST_TEST_EQ(res_1.imag(), T{1});
+}
+
+template <typename T>
+void test_log10()
+{
+    using std::complex;
+    using std::log10;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{-100}, T{0}};
+    lhs = log10(lhs);
+    complex_scalar rhs {T{2}, static_cast<T>(1.36438)};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+}
+
+template <typename T>
+void test_pow()
+{
+    using std::complex;
+    using std::pow;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{1}, T{2}};
+    lhs = pow(lhs, T{2});
+    complex_scalar rhs {T{-3}, T{4}};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real(), 100));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag(), 100));
+
+    lhs = {T{-1}, T{0}};
+    lhs = pow(lhs, T{1}/2);
+    rhs = {T{0}, T{1}};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+
+    // Check other side of the cut
+    // MSVC 14.0 gets this wrong with float and double
+    #if !defined(_MSC_VER) || (_MSC_VER > 1900)
+    BOOST_IF_CONSTEXPR (!std::is_same<T, cpp_dec_float_50>::value)
+    {
+        lhs = {T {-1}, -T {0}};
+        lhs = pow(lhs, T {1} / 2);
+        rhs = {T {0}, T {-1}};
+        BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+        BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+    }
+    #endif
+}
+
+template <typename T>
+void test_sqrt()
+{
+    using std::complex;
+    using std::sqrt;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{4}, T{0}};
+    lhs = sqrt(lhs);
+    complex_scalar rhs {T{2}, T{0}};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+
+    // Check other side of the cut
+    BOOST_IF_CONSTEXPR (!std::is_same<T, cpp_dec_float_50>::value)
+    {
+        lhs = {T {4}, -T {0}};
+        lhs = sqrt(lhs);
+        rhs = {T {2}, -T {0}};
+        BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+        BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+    }
+}
+
+template <typename T>
+void test_sinh()
+{
+    using std::complex;
+    using std::sinh;
+    using std::sin;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{1}, T{0}};
+    lhs = sinh(lhs);
+    complex_scalar rhs {sinh(T{1}), T{0}};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+
+    lhs = {T{0}, T{1}};
+    lhs = sinh(lhs);
+    rhs = {T{0}, sin(T{1})};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+}
+
+template <typename T>
+void test_cosh()
+{
+    using std::complex;
+    using std::cosh;
+    using std::cos;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{1}, T{0}};
+    lhs = cosh(lhs);
+    complex_scalar rhs {cosh(T{1}), T{0}};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+
+    lhs = {T{0}, T{1}};
+    lhs = cosh(lhs);
+    rhs = {cos(T{1}), T{0}};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+}
+
+template <typename T>
+void test_tanh()
+{
+    using std::complex;
+    using std::tanh;
+    using std::tan;
+    using complex_scalar = std::complex<T>;
+
+    complex_scalar lhs {T{1}, T{0}};
+    lhs = tanh(lhs);
+    complex_scalar rhs {tanh(T{1}), T{0}};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+
+    lhs = {T{0}, T{1}};
+    lhs = tanh(lhs);
+    rhs = {T{0}, tan(T{1})};
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+}
+
+int main()
+{
+    test_construction<float>();
+    test_construction<double>();
+    test_construction<cpp_bin_float_50>();
+    test_construction<cpp_bin_float_100>();
+
+    test_unary_operators<float>();
+    test_unary_operators<double>();
+    test_unary_operators<cpp_bin_float_50>();
+    test_unary_operators<cpp_bin_float_100>();
+
+    test_addition<float>();
+    test_addition<double>();
+    test_addition<cpp_bin_float_50>();
+    test_addition<cpp_bin_float_100>();
+
+    test_subtraction<float>();
+    test_subtraction<double>();
+    test_subtraction<cpp_bin_float_50>();
+    test_subtraction<cpp_bin_float_100>();
+
+    test_multiplication<float>();
+    test_multiplication<double>();
+    test_multiplication<cpp_bin_float_50>();
+    test_multiplication<cpp_bin_float_100>();
+
+    test_division<float>();
+    test_division<double>();
+    test_division<cpp_bin_float_50>();
+    test_division<cpp_bin_float_100>();
+
+    test_equality<float>();
+    test_equality<double>();
+    test_equality<cpp_bin_float_50>();
+    test_equality<cpp_bin_float_100>();
+
+    test_non_member_real_imag<float>();
+    test_non_member_real_imag<double>();
+    test_non_member_real_imag<cpp_bin_float_50>();
+    test_non_member_real_imag<cpp_bin_float_100>();
+
+    test_abs<float>();
+    test_abs<double>();
+    test_abs<cpp_bin_float_50>();
+    test_abs<cpp_bin_float_100>();
+
+    test_arg<float>();
+    test_arg<double>();
+    test_arg<cpp_bin_float_50>();
+    test_arg<cpp_bin_float_100>();
+
+    test_norm<float>();
+    test_norm<double>();
+    test_norm<cpp_bin_float_50>();
+    test_norm<cpp_bin_float_100>();
+
+    test_conj<float>();
+    test_conj<double>();
+    test_conj<cpp_bin_float_50>();
+    test_conj<cpp_bin_float_100>();
+
+    test_proj<float>();
+    test_proj<double>();
+    test_proj<cpp_bin_float_50>();
+    test_proj<cpp_bin_float_100>();
+
+    test_exp<float>();
+    test_exp<double>();
+    test_exp<cpp_bin_float_50>();
+    test_exp<cpp_bin_float_100>();
+
+    test_log<float>();
+    test_log<double>();
+    test_log<cpp_bin_float_50>();
+    test_log<cpp_bin_float_100>();
+
+    test_scalar_addition<float>();
+    test_scalar_addition<double>();
+    test_scalar_addition<cpp_bin_float_50>();
+    test_scalar_addition<cpp_bin_float_100>();
+
+    test_scalar_subtraction<float>();
+    test_scalar_subtraction<double>();
+    test_scalar_subtraction<cpp_bin_float_50>();
+    test_scalar_subtraction<cpp_bin_float_100>();
+
+    test_scalar_multiplication<float>();
+    test_scalar_multiplication<double>();
+    test_scalar_multiplication<cpp_bin_float_50>();
+    test_scalar_multiplication<cpp_bin_float_100>();
+
+    test_scalar_division<float>();
+    test_scalar_division<double>();
+    test_scalar_division<cpp_bin_float_50>();
+    test_scalar_division<cpp_bin_float_100>();
+
+    test_log10<float>();
+    test_log10<double>();
+    test_log10<cpp_bin_float_50>();
+    test_log10<cpp_bin_float_100>();
+
+    test_pow<float>();
+    test_pow<double>();
+    test_pow<cpp_bin_float_50>();
+    test_pow<cpp_bin_float_100>();
+
+    test_sqrt<float>();
+    test_sqrt<double>();
+    test_sqrt<cpp_bin_float_50>();
+    test_sqrt<cpp_bin_float_100>();
+
+    test_sinh<float>();
+    test_sinh<double>();
+    test_sinh<cpp_bin_float_50>();
+    test_sinh<cpp_bin_float_100>();
+
+    test_cosh<float>();
+    test_cosh<double>();
+    test_cosh<cpp_bin_float_50>();
+    test_cosh<cpp_bin_float_100>();
+
+    test_tanh<float>();
+    test_tanh<double>();
+    test_tanh<cpp_bin_float_50>();
+    test_tanh<cpp_bin_float_100>();
+
+    return boost::report_errors();
+}