CP-Algorithms Library

This documentation is automatically generated by competitive-verifier/competitive-verifier

View the Project on GitHub cp-algorithms/cp-algorithms-aux

:heavy_check_mark: Represent A Number As Two Square Sum (verify/number_theory/two_squares.test.cpp)

Depends on

Code

// @brief Represent A Number As Two Square Sum
#define PROBLEM "https://judge.yosupo.jp/problem/two_square_sum"
#pragma GCC optimize("O3,unroll-loops")
#pragma GCC target("avx2")
#include "cp-algo/number_theory/two_squares.hpp"
#include <bits/stdc++.h>

using namespace std;
using namespace cp_algo::math;

void solve() {
	int64_t n;
	cin >> n;
	auto res = two_squares_all(n);
	cout << size(res) << "\n";
	for(auto p: res) {
		cout << p.real() << ' ' << p.imag() << "\n";
	}
}

signed main() {
    //freopen("input.txt", "r", stdin);
    ios::sync_with_stdio(0);
    cin.tie(0);
    int t = 1;
    cin >> t;
	while(t--) {
		solve();
	}
}

#line 1 "verify/number_theory/two_squares.test.cpp"
// @brief Represent A Number As Two Square Sum
#define PROBLEM "https://judge.yosupo.jp/problem/two_square_sum"
#pragma GCC optimize("O3,unroll-loops")
#pragma GCC target("avx2")
#line 1 "cp-algo/number_theory/two_squares.hpp"


#line 1 "cp-algo/number_theory/euler.hpp"


#line 1 "cp-algo/number_theory/factorize.hpp"


#line 1 "cp-algo/number_theory/primality.hpp"


#line 1 "cp-algo/number_theory/modint.hpp"


#line 1 "cp-algo/math/common.hpp"


#include <functional>
#include <cstdint>
#include <cassert>
#include <bit>
#include <vector>
#include <algorithm>
namespace cp_algo::math {
#ifdef CP_ALGO_MAXN
    const int maxn = CP_ALGO_MAXN;
#else
    const int maxn = 1 << 19;
#endif
    const int magic = 64; // threshold for sizes to run the naive algo

    // Nonnegative 64-bit exponents, with an associative operation and its identity.
    // Windows >1 precompute odd powers only when that saves operations.
    template<int window = 1>
    auto bpow(auto const& x, auto n, auto const& one, auto op) {
        static_assert(window >= 1 && window <= 6);
        if constexpr(window > 1) {
            if(n == 0) {return one;}
            int bits = std::bit_width(uint64_t(n));
            auto low_bit = [&](int high) {
                int low = std::max(0, high - window + 1);
                while(!((n >> low) & 1)) {low++;}
                return low;
            };
            int first = low_bit(bits - 1);
            int cost = (1 << (window - 1)) + first;
            for(int j = first - 1; j >= 0;) {
                if(!((n >> j) & 1)) {j--;}
                else {cost++; j = low_bit(j) - 1;}
            }
            // Do not pay for the table when binary powering uses fewer operations.
            if(cost >= bits + std::popcount(uint64_t(n)) - 2) {return bpow<1>(x, n, one, op);}
            using T = std::decay_t<decltype(x)>;
            std::vector<T> odd;
            odd.reserve(1 << (window - 1));
            odd.push_back(x);
            auto square = op(x, x);
            while(odd.size() < size_t(1 << (window - 1))) {odd.push_back(op(odd.back(), square));}
            auto ans = odd[(n >> first) / 2];
            for(int j = first - 1; j >= 0;) {
                if(!((n >> j) & 1)) {ans = op(ans, ans); j--;}
                else {
                    int low = low_bit(j), length = j - low + 1;
                    auto digit = (n >> low) & ((1u << length) - 1);
                    for(int i = 0; i < length; i++) {ans = op(ans, ans);}
                    ans = op(ans, odd[digit / 2]);
                    j = low - 1;
                }
            }
            return ans;
        } else {
            if (n == 0) {
                return one;
            }
            auto ans = x;
            for(int j = std::bit_width<uint64_t>(n) - 2; ~j; j--) {
                ans = op(ans, ans);
                if((n >> j) & 1) {
                    ans = op(ans, x);
                }
            }
            return ans;
        }
    }
    template<int window = 1>
    auto bpow(auto x, auto n, auto ans) {
        return bpow<window>(x, n, ans, std::multiplies{});
    }
    template<typename T>
    T bpow(T const& x, auto n) {
        return bpow(x, n, T(1));
    }
    inline constexpr auto inv2(auto x) {
        assert(x % 2);
        std::make_unsigned_t<decltype(x)> y = 1;
        while(y * x != 1) {
            y *= 2 - x * y;
        }
        return y;
    }
}

