diff --git a/include/boost/decimal.hpp b/include/boost/decimal.hpp
index b746182af..af85af1ef 100644
--- a/include/boost/decimal.hpp
+++ b/include/boost/decimal.hpp
@@ -37,6 +37,7 @@
 #include <boost/decimal/type_traits.hpp>
 #include <boost/decimal/detail/io.hpp>
 #include <boost/decimal/cstdio.hpp>
+#include <boost/decimal/complex.hpp>
 
 #if defined(__clang__) && !defined(__GNUC__)
 #  pragma clang diagnostic pop
diff --git a/include/boost/decimal/complex.hpp b/include/boost/decimal/complex.hpp
new file mode 100644
index 000000000..6c240b390
--- /dev/null
+++ b/include/boost/decimal/complex.hpp
@@ -0,0 +1,394 @@
+// Copyright 2024 Matt Borland
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+
+#ifndef BOOST_DECIMAL_COMPLEX_HPP
+#define BOOST_DECIMAL_COMPLEX_HPP
+
+#include <boost/decimal/decimal32.hpp>
+#include <boost/decimal/decimal64.hpp>
+#include <boost/decimal/decimal128.hpp>
+#include <boost/decimal/cmath.hpp>
+#include <boost/decimal/detail/concepts.hpp>
+#include <boost/decimal/detail/type_traits.hpp>
+
+#ifndef BOOST_DECIMAL_BUILD_MODULE
+#include <type_traits>
+#include <complex>
+#include <limits>
+#endif
+
+namespace boost {
+namespace decimal {
+
+template <typename T>
+class complex
+{
+private:
+    T real_;
+    T imag_;
+
+public:
+    constexpr complex() noexcept = default;
+    constexpr explicit complex(T real) noexcept : real_ {real}, imag_ {T{0}} {}
+    constexpr complex(T real, T imag) noexcept : real_ {real}, imag_ {imag} {}
+
+    constexpr T real() const noexcept { return real_; }
+    constexpr T imag() const noexcept { return imag_; }
+
+    constexpr complex operator+() const { return *this; }
+    constexpr complex operator-() const { return {-real_, -imag_}; }
+
+    friend constexpr complex operator+(const complex& lhs, const complex& rhs) noexcept
+    {
+        return {lhs.real_ + rhs.real_, lhs.imag_ + rhs.imag_};
+    }
+
+    friend constexpr complex operator+(const complex& lhs, const T& rhs) noexcept
+    {
+        return {lhs.real_ + rhs, lhs.imag_};
+    }
+
+    friend constexpr complex operator+(const T& lhs, const complex& rhs) noexcept
+    {
+        return {lhs + rhs.real_, rhs.imag_};
+    }
+
+    friend constexpr complex operator-(const complex& lhs, const complex& rhs) noexcept
+    {
+        return {lhs.real_ - rhs.real_, lhs.imag_ - rhs.imag_};
+    }
+
+    friend constexpr complex operator-(const complex& lhs, const T& rhs) noexcept
+    {
+        return {lhs.real_ - rhs, lhs.imag_};
+    }
+
+    friend constexpr complex operator-(const T& lhs, const complex& rhs) noexcept
+    {
+        return {lhs - rhs.real_, -rhs.imag_};
+    }
+
+    friend constexpr complex operator*(const complex& lhs, const complex& rhs) noexcept
+    {
+        return {lhs.real_ * rhs.real_ - lhs.imag_ * rhs.imag_, lhs.imag_ * rhs.real_ + lhs.real_ * rhs.imag_};
+    }
+
+    friend constexpr complex operator*(const complex& lhs, const T& rhs) noexcept
+    {
+        return {lhs.real_ * rhs, lhs.imag_ * rhs};
+    }
+
+    friend constexpr complex operator*(const T& lhs, const complex& rhs) noexcept
+    {
+        return {lhs * rhs.real_, lhs * rhs.imag_};
+    }
+
+    friend constexpr complex operator/(const complex& lhs, const complex& rhs) noexcept
+    {
+        const T divisor = rhs.real_ * rhs.real_ + rhs.imag_ * rhs.imag_;
+        const T real_part = (lhs.real_ * rhs.real_ + lhs.imag_ * rhs.imag_) / divisor;
+        const T imag_part = (lhs.imag_ * rhs.real_ - lhs.real_ * rhs.imag_) / divisor;
+        return {real_part, imag_part};
+    }
+
+    friend constexpr complex operator/(const complex& lhs, const T& rhs) noexcept
+    {
+        const T divisor = rhs * rhs;
+        const T real_part = (lhs.real_ * rhs) / divisor;
+        const T imag_part = (lhs.imag_ * rhs) / divisor;
+        return {real_part, imag_part};
+    }
+
+    friend constexpr complex operator/(const T& lhs, const complex& rhs) noexcept
+    {
+        const T divisor = rhs.real_ * rhs.real_ + rhs.imag_ * rhs.imag_;
+        const T real_part = (lhs * rhs.real_) / divisor;
+        const T imag_part = -(lhs.real_ * rhs.imag_) / divisor;
+        return {real_part, imag_part};
+    }
+
+    constexpr complex& operator+=(const complex& rhs) noexcept
+    {
+        *this = *this + rhs;
+        return *this;
+    }
+
+    constexpr complex& operator+=(const T& rhs) noexcept
+    {
+        *this = *this + rhs;
+        return *this;
+    }
+
+    constexpr complex& operator-=(const complex& rhs) noexcept
+    {
+        *this = *this - rhs;
+        return *this;
+    }
+
+    constexpr complex& operator-=(const T& rhs) noexcept
+    {
+        *this = *this - rhs;
+        return *this;
+    }
+
+    constexpr complex& operator*=(const complex& rhs) noexcept
+    {
+        *this = *this * rhs;
+        return *this;
+    }
+
+    constexpr complex& operator*=(const T& rhs) noexcept
+    {
+        *this = *this * rhs;
+        return *this;
+    }
+
+    constexpr complex& operator/=(const complex& rhs) noexcept
+    {
+        *this = *this / rhs;
+        return *this;
+    }
+
+    constexpr complex& operator/=(const T& rhs) noexcept
+    {
+        *this = *this / rhs;
+        return *this;
+    }
+
+    constexpr bool operator==(const complex& rhs) const noexcept
+    {
+        return real_ == rhs.real_ && imag_ == rhs.imag_;
+    }
+
+    constexpr bool operator!=(const complex& rhs) const noexcept
+    {
+        return !(*this == rhs);
+    }
+
+    constexpr bool operator==(const T& rhs) const noexcept
+    {
+        return real_ == rhs && imag_ == T{0};
+    }
+
+    constexpr bool operator!=(const T& rhs) const noexcept
+    {
+        return !(*this == rhs);
+    }
+
+    template <typename CharT, typename Traits>
+    friend std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os, const complex& 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, complex& z)
+    {
+        CharT ch {};
+        T real = T{0};
+        T imag = T{0};
+
+        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.real_ = real;
+            z.imag_ = imag;
+        }
+
+        return is;
+    }
+};
+
+// Polar is specialized for each type since libc++ uses just template<typename T> which leads
+// to a compiler error
+
+constexpr complex<decimal32> polar(const decimal32& rho, const decimal32& theta) noexcept
+{
+    return {rho * cos(theta), rho * sin(theta)};
+}
+
+constexpr complex<decimal64> polar(const decimal64& rho, const decimal64& theta) noexcept
+{
+    return {rho * cos(theta), rho * sin(theta)};
+}
+
+constexpr complex<decimal128> polar(const decimal128& rho, const decimal128& theta) noexcept
+{
+    return {rho * cos(theta), rho * sin(theta)};
+}
+
+template <typename T>
+constexpr T real(const complex<T>& z) noexcept
+{
+    return z.real();
+}
+
+template <typename T>
+constexpr T imag(const complex<T>& z) noexcept
+{
+    return z.imag();
+}
+
+template <typename T>
+constexpr T abs(const complex<T>& z) noexcept
+{
+    return hypot(z.real(), z.imag());
+}
+
+template <typename T>
+constexpr T arg(const complex<T>& z) noexcept
+{
+    return atan2(z.imag(), z.real());
+}
+
+template <typename T>
+constexpr T norm(const complex<T>& z) noexcept
+{
+    return z.real() * z.real() + z.imag() * z.imag();
+}
+
+template <typename T>
+constexpr complex<T> conj(const complex<T>& z) noexcept
+{
+    return {z.real(), -z.imag()};
+}
+
+template <typename T>
+constexpr complex<T> proj(const complex<T>& z) noexcept
+{
+    if (isinf(z.real()) || isinf(z.imag()))
+    {
+        return {std::numeric_limits<T>::infinity(), copysign(T{0}, z.imag())};
+    }
+
+    return z;
+}
+
+template <typename T>
+constexpr complex<T> exp(const complex<T>& z) noexcept
+{
+    return polar(exp(z.real()), z.imag());
+}
+
+template <typename T>
+constexpr complex<T> log(const complex<T>& z) noexcept
+{
+    return {log(abs(z)), arg(z)};
+}
+
+template <typename T>
+constexpr complex<T> log10(const complex<T>& z) noexcept
+{
+    return log(z) / log(T{10});
+}
+
+template <typename T>
+constexpr complex<T> pow(const complex<T>& x, const complex<T>& y) noexcept
+{
+    return exp(y * log(x));
+}
+
+template <typename T>
+constexpr complex<T> pow(const T& x, const complex<T>& y) noexcept
+{
+    const complex<T> new_x {x};
+    return exp(y * log(new_x));
+}
+
+template <typename T>
+constexpr complex<T> pow(const complex<T>& x, const T& y) noexcept
+{
+    const complex<T> new_y {y};
+    return exp(new_y * log(x));
+}
+
+template <typename T>
+constexpr complex<T> sqrt(const complex<T>& z) noexcept
+{
+    return polar(sqrt(abs(z)), arg(z) / 2);
+}
+
+template <typename T>
+constexpr complex<T> sinh(const complex<T>& z) noexcept
+{
+    return {sinh(z.real()) * cos(z.imag()), cosh(z.real()) * sin(z.imag())};
+}
+
+template <typename T>
+constexpr complex<T> cosh(const complex<T>& z) noexcept
+{
+    return {cosh(z.real()) * cos(z.imag()), sinh(z.real()) * sin(z.imag())};
+}
+
+template <typename T>
+constexpr complex<T> tanh(const complex<T>& z) noexcept
+{
+    return {sinh(z) / cosh(z)};
+}
+
+} //namespace decimal
+} //namespace boost
+
+namespace std {
+
+template <>
+class complex<boost::decimal::decimal32> : public boost::decimal::complex<boost::decimal::decimal32>
+{
+public:
+    using boost::decimal::complex<boost::decimal::decimal32>::complex;
+};
+
+template <>
+class complex<boost::decimal::decimal64> : public boost::decimal::complex<boost::decimal::decimal64>
+{
+public:
+    using boost::decimal::complex<boost::decimal::decimal64>::complex;
+};
+
+template <>
+class complex<boost::decimal::decimal128> : public boost::decimal::complex<boost::decimal::decimal128>
+{
+    using boost::decimal::complex<boost::decimal::decimal128>::complex;
+};
+
+} // namespace std
+
+#endif // BOOST_DECIMAL_COMPLEX_HPP
diff --git a/test/Jamfile b/test/Jamfile
index cf017d74b..5416bb9d6 100644
--- a/test/Jamfile
+++ b/test/Jamfile
@@ -75,6 +75,7 @@ run test_big_uints.cpp ;
 run test_boost_math_univariate_stats.cpp ;
 run test_cbrt.cpp ;
 run test_cmath.cpp ;
