Number_Theory/Gauss-Integer.hpp

• View this file on GitHub
• Last update: 2023-04-03 20:21:35+08:00
• Include: #include "Number_Theory/Gauss-Integer.hpp"

Depends on

• Number_Theory/Binary-Gcd.hpp
• Number_Theory/Factorization.hpp
• Template/Template.hpp

Code

#pragma once

#include "../Template/Template.hpp"
#include "Factorization.hpp"

namespace Format_Fact {
using i128 = __int128;
struct G {
    i128 a, b;
    G(){};
    G(i128 a, i128 b) : a(a), b(b){};
    G friend operator+(const G &a, const G &b) {
        return {a.a + b.a, a.b + b.b};
    }
    G friend operator-(const G &a, const G &b) {
        return {a.a - b.a, a.b - b.b};
    }
    G friend operator*(const G &a, const G &b) {
        return {a.a * b.a - a.b * b.b, a.a * b.b + a.b * b.a};
    }
    bool operator==(const G &x) const { return x.a == a && x.b == b; }
    inline G operator*(const i128 &t) const { return {a * t, b * t}; }
    inline G operator/(const i128 &t) const { return {a / t, b / t}; }
    G friend operator/(const G &a, const G &b) {
        i128 len = b.a * b.a + b.b * b.b;
        G c = a * G(b.a, -b.b);
        auto div = [&](i128 a, i128 n) -> i128 {
            i128 now = (a % n + n) % n;
            return ((a - now) / n);
        };
        return {div(2 * c.a + len, 2 * len), div(2 * c.b + len, 2 * len)};
    }
    G power(i128 b) {
        G res(1, 0);
        G a = *this;
        for (; b; b /= 2, a = a * a) {
            if (b % 2) {
                res = res * a;
            }
        }
        return res;
    }
};

G solveprime(i128 p) {
    if (p == 2) return {1, 1};
    i128 t = 1;
    auto qpow = [](i128 a, i128 b, i128 p) {
        i128 res = 1;
        while (b) {
            if (b & 1) res = res * a % p;
            a = a * a % p;
            b = b / 2;
        }
        return res;
    };
    for (; qpow(t, (p - 1) / 2, p) != p - 1;) t++;
    i128 k = qpow(t, (p - 1) / 4, p);

    auto gcd = [&](auto &&self, G a, G b) -> G {
        if (b.a == 0 && b.b == 0)
            return a;
        else
            return self(self, b, a - (a / b) * b);
    };
    auto g = gcd(gcd, {p, 0}, {k, 1});

    if (g.a < 0) g.a = -g.a;
    if (g.b < 0) g.b = -g.b;
    if (g.a > g.b) std::swap(g.a, g.b);
    return g;
}
std::vector<G> solvecomposite(i128 n) {
    auto prm = Factor::factorSortedList<i128>(n);

    std::vector<G> v{{1, 0}};
    for (const auto &[p, tmp] : prm) {
        if (p % 4 == 1) {
            G A = solveprime(p);
            G B = {A.a, -A.b};
            auto now = A.power(2 * tmp);

            std::vector<G> res;
            for (i64 i = 0; i <= 2 * tmp; i++) {
                for (auto it : v) res.push_back(it * now);
                now = now * B / A;
            }
            swap(v, res);
        } else {
            G now(p, 0);
            now = now.power(tmp);
            for (auto &&it : v) it = it * now;
        }
    }
    for (auto &&[a, b] : v) {
        if (a < 0) a = -a;
        if (b < 0) b = -b;
    }
    std::sort(v.begin(), v.end(), [&](const G &a, const G &b) {
        return std::make_pair(a.a, a.b) < std::make_pair(b.a, b.b);
    });
    v.resize(unique(v.begin(), v.end()) - v.begin());

    std::vector<G> t;
    for (const auto &[a, b] : v)
        if (a != 0 && b != 0) t.emplace_back(a, b);
    return t;
}
}  // namespace Format_Fact

#line 2 "Number_Theory/Gauss-Integer.hpp"

#line 2 "Template/Template.hpp"

#include <bits/stdc++.h>

using i64 = std::int64_t;
#line 2 "Number_Theory/Factorization.hpp"

#line 2 "Number_Theory/Binary-Gcd.hpp"

#line 4 "Number_Theory/Binary-Gcd.hpp"

inline i64 binary_gcd(i64 a, i64 b) {
    if (a == 0 || b == 0) return a + b;
    char n = __builtin_ctzll(a);
    char m = __builtin_ctzll(b);
    a >>= n;
    b >>= m;
    n = std::min(n, m);
    while (a != b) {
        i64 d = a - b;
        char s = __builtin_ctzll(d);
        bool f = a > b;
        b = f ? b : a;
        a = (f ? d : -d) >> s;
    }
    return a << n;
}
#line 5 "Number_Theory/Factorization.hpp"