#line 4 "cp-algo/number_theory/modint.hpp"
#include <iostream>
#line 6 "cp-algo/number_theory/modint.hpp"
namespace cp_algo::math {

    template<typename modint, typename _Int>
    struct modint_base {
        using Int = _Int;
        using UInt = std::make_unsigned_t<Int>;
        static constexpr size_t bits = sizeof(Int) * 8;
        using Int2 = std::conditional_t<bits <= 32, int64_t, __int128_t>;
        using UInt2 = std::conditional_t<bits <= 32, uint64_t, __uint128_t>;
        constexpr static Int mod() {
            return modint::mod();
        }
        constexpr static UInt remod() {
            return modint::remod();
        }
        constexpr static UInt2 modmod() {
            return UInt2(mod()) * mod();
        }
        constexpr modint_base() = default;
        constexpr modint_base(Int2 rr) {
            to_modint().setr(UInt((rr + modmod()) % mod()));
        }
        constexpr modint inv() const {
            return bpow(to_modint(), mod() - 2);
        }
        modint operator-() const {
            modint neg;
            neg.r = std::min(-r, remod() - r);
            return neg;
        }
        modint& operator/=(const modint& t) {
            return to_modint() *= t.inv();
        }
        modint& operator*=(const modint& t) {
            r = UInt(UInt2(r) * t.r % mod());
            return to_modint();
        }
        modint& operator+=(const modint& t) {
            r += t.r;
            r = std::min(r, r - remod());
            return to_modint();
        }
        modint& operator-=(const modint& t) {
            r -= t.r;
            r = std::min(r, r + remod());
            return to_modint();
        }
        modint operator+(const modint& t) const { return modint(to_modint()) += t; }
        modint operator-(const modint& t) const { return modint(to_modint()) -= t; }
        modint operator*(const modint& t) const { return modint(to_modint()) *= t; }
        modint operator/(const modint& t) const { return modint(to_modint()) /= t; }
        // Why <=> doesn't work?..
        auto operator==(const modint& t) const { return to_modint().getr() == t.getr(); }
        auto operator!=(const modint& t) const { return to_modint().getr() != t.getr(); }
        auto operator<=(const modint& t) const { return to_modint().getr() <= t.getr(); }
        auto operator>=(const modint& t) const { return to_modint().getr() >= t.getr(); }
        auto operator<(const modint& t) const { return to_modint().getr() < t.getr(); }
        auto operator>(const modint& t) const { return to_modint().getr() > t.getr(); }
        Int rem() const {
            UInt R = to_modint().getr();
            return R - (R > (UInt)mod() / 2) * mod();
        }
        constexpr void setr(UInt rr) {
            r = rr;
        }
        constexpr UInt getr() const {
            return r;
        }

        // Only use these if you really know what you're doing!
        static uint64_t modmod8() { return uint64_t(8 * modmod()); }
        void add_unsafe(UInt t) { r += t; }
        void pseudonormalize() { r = std::min(r, r - modmod8()); }
        modint const& normalize() {
            if(r >= (UInt)mod()) {
                r %= mod();
            }
            return to_modint();
        }
        void setr_direct(UInt rr) { r = rr; }
        UInt getr_direct() const { return r; }

