cp-algo/util/complex.hpp

• View this file on GitHub
• Last update: 2025-12-17 02:23:31+01:00
• Include: #include "cp-algo/util/complex.hpp"

Depends on

• cp-algo/util/simd.hpp

Required by

• cp-algo/geometry/closest_pair.hpp
• cp-algo/geometry/convex_hull.hpp
• cp-algo/geometry/point.hpp
• cp-algo/linalg/frobenius.hpp
• cp-algo/math/bigint.hpp
• cp-algo/math/convolution.hpp
• cp-algo/math/cvector.hpp
• cp-algo/math/decimal.hpp
• cp-algo/math/fft.hpp
• cp-algo/math/fft64.hpp
• cp-algo/math/fft_simple.hpp
• cp-algo/math/laurent.hpp
• cp-algo/math/multivar.hpp
• cp-algo/math/poly.hpp
• cp-algo/math/poly/impl/div.hpp
• cp-algo/math/poly/impl/euclid.hpp
• cp-algo/number_theory.hpp
• cp-algo/number_theory/two_squares.hpp

Verified with

• Addition of Big Integers (verify/bigint/addition.test.cpp)
• Division of Big Integers (verify/bigint/divmod.test.cpp)
• Addition of Hex Big Integers (verify/bigint/hex_addition.test.cpp)
• Division of Hex Big Integers (verify/bigint/hex_divmod.test.cpp)
• Multiplication of Hex Big Integers (verify/bigint/hex_multiplication.test.cpp)
• Multiplication of Big Integers (verify/bigint/multiplication.test.cpp)
• Sort Points by Argument (verify/geometry/arg.test.cpp)
• Closest Pair of Points (verify/geometry/closest_pair.test.cpp)
• Static Convex Hull (verify/geometry/hull.test.cpp)
• Characteristic Polynomial (verify/linalg/characteristic.test.cpp)
• Pow of Matrix (Frobenius) (verify/linalg/pow_fast.test.cpp)
• Represent A Number As Two Square Sum (verify/number_theory/two_squares.test.cpp)
• Bell Number (verify/poly/bell.test.cpp)
• Convolution mod $10^9+7$ (verify/poly/convolution107.test.cpp)
• Convolution (Mod $2^{64}$) (verify/poly/convolution64.test.cpp)
• Convolution (Large) (verify/poly/convolution_large.test.cpp)
• Convolution on the Multiplicative Monoid of $\mathbb Z/p\mathbb{Z}$ (verify/poly/convolution_mul.test.cpp)
• Division of Polynomials (verify/poly/div.test.cpp)
• Exp of Power Series (verify/poly/exp.test.cpp)
• Find Linear Recurrence (verify/poly/find_linrec.test.cpp)
• Polynomial Interpolation (verify/poly/inter.test.cpp)
• Polynomial Interpolation (Geometric Sequence) (verify/poly/inter_geo.test.cpp)
• Inv of Power Series (verify/poly/inv.test.cpp)
• Inv of Polynomials (verify/poly/inv_mod.test.cpp)
• Consecutive Terms of Linear Recursion (verify/poly/linrec_consecutive.test.cpp)
• Kth term of Linearly Recurrent Sequence (verify/poly/linrec_single.test.cpp)
• Log of Power Series (verify/poly/log.test.cpp)
• Multipoint Evaluation (verify/poly/multieval.test.cpp)
• Multipoint Evaluation (Geometric Sequence) (verify/poly/multieval_geo.test.cpp)
• Multidimensional Convolution (Truncated) (verify/poly/multivar.test.cpp)
• Conversion to Newton Basis (verify/poly/newton.test.cpp)
• Pow of Power Series (verify/poly/pow.test.cpp)
• Product of Polynomial Sequence (verify/poly/prod_sequence.test.cpp)
• Polynomial Root Finding (verify/poly/roots.test.cpp)
• Shift of Sampling Points of Polynomial (verify/poly/sampling.test.cpp)
• Sqrt of Power Series (verify/poly/sqrt.test.cpp)
• $\sum\limits_{i=0}^\infty r^i i^d$ (verify/poly/sum_exp_poly.test.cpp)
• Polynomial Taylor Shift (verify/poly/taylor.test.cpp)
• Wildcard Pattern Matching (verify/poly/wildcard.test.cpp)