namespace Factor {
using u64 = std::uint64_t;

u64 modmul(u64 a, u64 b, u64 M) {
    i64 ret = a * b - M * u64(1.L / M * a * b);
    return ret + M * (ret < 0) - M * (ret >= (i64)M);
}

u64 modpow(u64 b, u64 e, u64 mod) {
    u64 ans = 1;
    for (; e; b = modmul(b, b, mod), e /= 2)
        if (e & 1) ans = modmul(ans, b, mod);
    return ans;
}

bool isPrime(u64 n) {
    if (n < 2 || n % 6 % 4 != 1) return (n | 1) == 3;
    std::vector<u64> A = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
    u64 s = __builtin_ctzll(n - 1), d = n >> s;
    for (u64 a : A) {  // ^ count trailing zeroes
        u64 p = modpow(a % n, d, n), i = s;
        while (p != 1 && p != n - 1 && a % n && i--) p = modmul(p, p, n);
        if (p != n - 1 && i != s) return 0;
    }
    return 1;
}
u64 pollard(u64 n) {
    auto f = [n](u64 x, u64 k) { return modmul(x, x, n) + k; };
    u64 x = 0, y = 0, t = 30, prd = 2, i = 1, q;
    while (t++ % 40 || binary_gcd(prd, n) == 1) {
        if (x == y) x = ++i, y = f(x, i);
        if ((q = modmul(prd, std::max(x, y) - std::min(x, y), n))) prd = q;
        x = f(x, i), y = f(f(y, i), i);
    }
    return std::gcd(prd, n);
}
std::vector<u64> factor(u64 n) {
    if (n == 1) return {};
    if (isPrime(n)) return {n};
    u64 x = pollard(n);
    auto l = factor(x), r = factor(n / x);
    l.insert(l.end(), r.begin(), r.end());
    return l;
}

template <class T = u64>
std::vector<std::pair<T, int>> factorSortedList(u64 n) {
    // \prid x_i^p_i
    auto fac = factor(n);
    std::sort(fac.begin(), fac.end());

    std::vector<std::pair<T, int>> lt;
    for (int i = 0, j; i < int(fac.size()); i = j) {
        j = i;
        while (j < static_cast<int>(fac.size()) && fac[i] == fac[j]) j++;
        lt.emplace_back(fac[i], j - i);
    }

    return lt;
}
}  // namespace Factor
#line 5 "Number_Theory/Gauss-Integer.hpp"

namespace Format_Fact {
using i128 = __int128;
struct G {
    i128 a, b;
    G(){};
    G(i128 a, i128 b) : a(a), b(b){};
    G friend operator+(const G &a, const G &b) {
        return {a.a + b.a, a.b + b.b};
    }
    G friend operator-(const G &a, const G &b) {
        return {a.a - b.a, a.b - b.b};
    }
    G friend operator*(const G &a, const G &b) {
        return {a.a * b.a - a.b * b.b, a.a * b.b + a.b * b.a};
    }
    bool operator==(const G &x) const { return x.a == a && x.b == b; }
    inline G operator*(const i128 &t) const { return {a * t, b * t}; }
    inline G operator/(const i128 &t) const { return {a / t, b / t}; }
    G friend operator/(const G &a, const G &b) {
        i128 len = b.a * b.a + b.b * b.b;
        G c = a * G(b.a, -b.b);
        auto div = [&](i128 a, i128 n) -> i128 {
            i128 now = (a % n + n) % n;
            return ((a - now) / n);
        };
        return {div(2 * c.a + len, 2 * len), div(2 * c.b + len, 2 * len)};
    }
    G power(i128 b) {
        G res(1, 0);
        G a = *this;
        for (; b; b /= 2, a = a * a) {
            if (b % 2) {
                res = res * a;
            }
        }
        return res;
    }
};

G solveprime(i128 p) {
    if (p == 2) return {1, 1};
    i128 t = 1;
    auto qpow = [](i128 a, i128 b, i128 p) {
        i128 res = 1;
        while (b) {
            if (b & 1) res = res * a % p;
            a = a * a % p;
            b = b / 2;
        }
        return res;
    };
    for (; qpow(t, (p - 1) / 2, p) != p - 1;) t++;
    i128 k = qpow(t, (p - 1) / 4, p);

    auto gcd = [&](auto &&self, G a, G b) -> G {
        if (b.a == 0 && b.b == 0)
            return a;
        else
            return self(self, b, a - (a / b) * b);
    };
    auto g = gcd(gcd, {p, 0}, {k, 1});

    if (g.a < 0) g.a = -g.a;
    if (g.b < 0) g.b = -g.b;
    if (g.a > g.b) std::swap(g.a, g.b);
    return g;
}
std::vector<G> solvecomposite(i128 n) {
    auto prm = Factor::factorSortedList<i128>(n);

    std::vector<G> v{{1, 0}};
    for (const auto &[p, tmp] : prm) {
        if (p % 4 == 1) {
            G A = solveprime(p);
            G B = {A.a, -A.b};
            auto now = A.power(2 * tmp);

            std::vector<G> res;
            for (i64 i = 0; i <= 2 * tmp; i++) {
                for (auto it : v) res.push_back(it * now);
                now = now * B / A;
            }
            swap(v, res);
        } else {
            G now(p, 0);
            now = now.power(tmp);
            for (auto &&it : v) it = it * now;
        }
    }
    for (auto &&[a, b] : v) {
        if (a < 0) a = -a;
        if (b < 0) b = -b;
    }
    std::sort(v.begin(), v.end(), [&](const G &a, const G &b) {
        return std::make_pair(a.a, a.b) < std::make_pair(b.a, b.b);
    });
    v.resize(unique(v.begin(), v.end()) - v.begin());

    std::vector<G> t;
    for (const auto &[a, b] : v)
        if (a != 0 && b != 0) t.emplace_back(a, b);
    return t;
}
}  // namespace Format_Fact

Back to top page