    protected:
        UInt r;

    private:
        constexpr modint& to_modint() { return static_cast<modint&>(*this); }
        constexpr modint const& to_modint() const { return static_cast<modint const&>(*this); }
    };
    template<typename modint>
    concept modint_type = std::is_base_of_v<modint_base<modint, typename modint::Int>, modint>;
    template<modint_type modint>
    decltype(std::cin)& operator>>(decltype(std::cin)& in, modint& x) {
        typename modint::UInt r;
        auto &res = in >> r;
        x.setr(r);
        return res;
    }
    template<modint_type modint>
    decltype(std::cout)& operator<<(decltype(std::cout)& out, modint const& x) {
        return out << x.getr();
    }

    template<auto m>
    struct modint: modint_base<modint<m>, decltype(m)> {
        using Base = modint_base<modint<m>, decltype(m)>;
        using Base::Base;
        static constexpr Base::Int mod() { return m; }
        static constexpr Base::UInt remod() { return m; }
        auto getr() const { return Base::r; }
    };

    // Odd moduli up to a quarter of the unsigned word keep Montgomery residues lazily in
    // [0, 2 mod): remod() = 2 mod, and both 4 mod and ab + q * mod fit. Any other modulus keeps
    // fully reduced residues, remod() = mod, which is what the sums of modint_base need then:
    // an even one without the Montgomery form, a wide odd one with a reduction that subtracts
    // high words instead of adding double words that would overflow.
    template<typename Int = int>
    struct dynamic_modint: modint_base<dynamic_modint<Int>, Int> {
        using Base = modint_base<dynamic_modint<Int>, Int>;
        using Base::Base;

        // Out of line, so that the hot path stays as small as it was.
        [[gnu::noinline, gnu::cold]] static Base::UInt m_reduce_reduced(Base::UInt2 ab) {
            if(mod() % 2 == 0) {
                return typename Base::UInt(ab % mod());
            }
            // q * mod has the low word of ab, so the difference of the high words is exact.
            typename Base::UInt q = -(typename Base::UInt(ab) * inverse);
            auto high = typename Base::UInt(ab >> Base::bits);
            auto low = typename Base::UInt(typename Base::UInt2(q) * typename Base::UInt(mod()) >> Base::bits);
            return high >= low ? high - low : high - low + mod();
        }
        static Base::UInt m_reduce(Base::UInt2 ab) {
            if(imod() == 0) [[unlikely]] {
                return m_reduce_reduced(ab);
            }
            typename Base::UInt2 m = typename Base::UInt(ab) * imod();
            return typename Base::UInt((ab + m * mod()) >> Base::bits);
        }
        static Base::UInt m_transform(Base::UInt a) {
            if(mod() % 2 == 0) [[unlikely]] {
                return a;
            } else {
                return m_reduce(a * pw128());
            }
        }
        dynamic_modint& operator*=(const dynamic_modint& t) {
            Base::r = m_reduce(typename Base::UInt2(Base::r) * t.r);
            return *this;
        }
        void setr(Base::UInt rr) {
            Base::r = m_transform(rr);
        }
        Base::UInt getr() const {
            typename Base::UInt res = m_reduce(Base::r);
            return std::min(res, res - mod());
        }
        static Int mod() { return m; }
        static Base::UInt remod() { return rm; }
        static Base::UInt imod() { return im; }
        static Base::UInt2 pw128() { return r2; }
        static void switch_mod(Int nm) {
            m = nm;
            bool lazy = m % 2 && typename Base::UInt(m) <= typename Base::UInt(-1) / 4;
            rm = typename Base::UInt(m) * (lazy ? 2 : 1);
            inverse = m % 2 ? inv2(-m) : 0;
            im = lazy ? inverse : 0;
            r2 = static_cast<Base::UInt>(static_cast<Base::UInt2>(-1) % m + 1);
        }