Code

#ifndef CP_ALGO_UTIL_COMPLEX_HPP
#define CP_ALGO_UTIL_COMPLEX_HPP
#include <iostream>
#include <cmath>
#include <type_traits>
#include "simd.hpp"
CP_ALGO_SIMD_PRAGMA_PUSH
namespace cp_algo {
    // Custom implementation, since std::complex is UB on non-floating types
    template<typename T>
    struct complex {
        using value_type = T;
        T x, y;
        inline constexpr complex(): x(), y() {}
        inline constexpr complex(T const& x): x(x), y() {}
        inline constexpr complex(T const& x, T const& y): x(x), y(y) {}
        inline complex& operator *= (T const& t) {x *= t; y *= t; return *this;}
        inline complex& operator /= (T const& t) {x /= t; y /= t; return *this;}
        inline complex operator * (T const& t) const {return complex(*this) *= t;}
        inline complex operator / (T const& t) const {return complex(*this) /= t;}
        inline complex& operator += (complex const& t) {x += t.x; y += t.y; return *this;}
        inline complex& operator -= (complex const& t) {x -= t.x; y -= t.y; return *this;}
        inline complex operator * (complex const& t) const {return {x * t.x - y * t.y, x * t.y + y * t.x};}
        inline complex operator / (complex const& t) const {return *this * t.conj() / t.norm();}
        inline complex operator + (complex const& t) const {return complex(*this) += t;}
        inline complex operator - (complex const& t) const {return complex(*this) -= t;}
        inline complex& operator *= (complex const& t) {return *this = *this * t;}
        inline complex& operator /= (complex const& t) {return *this = *this / t;}
        inline complex operator - () const {return {-x, -y};}
        inline complex conj() const {return {x, -y};}
        inline T norm() const {return x * x + y * y;}
        inline T abs() const {return std::sqrt(norm());}
        inline T const real() const {return x;}
        inline T const imag() const {return y;}
        inline T& real() {return x;}
        inline T& imag() {return y;}
        inline static constexpr complex polar(T r, T theta) {return {T(r * cos(theta)), T(r * sin(theta))};}
        inline auto operator <=> (complex const& t) const = default;
    };
    template<typename T> inline complex<T> conj(complex<T> const& x) {return x.conj();}
    template<typename T> inline T norm(complex<T> const& x) {return x.norm();}
    template<typename T> inline T abs(complex<T> const& x) {return x.abs();}
    template<typename T> inline T& real(complex<T> &x) {return x.real();}
    template<typename T> inline T& imag(complex<T> &x) {return x.imag();}
    template<typename T> inline T const real(complex<T> const& x) {return x.real();}
    template<typename T> inline T const imag(complex<T> const& x) {return x.imag();}
    template<typename T>
    inline constexpr complex<T> polar(T r, T theta) {
        return complex<T>::polar(r, theta);
    }
    template<typename T>
    inline std::ostream& operator << (std::ostream &out, complex<T> const& x) {
        return out << x.real() << ' ' << x.imag();
    }
}
#pragma GCC pop_options
#endif // CP_ALGO_UTIL_COMPLEX_HPP

#line 1 "cp-algo/util/complex.hpp"


#include <iostream>
#include <cmath>
#include <type_traits>
#line 1 "cp-algo/util/simd.hpp"


#include <experimental/simd>
#include <cstdint>
#include <cstddef>
#include <memory>

#if defined(__x86_64__) && !defined(CP_ALGO_DISABLE_AVX2)
#define CP_ALGO_SIMD_AVX2_TARGET _Pragma("GCC target(\"avx2\")")
#else
#define CP_ALGO_SIMD_AVX2_TARGET
#endif

#define CP_ALGO_SIMD_PRAGMA_PUSH \
    _Pragma("GCC push_options") \
    CP_ALGO_SIMD_AVX2_TARGET