+run test_complex.cpp ;
 run test_constants.cpp ;
 run test_cosh.cpp ;
 run test_decimal32.cpp ;
diff --git a/test/test_complex.cpp b/test/test_complex.cpp
new file mode 100644
index 000000000..ef5759a29
--- /dev/null
+++ b/test/test_complex.cpp
@@ -0,0 +1,708 @@
+// Copyright 2024 Matt Borland
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/decimal.hpp>
+
+// Propogates up from boost.math
+#define _SILENCE_CXX23_DENORM_DEPRECATION_WARNING
+
+#if defined(__clang__)
+#  pragma clang diagnostic push
+#  pragma clang diagnostic ignored "-Wold-style-cast"
+#  pragma clang diagnostic ignored "-Wundef"
+#  pragma clang diagnostic ignored "-Wconversion"
+#  pragma clang diagnostic ignored "-Wsign-conversion"
+#  pragma clang diagnostic ignored "-Wfloat-equal"
+#  pragma clang diagnostic ignored "-Wfloat-conversion"
+
+#  if (__clang_major__ >= 10 && !defined(__APPLE__)) || __clang_major__ >= 13
+#    pragma clang diagnostic ignored "-Wdeprecated-copy"
+#  endif
+
+#elif defined(__GNUC__)
+#  pragma GCC diagnostic push
+#  pragma GCC diagnostic ignored "-Wold-style-cast"
+#  pragma GCC diagnostic ignored "-Wundef"
+#  pragma GCC diagnostic ignored "-Wconversion"
+#  pragma GCC diagnostic ignored "-Wsign-conversion"
+#  pragma GCC diagnostic ignored "-Wfloat-equal"
+#  pragma GCC diagnostic ignored "-Wfloat-conversion"
+#endif
+
+#include <boost/math/constants/constants.hpp>
+#include <boost/core/lightweight_test.hpp>
+#include <limits>
+#include <complex>
+#include <iostream>
+#include <type_traits>
+#include <cmath>
+
+using namespace boost::decimal;
+
+template <typename T>
+bool test_equal(T lhs, T rhs, int tol = 15) 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>
+void test_construction()
+{
+    using std::complex;
+    using std::polar;
+    using complex_scalar = decltype(polar(T(), T()));
+    std::cerr << typeid(complex_scalar).name() << std::endl;
+
+    complex_scalar v {};
+    BOOST_TEST(test_equal(v.real(), T{0}));
+    BOOST_TEST(test_equal(v.imag(), T{0}));
+
+    complex_scalar v1 {T{1}};
+    BOOST_TEST(test_equal(v1.real(), T{1}));
+    BOOST_TEST(test_equal(v1.imag(), T{0}));
+
+    complex_scalar v2 {T{2}, T{2}};
+    BOOST_TEST(test_equal(v2.real(), T{2}));
+    BOOST_TEST(test_equal(v2.imag(), T{2}));
+}
+
+template <typename T>
+void test_unary_operators()
+{
+    using std::complex;
+    using std::polar;
+    using complex_scalar = decltype(polar(T(), T()));
+
+    const complex_scalar val {T{2}, T{-2}};
+    complex_scalar pos_val = +val;
+    BOOST_TEST(test_equal(val.real(), pos_val.real()));
+    BOOST_TEST(test_equal(val.imag(), pos_val.imag()));
+
+    complex_scalar neg_val = -val;
+    BOOST_TEST(test_equal(neg_val.real(), T{-2}));
+    BOOST_TEST(test_equal(neg_val.imag(), T{2}));
+}
+
+template <typename T>
+void test_addition()
+{
+    using std::complex;
+    using std::polar;
+    using complex_scalar = decltype(polar(T(), 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(test_equal(res_1.real(), T{3}));
+    BOOST_TEST(test_equal(res_1.imag(), T{3}));
+}
+
+template <typename T>
+void test_subtraction()
+{
+    using std::complex;
+    using std::polar;
+    using complex_scalar = decltype(polar(T(), 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(test_equal(res_1.real(), T{-1}));
+    BOOST_TEST(test_equal(res_1.imag(), T{-1}));
+}
+
+template <typename T>
+void test_multiplication()
+{
+    using std::complex;
+    using std::polar;
+    using complex_scalar = decltype(polar(T(), 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(test_equal(res_1.real(), T{-12}));
+    BOOST_TEST(test_equal(res_1.imag(), T{0}));
+}
+
+template <typename T>
+void test_division()
+{
+    using std::complex;
+    using std::polar;
+    using complex_scalar = decltype(polar(T(), 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(test_equal(res_1.real(), T{2}));
+    BOOST_TEST(test_equal(res_1.imag(), T{-1}));
+}
+
+template <typename T>
+void test_equality()
+{
+    using std::complex;
+    using std::polar;
+    using complex_scalar = decltype(polar(T(), 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::polar;
+    using std::real;
+    using complex_scalar = decltype(polar(T(), T()));
+
+    complex_scalar lhs {T{2}, T{-1}};
+
+    BOOST_TEST(test_equal(real(lhs), lhs.real()));
+    BOOST_TEST(test_equal(imag(lhs), lhs.imag()));
+}
+
+template <typename T>
+void test_abs()
+{
+    using std::complex;
+    using std::polar;
+    using std::abs;
+    using complex_scalar = decltype(polar(T(), 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, std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
+void test_arg()
+{
+    using std::complex;
+    using std::polar;
+    using std::arg;
+    using boost::math::constants::pi;
+    using boost::math::constants::half_pi;
+    using complex_scalar = decltype(polar(T(), 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, std::enable_if_t<!std::is_floating_point<T>::value, bool> = true>
+void test_arg()
+{
+    using std::complex;
+    using std::polar;
+    using std::arg;
+    using complex_scalar = decltype(polar(T(), 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}}), numbers::pi_v<T> / 2));
+    BOOST_TEST(test_equal(arg(complex_scalar{T{-1}, T{0}}), numbers::pi_v<T>));
+}
+
+
+template <typename T>
+void test_norm()
+{
+    using std::complex;
+    using std::polar;
+    using std::norm;
+    using complex_scalar = decltype(polar(T(), 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::polar;
+    using std::conj;
+    using complex_scalar = decltype(polar(T(), 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::polar;
+    using std::proj;
+    using complex_scalar = decltype(polar(T(), 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);
+}
+
+template <typename T, std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
+void test_exp()
+{
+    using std::complex;
+    using std::polar;
+    using std::exp;
+    using boost::math::constants::pi;
+    using complex_scalar = decltype(polar(T(), 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, std::enable_if_t<!std::is_floating_point<T>::value, bool> = true>
+void test_exp()
+{
+    using std::complex;
+    using std::polar;
+    using std::exp;
+    using complex_scalar = decltype(polar(T(), T()));
+
+    complex_scalar lhs {T{0}, numbers::pi_v<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, std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
+void test_log()
+{
+    using std::complex;
+    using std::polar;
+    using std::log;
+    using boost::math::constants::half_pi;
+    using boost::math::constants::pi;
+    using complex_scalar = decltype(polar(T(), 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()));
+
+    // 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, std::enable_if_t<!std::is_floating_point<T>::value, bool> = true>
+void test_log()
+{
+    using std::complex;
+    using std::polar;
+    using std::log;
+    using boost::math::constants::half_pi;
+    using boost::math::constants::pi;
+    using complex_scalar = decltype(polar(T(), T()));
+
+    complex_scalar lhs {T{0}, T{1}};
+    lhs = log(lhs);
+    complex_scalar rhs {T{0}, numbers::pi_v<T> / 2};
+
+    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}, numbers::pi_v<T>};
+
+    BOOST_TEST(test_equal(lhs.real(), rhs.real()));
+    BOOST_TEST(test_equal(lhs.imag(), rhs.imag()));
+
+    // Other side of the cut line
+    lhs = {T {-1}, -T {0}};
+    lhs = log(lhs);
+    rhs = {T {0}, -numbers::pi_v<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 std::polar;
+    using complex_scalar = decltype(polar(T(), T()));
+
+    complex_scalar lhs_1 {T{1}, T{1}};
+    T rhs_1 {T{2}};
+    complex_scalar res_1 = lhs_1 + rhs_1;
+
+    BOOST_TEST(test_equal(res_1.real(), T{3}));
+    BOOST_TEST(test_equal(res_1.imag(), T{1}));
+}
+
+template <typename T>
+void test_scalar_subtraction()
+{
+    using std::complex;
+    using std::polar;
+    using complex_scalar = decltype(polar(T(), T()));
+
+    complex_scalar lhs_1 {T{1}, T{1}};
+    T rhs_1 {T{2}};
+    complex_scalar res_1 = lhs_1 - rhs_1;
+
+    BOOST_TEST(test_equal(res_1.real(), T{-1}));
+    BOOST_TEST(test_equal(res_1.imag(), T{1}));
+}
+
+template <typename T>
+void test_scalar_multiplication()
+{
+    using std::complex;
+    using std::polar;
+    using complex_scalar = decltype(polar(T(), T()));
+
+    complex_scalar lhs_1 {T{3}, T{2}};
+    T rhs_1 {T{2}};
+    complex_scalar res_1 = lhs_1 * rhs_1;
+
+    BOOST_TEST(test_equal(res_1.real(), T{6}));
+    BOOST_TEST(test_equal(res_1.imag(), T{4}));
+}
+
+template <typename T>
+void test_scalar_division()
+{
+    using std::complex;
+    using std::polar;
+    using complex_scalar = decltype(polar(T(), T()));
+
+    complex_scalar lhs_1 {T{4}, T{2}};
+    T rhs_1 {T{2}};
+    complex_scalar res_1 = lhs_1 / rhs_1;
+
+    BOOST_TEST(test_equal(res_1.real(), T{2}));
+    BOOST_TEST(test_equal(res_1.imag(), T{1}));
+}
+
+/*
+template <typename T>
+void test_log10()
+{
+    using std::complex;
+    using std::polar;
+    using std::log10;
+    using complex_scalar = decltype(polar(T(), 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::polar;
+    using std::pow;
+    using complex_scalar = decltype(polar(T(), 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)
+    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::polar;
+    using std::sqrt;
+    using complex_scalar = decltype(polar(T(), 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
+    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::polar;
+    using std::sinh;
+    using complex_scalar = decltype(polar(T(), 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::polar;
+    using std::sinh;
+    using complex_scalar = decltype(polar(T(), 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::polar;
+    using std::sinh;
+    using complex_scalar = decltype(polar(T(), 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<decimal32>();
+    test_construction<decimal64>();
+
+    test_unary_operators<float>();
+    test_unary_operators<double>();
+    test_unary_operators<decimal32>();
+    test_unary_operators<decimal64>();
+
+    test_addition<float>();
+    test_addition<double>();
+    test_addition<decimal32>();
+    test_addition<decimal64>();
+
+    test_subtraction<float>();
+    test_subtraction<double>();
+    test_subtraction<decimal32>();
+    test_subtraction<decimal64>();
+
+    test_multiplication<float>();
+    test_multiplication<double>();
+    test_multiplication<decimal32>();
+    test_multiplication<decimal64>();
+
+    test_division<float>();
+    test_division<double>();
+    test_division<decimal32>();
+    test_division<decimal64>();
+
+    test_equality<float>();
+    test_equality<double>();
+    test_equality<decimal32>();
+    test_equality<decimal64>();
+
+    test_non_member_real_imag<float>();
+    test_non_member_real_imag<double>();
+    test_non_member_real_imag<decimal32>();
+    test_non_member_real_imag<decimal64>();
+
+    test_abs<float>();
+    test_abs<double>();
+    test_abs<decimal32>();
+    test_abs<decimal64>();
+
+    test_arg<float>();
+    test_arg<double>();
+    test_arg<decimal32>();
+    test_arg<decimal64>();
+
+    test_norm<float>();
+    test_norm<double>();
+    test_norm<decimal32>();
+    test_norm<decimal64>();
+
+    test_conj<float>();
+    test_conj<double>();
+    test_conj<decimal32>();
+    test_conj<decimal64>();
+
+    test_proj<float>();
+    test_proj<double>();
+    test_proj<decimal32>();
+    test_proj<decimal64>();
+
+    test_exp<float>();
+    test_exp<double>();
+    test_exp<decimal32>();
+    test_exp<decimal64>();
+
+    /*
+    test_log<float>();
+    test_log<double>();
+    test_log<decimal32>();
+    test_log<decimal64>();
+    */
+
+    test_scalar_addition<float>();
+    test_scalar_addition<double>();
+    test_scalar_addition<decimal32>();
+    test_scalar_addition<decimal64>();
+
+    test_scalar_subtraction<float>();
+    test_scalar_subtraction<double>();
+    test_scalar_subtraction<decimal32>();
+    test_scalar_subtraction<decimal64>();
+
+    test_scalar_multiplication<float>();
+    test_scalar_multiplication<double>();
+    test_scalar_multiplication<decimal32>();
+    test_scalar_multiplication<decimal64>();
+
+    test_scalar_division<float>();
+    test_scalar_division<double>();
+    test_scalar_division<decimal32>();
+    test_scalar_division<decimal64>();
+
+    /*
+    test_log10<float>();
+    test_log10<double>();
+    test_log10<decimal32>();
+    test_log10<decimal64>();
+
+    test_pow<float>();
+    test_pow<double>();
+    test_pow<decimal32>();
+    test_pow<decimal64>();
+    */
+
+    test_sqrt<float>();
+    test_sqrt<double>();
+    test_sqrt<decimal32>();
+    test_sqrt<decimal64>();
+
+    test_sinh<float>();
+    test_sinh<double>();
+    test_sinh<decimal32>();
+    test_sinh<decimal64>();
+
+    test_cosh<float>();
+    test_cosh<double>();
+    test_cosh<decimal32>();
+    test_cosh<decimal64>();
+
+    test_tanh<float>();
+    test_tanh<double>();
+    test_tanh<decimal32>();
+    test_tanh<decimal64>();
+
+    return boost::report_errors();
+}