        // Wrapper for temp switching
        auto static with_mod(Int tmp, auto callback) {
            struct scoped {
                Int prev = mod();
                ~scoped() { switch_mod(prev); }
            } _;
            switch_mod(tmp);
            return callback();
        }

    private:
        static inline thread_local Int m = 1;
        // im: -1 / mod modulo 2^bits for lazy residues and 0 for reduced ones; inverse: the same
        // for every odd mod; rm: the value of remod().
        static inline thread_local Base::UInt im = -1, r2 = 0, inverse = -1, rm = 2;
    };
}

#line 6 "cp-algo/number_theory/primality.hpp"
namespace cp_algo::math {
    // https://en.wikipedia.org/wiki/Miller–Rabin_primality_test
    template<typename _Int>
    bool is_prime(_Int m) {
        using Int = std::make_signed_t<_Int>;
        using UInt = std::make_unsigned_t<Int>;
        if(m == 1 || m % 2 == 0) {
            return m == 2;
        }
        // m - 1 = 2^s * d
        int s = std::countr_zero(UInt(m - 1));
        auto d = (m - 1) >> s;
        using base = dynamic_modint<Int>;
        auto test = [&](base x) {
            x = bpow(x, d);
            if(std::abs(x.rem()) <= 1) {
                return true;
            }
            for(int i = 1; i < s && x != -1; i++) {
                x *= x;
            }
            return x == -1;
        };
        return base::with_mod(m, [&]() {
#ifdef CP_ALGO_NUMBER_THEORY_PRIMALITY_BASES_HPP
            uint16_t base2 = 7, base3 = 61;
            if (m != uint32_t(m)) {
                base2 = base_table1[uint32_t(m * 0xAD625B89) >> 18];
                base3 = base_table2[base2 >> 13];
            }
            return test(2) && test(base2) && test(base3);
#else
            return std::ranges::all_of(std::array{2, 325, 9375, 28178, 450775, 9780504, 1795265022}, test);
#endif
        });
    }
}

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



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



#include <set>
#include <map>
#include <deque>
#include <stack>
#include <queue>
#line 10 "cp-algo/util/big_alloc.hpp"
#include <string>
#include <cstddef>
#line 14 "cp-algo/util/big_alloc.hpp"
#include <forward_list>

// Single macro to detect POSIX platforms (Linux, Unix, macOS)
#if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
#define CP_ALGO_USE_MMAP 1
#include <sys/mman.h>
#else
#define CP_ALGO_USE_MMAP 0
#endif

namespace cp_algo {
    template<typename T, size_t Align = 32>
    class big_alloc {
        static_assert(Align >= alignof(void*), "Align must be at least pointer-size");
        static_assert(std::has_single_bit(Align), "Align must be a power of two");

    public:
        using value_type = T;
        template<class U>
        struct rebind {
            using other = big_alloc<U, Align>;
        };
        constexpr bool operator==(const big_alloc&) const = default;

        big_alloc() noexcept = default;
        template<typename U, std::size_t A>
        big_alloc(const big_alloc<U, A>&) noexcept {}

        [[nodiscard]] T* allocate(std::size_t n) {
            std::size_t padded = round_up(n * sizeof(T));
            std::size_t align = std::max<std::size_t>(alignof(T), Align);
#if CP_ALGO_USE_MMAP
            if(padded >= MEGABYTE) {
                void* raw = mmap(nullptr, padded, PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
                madvise(raw, padded, MADV_HUGEPAGE);
                return static_cast<T*>(raw);
            }
#endif
            return static_cast<T*>(::operator new(padded, std::align_val_t(align)));
        }

        void deallocate(T* p, std::size_t n) noexcept {
            if(!p) {
                return;
            }
            std::size_t padded = round_up(n * sizeof(T));
            std::size_t align = std::max<std::size_t>(alignof(T), Align);
#if CP_ALGO_USE_MMAP
            if(padded >= MEGABYTE) {
                munmap(p, padded);
                return;
            }
#endif
            ::operator delete(p, padded, std::align_val_t(align));
        }

    private:
        static constexpr std::size_t MEGABYTE = 1 << 20;
        static constexpr std::size_t round_up(std::size_t x) noexcept {
            return (x + Align - 1) / Align * Align;
        }
    };