CP_ALGO_SIMD_PRAGMA_PUSH
namespace cp_algo {
    template<typename T, size_t len>
    using simd [[gnu::vector_size(len * sizeof(T))]] = T;
    using u64x8 = simd<uint64_t, 8>;
    using u32x16 = simd<uint32_t, 16>;
    using i64x4 = simd<int64_t, 4>;
    using u64x4 = simd<uint64_t, 4>;
    using u32x8 = simd<uint32_t, 8>;
    using u16x16 = simd<uint16_t, 16>;
    using i32x4 = simd<int32_t, 4>;
    using u32x4 = simd<uint32_t, 4>;
    using u16x8 = simd<uint16_t, 8>;
    using u16x4 = simd<uint16_t, 4>;
    using i16x4 = simd<int16_t, 4>;
    using u8x32 = simd<uint8_t, 32>;
    using u8x8 = simd<uint8_t, 8>;
    using u8x4 = simd<uint8_t, 4>;
    using dx4 = simd<double, 4>;

    inline dx4 abs(dx4 a) {
        return dx4{
            std::abs(a[0]),
            std::abs(a[1]),
            std::abs(a[2]),
            std::abs(a[3])
        };
    }

    // https://stackoverflow.com/a/77376595
    // works for ints in (-2^51, 2^51)
    static constexpr dx4 magic = dx4() + (3ULL << 51);
    inline i64x4 lround(dx4 x) {
        return i64x4(x + magic) - i64x4(magic);
    }
    inline dx4 to_double(i64x4 x) {
        return dx4(x + i64x4(magic)) - magic;
    }

    inline dx4 round(dx4 a) {
        return dx4{
            std::nearbyint(a[0]),
            std::nearbyint(a[1]),
            std::nearbyint(a[2]),
            std::nearbyint(a[3])
        };
    }

    inline u64x4 low32(u64x4 x) {
        return x & uint32_t(-1);
    }
    inline auto swap_bytes(auto x) {
        return decltype(x)(__builtin_shufflevector(u32x8(x), u32x8(x), 1, 0, 3, 2, 5, 4, 7, 6));
    }
    inline u64x4 montgomery_reduce(u64x4 x, uint32_t mod, uint32_t imod) {
#ifdef __AVX2__
        auto x_ninv = u64x4(_mm256_mul_epu32(__m256i(x), __m256i() + imod));
        x += u64x4(_mm256_mul_epu32(__m256i(x_ninv), __m256i() + mod));
#else
        auto x_ninv = u64x4(u32x8(low32(x)) * imod);
        x += x_ninv * uint64_t(mod);
#endif
        return swap_bytes(x);
    }

    inline u64x4 montgomery_mul(u64x4 x, u64x4 y, uint32_t mod, uint32_t imod) {
#ifdef __AVX2__
        return montgomery_reduce(u64x4(_mm256_mul_epu32(__m256i(x), __m256i(y))), mod, imod);
#else
        return montgomery_reduce(x * y, mod, imod);
#endif
    }
    inline u32x8 montgomery_mul(u32x8 x, u32x8 y, uint32_t mod, uint32_t imod) {
        return u32x8(montgomery_mul(u64x4(x), u64x4(y), mod, imod)) |
               u32x8(swap_bytes(montgomery_mul(u64x4(swap_bytes(x)), u64x4(swap_bytes(y)), mod, imod)));
    }
    inline dx4 rotate_right(dx4 x) {
        static constexpr u64x4 shuffler = {3, 0, 1, 2};
        return __builtin_shuffle(x, shuffler);
    }

    template<std::size_t Align = 32>
    inline bool is_aligned(const auto* p) noexcept {
        return (reinterpret_cast<std::uintptr_t>(p) % Align) == 0;
    }

    template<class Target>
    inline Target& vector_cast(auto &&p) {
        return *reinterpret_cast<Target*>(std::assume_aligned<alignof(Target)>(&p));
    }
}
#pragma GCC pop_options