    template<typename T>
    using big_vector = std::vector<T, big_alloc<T>>;
    template<typename T>
    using big_basic_string = std::basic_string<T, std::char_traits<T>, big_alloc<T>>;
    template<typename T>
    using big_deque = std::deque<T, big_alloc<T>>;
    template<typename T>
    using big_stack = std::stack<T, big_deque<T>>;
    template<typename T>
    using big_queue = std::queue<T, big_deque<T>>;
    template<typename T>
    using big_priority_queue = std::priority_queue<T, big_vector<T>>;
    template<typename T>
    using big_forward_list = std::forward_list<T, big_alloc<T>>;
    using big_string = big_basic_string<char>;

    template<typename Key, typename Value, typename Compare = std::less<Key>>
    using big_map = std::map<Key, Value, Compare, big_alloc<std::pair<const Key, Value>>>;
    template<typename T, typename Compare = std::less<T>>
    using big_multiset = std::multiset<T, Compare, big_alloc<T>>;
    template<typename T, typename Compare = std::less<T>>
    using big_set = std::set<T, Compare, big_alloc<T>>;
}


#line 5 "cp-algo/util/generator.hpp"
#include <generator>

namespace cp_algo {
    template<typename Ref, typename V = void>
    using big_generator = std::generator<Ref, V, big_alloc<std::byte>>;
}

namespace std::ranges {
    template<typename Ref, typename V>
    elements_of(cp_algo::big_generator<Ref, V>&&) -> elements_of<cp_algo::big_generator<Ref, V>&&, cp_algo::big_alloc<std::byte>>;
}


#line 1 "cp-algo/random/rng.hpp"


#include <chrono>
#include <random>
namespace cp_algo::random {
    std::mt19937_64 gen(
        std::chrono::steady_clock::now().time_since_epoch().count()
    );
    uint64_t rng() {
        return gen();
    }
}

#line 6 "cp-algo/number_theory/factorize.hpp"
namespace cp_algo::math {
    // https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
    template<typename _Int>
    auto proper_divisor(_Int m) {
        using Int = std::make_signed_t<_Int>;
        using base = dynamic_modint<Int>;
        return m % 2 == 0 ? 2 : base::with_mod(m, [&]() {
            base t = random::rng();
            auto f = [&](auto x) {
                return x * x + t;
            };
            const base zero = 0, one = 1;
            base x = zero, y = zero;
            base g = one;
            while(g == one) {
                for(int i = 0; i < 64; i++) {
                    x = f(x);
                    y = f(f(y));
                    if(x == y) [[unlikely]] {
                        t = random::rng();
                        x = y = zero;
                    } else {
                        base t = g * (x - y);
                        g = t == zero ? g : t;
                    }
                }
                g = std::gcd(g.getr(), m);
            }
            return g.getr();
        });
    }
    template<typename Int>
    big_generator<Int> factorize(Int m) {
        if(is_prime(m)) {
            co_yield m;
        } else if(m > 1) {
            auto g = proper_divisor(m);
            co_yield std::ranges::elements_of(factorize(g));
            co_yield std::ranges::elements_of(factorize(m / g));
        }
    }
    template<typename Int>
    big_generator<Int> divisors_sqrt(Int m) {
        for(Int i = 1; i * i <= m; i++) {
            if(m % i == 0) {
                co_yield i;
                if(i * i != m) {
                    co_yield m / i;
                }
            }
        }
    }
}

#line 5 "cp-algo/number_theory/euler.hpp"
namespace cp_algo::math {
    auto euler_phi(auto m) {
        using T = std::decay_t<decltype(m)>;
        auto primes = to<big_vector<T>>(factorize(m));
        std::ranges::sort(primes);
        auto [from, to] = std::ranges::unique(primes);
        primes.erase(from, to);
        auto ans = m;
        for(auto it: primes) {
            ans -= ans / it;
        }
        return ans;
    }
    template<modint_type base>
    auto period(base x) {
        auto ans = euler_phi(base::mod());
        base x0 = bpow(x, ans);
        for(auto t: factorize(ans)) {
            while(ans % t == 0 && x0 * bpow(x, ans / t) == x0) {
                ans /= t;
            }
        }
        return ans;
    }
    template<typename _Int>
    _Int primitive_root(_Int p) {
        using Int = std::make_signed_t<_Int>;
        using base = dynamic_modint<Int>;
        return base::with_mod(p, [p](){
            base t = 1;
            while(period(t) != p - 1) {
                t = random::rng();
            }
            return t.getr();
        });
    }
}

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


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


#include <experimental/simd>
#line 6 "cp-algo/util/simd.hpp"
#include <memory>

#if defined(__x86_64__) && !defined(CP_ALGO_DISABLE_AVX2)
#define CP_ALGO_SIMD_AVX2_TARGET _Pragma("GCC target(\"avx2,fma\")")
#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 u8x16 = simd<uint8_t, 16>;
    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);
    }

    // x - mod where x >= mod, for x in [0, 2 mod). Unsigned, so it holds up to mod < 2^31;
    // zero upper halves of 64-bit lanes stay zero.
    inline u32x8 reduce_once(u32x8 x, uint32_t mod) {
        auto y = x - mod;
        return x < y ? x : y;
    }
    inline u64x4 reduce_once(u64x4 x, uint32_t mod) {
        return u64x4(reduce_once(u32x8(x), mod));
    }

    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

#line 7 "cp-algo/number_theory/two_squares.hpp"
#include <utility>
#line 10 "cp-algo/number_theory/two_squares.hpp"
namespace cp_algo::math {
    template<typename T>
    using gaussint = complex<T>;
    template<typename _Int>
    auto two_squares_prime_any(_Int p) {
        if(p == 2) {
            return gaussint<_Int>{1, 1};
        }
        assert(p % 4 == 1);
        using Int = std::make_signed_t<_Int>;
        using base = dynamic_modint<Int>;
        return base::with_mod(p, [&](){
            base g = primitive_root(p);
            int64_t i = bpow(g, (p - 1) / 4).getr();
            int64_t q0 = 1, q1 = 0;
            int64_t r = i, m = p;
            // TODO: Use library contfrac?
            do {
                int64_t d = r / m;
                q0 = std::exchange(q1, q0 + d * q1);
                r = std::exchange(m, r % m);
            } while(q1 < p / q1);
            return gaussint<_Int>{q0, (base(i) * base(q0)).rem()};
        });
    }

    template<typename Int>
    big_vector<gaussint<Int>> two_squares_all(Int n) {
        if(n == 0) {
            return {0};
        }
        auto primes = factorize(n);
        big_map<Int, int> cnt;
        for(auto p: primes) {
            cnt[p]++;
        }
        big_vector<gaussint<Int>> res = {1};
        for(auto [p, c]: cnt) {
            if(p % 4 == 3 && c % 2) {
                return {};
            }
            if(p % 4 != 1) {
                auto mul = p == 2 ? bpow(gaussint<Int>(1, 1), c) : bpow(gaussint<Int>(p), c / 2);
                for(auto &value: res) {
                    value *= mul;
                }
            } else {
                big_vector<gaussint<Int>> nres;
                auto base = two_squares_prime_any(p);
                for(int i = 0; i <= c; i++) {
                    auto mul = bpow(base, i) * bpow(conj(base), c - i);
                    for(auto value: res) {
                        nres.push_back(value * mul);
                    }
                }
                res = std::move(nres);
            }
        }
        big_vector<gaussint<Int>> nres;
        for(auto p: res) {
            while(p.real() < 0 || p.imag() < 0) {
                p *= gaussint<Int>(0, 1);
            }
            nres.push_back(p);
            if(!p.real() || !p.imag()) {
                nres.emplace_back(p.imag(), p.real());
            }
        }
        return nres;
    }
}