#line 7 "cp-algo/util/complex.hpp"
CP_ALGO_SIMD_PRAGMA_PUSH
namespace cp_algo {
    // Custom implementation, since std::complex is UB on non-floating types
    template<typename T>
    struct complex {
        using value_type = T;
        T x, y;
        inline constexpr complex(): x(), y() {}
        inline constexpr complex(T const& x): x(x), y() {}
        inline constexpr complex(T const& x, T const& y): x(x), y(y) {}
        inline complex& operator *= (T const& t) {x *= t; y *= t; return *this;}
        inline complex& operator /= (T const& t) {x /= t; y /= t; return *this;}
        inline complex operator * (T const& t) const {return complex(*this) *= t;}
        inline complex operator / (T const& t) const {return complex(*this) /= t;}
        inline complex& operator += (complex const& t) {x += t.x; y += t.y; return *this;}
        inline complex& operator -= (complex const& t) {x -= t.x; y -= t.y; return *this;}
        inline complex operator * (complex const& t) const {return {x * t.x - y * t.y, x * t.y + y * t.x};}
        inline complex operator / (complex const& t) const {return *this * t.conj() / t.norm();}
        inline complex operator + (complex const& t) const {return complex(*this) += t;}
        inline complex operator - (complex const& t) const {return complex(*this) -= t;}
        inline complex& operator *= (complex const& t) {return *this = *this * t;}
        inline complex& operator /= (complex const& t) {return *this = *this / t;}
        inline complex operator - () const {return {-x, -y};}
        inline complex conj() const {return {x, -y};}
        inline T norm() const {return x * x + y * y;}
        inline T abs() const {return std::sqrt(norm());}
        inline T const real() const {return x;}
        inline T const imag() const {return y;}
        inline T& real() {return x;}
        inline T& imag() {return y;}
        inline static constexpr complex polar(T r, T theta) {return {T(r * cos(theta)), T(r * sin(theta))};}
        inline auto operator <=> (complex const& t) const = default;
    };
    template<typename T> inline complex<T> conj(complex<T> const& x) {return x.conj();}
    template<typename T> inline T norm(complex<T> const& x) {return x.norm();}
    template<typename T> inline T abs(complex<T> const& x) {return x.abs();}
    template<typename T> inline T& real(complex<T> &x) {return x.real();}
    template<typename T> inline T& imag(complex<T> &x) {return x.imag();}
    template<typename T> inline T const real(complex<T> const& x) {return x.real();}
    template<typename T> inline T const imag(complex<T> const& x) {return x.imag();}
    template<typename T>
    inline constexpr complex<T> polar(T r, T theta) {
        return complex<T>::polar(r, theta);
    }
    template<typename T>
    inline std::ostream& operator << (std::ostream &out, complex<T> const& x) {
        return out << x.real() << ' ' << x.imag();
    }
}
#pragma GCC pop_options