#line 6 "verify/number_theory/two_squares.test.cpp"
#include <bits/stdc++.h>

using namespace std;
using namespace cp_algo::math;

void solve() {
	int64_t n;
	cin >> n;
	auto res = two_squares_all(n);
	cout << size(res) << "\n";
	for(auto p: res) {
		cout << p.real() << ' ' << p.imag() << "\n";
	}
}

signed main() {
    //freopen("input.txt", "r", stdin);
    ios::sync_with_stdio(0);
    cin.tie(0);
    int t = 1;
    cin >> t;
	while(t--) {
		solve();
	}
}

Test cases

Env Name Status Elapsed Memory
g++ big2_00 :heavy_check_mark: AC 43 ms 4 MB
g++ big2_01 :heavy_check_mark: AC 47 ms 4 MB
g++ big2_02 :heavy_check_mark: AC 50 ms 4 MB
g++ big_semiprime_random_00 :heavy_check_mark: AC 35 ms 4 MB
g++ big_semiprime_random_01 :heavy_check_mark: AC 36 ms 4 MB
g++ example_00 :heavy_check_mark: AC 2 ms 4 MB
g++ example_01 :heavy_check_mark: AC 2 ms 4 MB
g++ fixed_RNG_buster_00 :heavy_check_mark: AC 2 ms 4 MB
g++ max_00 :heavy_check_mark: AC 6 ms 4 MB
g++ max_01 :heavy_check_mark: AC 5 ms 4 MB
g++ max_02 :heavy_check_mark: AC 6 ms 4 MB
g++ max_03 :heavy_check_mark: AC 6 ms 4 MB
g++ max_04 :heavy_check_mark: AC 6 ms 4 MB
g++ num_of_solution_max_00 :heavy_check_mark: AC 65 ms 4 MB
g++ num_of_solution_max_01 :heavy_check_mark: AC 53 ms 4 MB
g++ power_of_2_3_5_13_00 :heavy_check_mark: AC 3 ms 4 MB
g++ power_of_2_3_5_13_01 :heavy_check_mark: AC 5 ms 4 MB
g++ power_of_2_3_5_13_02 :heavy_check_mark: AC 5 ms 4 MB
g++ power_of_2_3_5_13_03 :heavy_check_mark: AC 4 ms 4 MB
g++ prime_1_mod_4_00 :heavy_check_mark: AC 12 ms 4 MB
g++ prime_1_mod_4_01 :heavy_check_mark: AC 10 ms 4 MB
g++ prime_1_mod_4_02 :heavy_check_mark: AC 10 ms 4 MB
g++ prime_3_mod_4_00 :heavy_check_mark: AC 2 ms 4 MB
g++ random_00 :heavy_check_mark: AC 4 ms 4 MB
g++ random_01 :heavy_check_mark: AC 4 ms 4 MB
g++ random_02 :heavy_check_mark: AC 4 ms 4 MB
g++ random_03 :heavy_check_mark: AC 4 ms 4 MB
g++ random_04 :heavy_check_mark: AC 3 ms 4 MB
g++ small_00 :heavy_check_mark: AC 2 ms 4 MB
g++ small_01 :heavy_check_mark: AC 3 ms 4 MB
g++ small_02 :heavy_check_mark: AC 2 ms 4 MB
g++ small_03 :heavy_check_mark: AC 3 ms 4 MB
g++ small_04 :heavy_check_mark: AC 3 ms 4 MB
g++ smooth_00 :heavy_check_mark: AC 6 ms 4 MB
g++ smooth_01 :heavy_check_mark: AC 7 ms 4 MB
g++ smooth_02 :heavy_check_mark: AC 11 ms 4 MB
g++ smooth_03 :heavy_check_mark: AC 20 ms 4 MB
g++ smooth_04 :heavy_check_mark: AC 10 ms 4 MB
Back to top page