#ifndef CP_ALGO_UTIL_COMPLEX_HPP
#define CP_ALGO_UTIL_COMPLEX_HPP
#include <iostream>
#include <cmath>
#include <type_traits>
#include "simd.hpp"
CP_ALGO_SIMD_PRAGMA_PUSHnamespace cp_algo{template<typename T>struct complex{using value_type=T;T x,y;inline constexpr complex():x(),y(){}inline constexpr complex(T const&x):x(x),y(){}inline constexpr complex(T const&x,T const&y):x(x),y(y){}inline complex&operator*=(T const&t){x*=t;y*=t;return*this;}inline complex&operator/=(T const&t){x/=t;y/=t;return*this;}inline complex operator*(T const&t)const{return complex(*this)*=t;}inline complex operator/(T const&t)const{return complex(*this)/=t;}inline complex&operator+=(complex const&t){x+=t.x;y+=t.y;return*this;}inline complex&operator-=(complex const&t){x-=t.x;y-=t.y;return*this;}inline complex operator*(complex const&t)const{return{x*t.x-y*t.y,x*t.y+y*t.x};}inline complex operator/(complex const&t)const{return*this*t.conj()/t.norm();}inline complex operator+(complex const&t)const{return complex(*this)+=t;}inline complex operator-(complex const&t)const{return complex(*this)-=t;}inline complex&operator*=(complex const&t){return*this=*this*t;}inline complex&operator/=(complex const&t){return*this=*this/t;}inline complex operator-()const{return{-x,-y};}inline complex conj()const{return{x,-y};}inline T norm()const{return x*x+y*y;}inline T abs()const{return std::sqrt(norm());}inline T const real()const{return x;}inline T const imag()const{return y;}inline T&real(){return x;}inline T&imag(){return y;}inline static constexpr complex polar(T r,T theta){return{T(r*cos(theta)),T(r*sin(theta))};}inline auto operator<=>(complex const&t)const=default;};template<typename T>inline complex<T>conj(complex<T>const&x){return x.conj();}template<typename T>inline T norm(complex<T>const&x){return x.norm();}template<typename T>inline T abs(complex<T>const&x){return x.abs();}template<typename T>inline T&real(complex<T>&x){return x.real();}template<typename T>inline T&imag(complex<T>&x){return x.imag();}template<typename T>inline T const real(complex<T>const&x){return x.real();}template<typename T>inline T const imag(complex<T>const&x){return x.imag();}template<typename T>inline constexpr complex<T>polar(T r,T theta){return complex<T>::polar(r,theta);}template<typename T>inline std::ostream&operator<<(std::ostream&out,complex<T>const&x){return out<<x.real()<<' '<<x.imag();}}
#pragma GCC pop_options
#endif

#line 1 "cp-algo/util/complex.hpp"
#include <iostream>
#include <cmath>
#include <type_traits>
#line 1 "cp-algo/util/simd.hpp"
#include <experimental/simd>
#include <cstdint>
#include <cstddef>
#include <memory>
#if defined(__x86_64__) && !defined(CP_ALGO_DISABLE_AVX2)
#define CP_ALGO_SIMD_AVX2_TARGET _Pragma("GCC target(\"avx2\")")
#else
#define CP_ALGO_SIMD_AVX2_TARGET
#endif
#define CP_ALGO_SIMD_PRAGMA_PUSH \
_Pragma("GCC push_options")\CP_ALGO_SIMD_AVX2_TARGETCP_ALGO_SIMD_PRAGMA_PUSHnamespace cp_algo{template<typename T,size_t len>using simd[[gnu::vector_size(len*sizeof(T))]]=T;using u64x8=simd<uint64_t,8>;using u32x16=simd<uint32_t,16>;using i64x4=simd<int64_t,4>;using u64x4=simd<uint64_t,4>;using u32x8=simd<uint32_t,8>;using u16x16=simd<uint16_t,16>;using i32x4=simd<int32_t,4>;using u32x4=simd<uint32_t,4>;using u16x8=simd<uint16_t,8>;using u16x4=simd<uint16_t,4>;using i16x4=simd<int16_t,4>;using u8x32=simd<uint8_t,32>;using u8x8=simd<uint8_t,8>;using u8x4=simd<uint8_t,4>;using dx4=simd<double,4>;inline dx4 abs(dx4 a){return dx4{std::abs(a[0]),std::abs(a[1]),std::abs(a[2]),std::abs(a[3])};}static constexpr dx4 magic=dx4()+(3ULL<<51);inline i64x4 lround(dx4 x){return i64x4(x+magic)-i64x4(magic);}inline dx4 to_double(i64x4 x){return dx4(x+i64x4(magic))-magic;}inline dx4 round(dx4 a){return dx4{std::nearbyint(a[0]),std::nearbyint(a[1]),std::nearbyint(a[2]),std::nearbyint(a[3])};}inline u64x4 low32(u64x4 x){return x&uint32_t(-1);}inline auto swap_bytes(auto x){return decltype(x)(__builtin_shufflevector(u32x8(x),u32x8(x),1,0,3,2,5,4,7,6));}inline u64x4 montgomery_reduce(u64x4 x,uint32_t mod,uint32_t imod){
#ifdef __AVX2__
auto x_ninv=u64x4(_mm256_mul_epu32(__m256i(x),__m256i()+imod));x+=u64x4(_mm256_mul_epu32(__m256i(x_ninv),__m256i()+mod));
#else
auto x_ninv=u64x4(u32x8(low32(x))*imod);x+=x_ninv*uint64_t(mod);
#endif
return swap_bytes(x);}inline u64x4 montgomery_mul(u64x4 x,u64x4 y,uint32_t mod,uint32_t imod){
#ifdef __AVX2__
return montgomery_reduce(u64x4(_mm256_mul_epu32(__m256i(x),__m256i(y))),mod,imod);
#else
return montgomery_reduce(x*y,mod,imod);
#endif
}inline u32x8 montgomery_mul(u32x8 x,u32x8 y,uint32_t mod,uint32_t imod){return u32x8(montgomery_mul(u64x4(x),u64x4(y),mod,imod))|u32x8(swap_bytes(montgomery_mul(u64x4(swap_bytes(x)),u64x4(swap_bytes(y)),mod,imod)));}inline dx4 rotate_right(dx4 x){static constexpr u64x4 shuffler={3,0,1,2};return __builtin_shuffle(x,shuffler);}template<std::size_t Align=32>inline bool is_aligned(const auto*p)noexcept{return(reinterpret_cast<std::uintptr_t>(p)%Align)==0;}template<class Target>inline Target&vector_cast(auto&&p){return*reinterpret_cast<Target*>(std::assume_aligned<alignof(Target)>(&p));}}
#pragma GCC pop_options
#line 7 "cp-algo/util/complex.hpp"
CP_ALGO_SIMD_PRAGMA_PUSHnamespace cp_algo{template<typename T>struct complex{using value_type=T;T x,y;inline constexpr complex():x(),y(){}inline constexpr complex(T const&x):x(x),y(){}inline constexpr complex(T const&x,T const&y):x(x),y(y){}inline complex&operator*=(T const&t){x*=t;y*=t;return*this;}inline complex&operator/=(T const&t){x/=t;y/=t;return*this;}inline complex operator*(T const&t)const{return complex(*this)*=t;}inline complex operator/(T const&t)const{return complex(*this)/=t;}inline complex&operator+=(complex const&t){x+=t.x;y+=t.y;return*this;}inline complex&operator-=(complex const&t){x-=t.x;y-=t.y;return*this;}inline complex operator*(complex const&t)const{return{x*t.x-y*t.y,x*t.y+y*t.x};}inline complex operator/(complex const&t)const{return*this*t.conj()/t.norm();}inline complex operator+(complex const&t)const{return complex(*this)+=t;}inline complex operator-(complex const&t)const{return complex(*this)-=t;}inline complex&operator*=(complex const&t){return*this=*this*t;}inline complex&operator/=(complex const&t){return*this=*this/t;}inline complex operator-()const{return{-x,-y};}inline complex conj()const{return{x,-y};}inline T norm()const{return x*x+y*y;}inline T abs()const{return std::sqrt(norm());}inline T const real()const{return x;}inline T const imag()const{return y;}inline T&real(){return x;}inline T&imag(){return y;}inline static constexpr complex polar(T r,T theta){return{T(r*cos(theta)),T(r*sin(theta))};}inline auto operator<=>(complex const&t)const=default;};template<typename T>inline complex<T>conj(complex<T>const&x){return x.conj();}template<typename T>inline T norm(complex<T>const&x){return x.norm();}template<typename T>inline T abs(complex<T>const&x){return x.abs();}template<typename T>inline T&real(complex<T>&x){return x.real();}template<typename T>inline T&imag(complex<T>&x){return x.imag();}template<typename T>inline T const real(complex<T>const&x){return x.real();}template<typename T>inline T const imag(complex<T>const&x){return x.imag();}template<typename T>inline constexpr complex<T>polar(T r,T theta){return complex<T>::polar(r,theta);}template<typename T>inline std::ostream&operator<<(std::ostream&out,complex<T>const&x){return out<<x.real()<<' '<<x.imag();}}
#pragma GCC pop_options

Back to top page