```cpp
#define __USE_MY_RAND__
#define __USE_MY_GCD__
#include <iostream>
#include <string>
#include <bitset>
#include <type_traits>
#include <cstdint>
#include <cstdlib>
#include <cstddef>
#include <algorithm>
#include <chrono>
#include <set>
#include <cstddef>
#include <cstdint>
#include <random>
#include <climits>
#include <cstring>
#include <unordered_map>
#if !__has_include(<bit>)
namespace std {
template <typename _Tp>
constexpr int __countr_zero(_Tp __x) noexcept {
constexpr auto _Nd = numeric_limits<_Tp>::digits;
if (__x == 0) return _Nd;
constexpr auto _Nd_ull = numeric_limits<unsigned long long>::digits;
constexpr auto _Nd_ul = numeric_limits<unsigned long>::digits;
constexpr auto _Nd_u = numeric_limits<unsigned>::digits;
if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_u) return __builtin_ctz(__x);
else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ul) return __builtin_ctzl(__x);
else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ull) return __builtin_ctzll(__x);
else { // (_Nd > _Nd_ull)
static_assert(_Nd <= (2 * _Nd_ull), "Maximum supported integer size is 128-bit");
constexpr auto __max_ull = numeric_limits<unsigned long long>::max();
unsigned long long __low = __x & __max_ull;
if (__low) return __builtin_ctzll(__low);
unsigned long long __high = __x >> _Nd_ull;
return __builtin_ctzll(__high) + _Nd_ull;
}
}
}
#else
#include <bit>
#endif
#if ( _WIN32 || __WIN32__ || _WIN64 || __WIN64__ )
# if !defined(_MSC_VER) || _MSC_VER>1400
# define NOMINMAX 1
# include <windows.h>
# else
# define WORD unsigned short
# include <unistd.h>
# endif
# include <io.h>
# define ON_WINDOWS
#else
# define WORD unsigned short
#endif
void set_text_color(
#if !(defined(ON_WINDOWS) && (!defined(_MSC_VER) || _MSC_VER>1400)) && defined(__GNUC__)
__attribute__((unused))
#endif
WORD color
)
{
#if defined(ON_WINDOWS) && (!defined(_MSC_VER) || _MSC_VER>1400)
static HANDLE handle = GetStdHandle(STD_OUTPUT_HANDLE);
SetConsoleTextAttribute(handle, color);
#endif
#if !defined(ON_WINDOWS) && defined(__CNUC__)
if (isatty(2))
{
switch (color)
{
case 0x0c:
fprintf(stderr, "\033[1;31m");
break;
case 0x0b:
fprintf(stderr, "\033[1;36m");
break;
case 0x0a:
fprintf(stderr, "\033[1;32m");
break;
case 0x0e:
fprintf(stderr, "\033[1;33m");
break;
case 0x07:
default:
fprintf(stderr, "\033[0m");
}
}
#endif
}
template <typename _Clock, typename = std::void_t<typename _Clock::rep, typename _Clock::period, typename _Clock::duration,
typename _Clock::time_point, decltype(_Clock::is_steady), decltype(_Clock::now())>>
struct __check_time_helper {
typename _Clock::time_point t;
double used;
void start() { t = _Clock::now(); }
void stop() { used += std::chrono::duration_cast<std::chrono::nanoseconds>
(_Clock::now() - t).count() / 1e6; }
~__check_time_helper() {
set_text_color(0x0c);
fprintf(stderr, "time used: %.2lfms\n", used);
set_text_color(0x07);
}
};
inline __uint128_t readu128() {
__uint128_t ret = 0;
char c = getchar();
while (!isdigit(c)) c = getchar();
while (isdigit(c)) ret = ret * 10 + (c & 15), c = getchar();
return ret;
}
inline void put128(__uint128_t x) {
static char buf[45];
static int curpos;
curpos = 0;
do {
buf[curpos++] = x % 10;
x /= 10;
} while (x);
while (curpos--) putchar(buf[curpos] | 48);
}
template <__uint128_t __mod, bool __isprime = true>
class DynamicModInt {
public:
using u32 = uint64_t;
using i64 = __int128_t;
using u64 = __uint128_t;
using m64 = DynamicModInt;
using value_type = u64;
static inline u64 mod() { return s_mod; }
static inline u64 get_primitive_root_prime() {
if (!s_isPrime) return 0;
u64 tmp[500] = {};
u64 cnt = 0;
const u64 phi = s_mod - 1;
u64 m = phi;
for (u64 i = 2; i * i <= m; ++i) {
if (m % i == 0) {
tmp[cnt++] = i;
do m /= i;
while (m % i == 0);
}
}
if (m != 1) tmp[cnt++] = m;
for (m64 res = 2;; res += 1) {
bool f = true;
for (int i = 0; i < cnt && f; ++i) f &= res.pow(phi / tmp[i]) != 1;
if (f) return u64(res);
}
}
constexpr DynamicModInt() : v_() {}
~DynamicModInt() = default;
template <typename T, typename std::enable_if<std::is_arithmetic<T>::value || std::is_same<T, i64>::value || std::is_same<T, u64>::value, int>::type = 0>
inline DynamicModInt(T v) : v_(reduce(mul(norm(v % i64(s_mod)), r2))) {}
constexpr DynamicModInt(const m64 &) = default;
inline u64 val() const { return reduce({0, v_}); }
template <typename T, typename std::enable_if<std::is_arithmetic<T>::value || std::is_same<T, i64>::value || std::is_same<T, u64>::value, int>::type = 0>
explicit constexpr operator T() const { return T(val()); }
inline m64 operator-() const {
m64 res;
res.v_ = (s_mod & -(v_ != 0)) - v_;
return res;
}
inline m64 inv_exgcd() const {
i64 x1 = 1, x3 = 0, a = val(), b = s_mod;
while (b != 0) {
i64 q = a / b, x1_old = x1, a_old = a;
x1 = x3, x3 = x1_old - x3 * q, a = b, b = a_old - b * q;
}
return m64(x1);
}
inline m64 inv_fermat() const { return pow(s_mod - 2); }
inline m64 inv() const { return s_isPrime ? inv_fermat() : inv_exgcd(); }
inline m64 &operator=(const m64 &) = default;
inline m64 &operator+=(const m64 &rhs) {
v_ += rhs.v_ - s_mod;
v_ += s_mod & -(v_ >> 127);
return *this;
}
inline m64 &operator-=(const m64 &rhs) {
v_ -= rhs.v_;
v_ += s_mod & -(v_ >> 127);
return *this;
}
inline m64 &operator*=(const m64 &rhs) {
v_ = reduce(mul(v_, rhs.v_));
return *this;
}
inline m64 &operator/=(const m64 &rhs) { return operator*=(rhs.inv()); }
friend inline m64 operator+(const m64 &lhs, const m64 &rhs) { return m64(lhs) += rhs; }
friend inline m64 operator-(const m64 &lhs, const m64 &rhs) { return m64(lhs) -= rhs; }
friend inline m64 operator*(const m64 &lhs, const m64 &rhs) { return m64(lhs) *= rhs; }
friend inline m64 operator/(const m64 &lhs, const m64 &rhs) { return m64(lhs) /= rhs; }
friend inline bool operator==(const m64 &lhs, const m64 &rhs) { return lhs.v_ == rhs.v_; }
friend inline bool operator!=(const m64 &lhs, const m64 &rhs) { return lhs.v_ != rhs.v_; }
template <typename _Istream>
friend inline _Istream &operator>>(_Istream &is, m64 &rhs) {
i64 x;
is >> x;
rhs = m64(x);
return is;
}
template <typename _Ostream>
friend _Ostream &operator<<(_Ostream &os, const m64 &rhs) { return os << rhs.val(); }
inline m64 pow(u64 y) const {
m64 res(1), x(*this);
for (; y != 0; y >>= 1, x *= x) if (y & 1) res *= x;
return res;
}
static inline void set_mod(u64 mod, bool prime = true) {
s_mod = mod;
s_isPrime = prime;
r = get_r();
r2 = get_r2();
}
static inline u64 init(u64 val) { return reduce(mul(val, r2)); }
private:
static constexpr std::pair<u64, u64> mul(u64 x, u64 y) {
u64 a = x >> 64, b = u32(x), c = y >> 64, d = u32(y), ad = a * d, bc = b * c;
return {a * c + (ad >> 64) + (bc >> 64) + (((ad & ~UINT64_C(0)) + (bc & ~UINT64_C(0)) + (b * d >> 64)) >> 64), x * y};
}
static constexpr u64 mulh(u64 x, u64 y) {
u64 a = x >> 64, b = u32(x), c = y >> 64, d = u32(y), ad = a * d, bc = b * c;
return a * c + (ad >> 64) + (bc >> 64) + (((ad & ~UINT64_C(0)) + (bc & ~UINT64_C(0)) + (b * d >> 64)) >> 64);
}
static inline u64 get_r() {
u64 two = 2, iv = s_mod * (two - s_mod * s_mod);
iv *= two - s_mod * iv;
iv *= two - s_mod * iv;
iv *= two - s_mod * iv;
iv *= two - s_mod * iv;
return iv * (two - s_mod * iv);
}
static inline u64 get_r2() {
u64 iv = -u64(s_mod) % s_mod;
for (int i = 0; i != 128; ++i) ((iv <<= 1) >= s_mod) && (iv -= s_mod);
return iv;
}
static constexpr u64 reduce(const std::pair<u64, u64> &x) {
u64 res = x.first - mulh(x.second * r, s_mod);
return res + (s_mod & -(res >> 127));
}
static constexpr u64 norm(i64 x) { return x + (s_mod & -(x < 0)); }
u64 v_;
static inline u64 s_mod = __mod;
static inline u64 s_isPrime = __isprime;
static inline u64 r = get_r();
static inline u64 r2 = get_r2();
};
using mint = DynamicModInt<3, true>;
inline bool prime(const __uint128_t& x) {
static constexpr __uint128_t first_prime[25] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97};
if (x <= 100) return std::binary_search(first_prime, first_prime + 25, x);
std::vector<__uint128_t> bases;
mint::set_mod(x, false);
const mint one(1), none(-1);
#define DEFINE(threshold, cnt) if (x >= threshold) bases = std::vector<__uint128_t>(first_prime, first_prime + cnt)
DEFINE(1543267864443420616877677640751301, 20);
else DEFINE(564132928021909221014087501701, 18);
else DEFINE(59276361075595573263446330101, 16);
else DEFINE(6003094289670105800312596501, 15);
else DEFINE(3317044064679887385961981, 14);
else DEFINE(318665857834031151167461, 13);
else DEFINE(3825123056546413051, 12);
else DEFINE(341550071728321, 9);
else DEFINE(3474749660383, 7);
else DEFINE(2152302898747, 6);
else DEFINE(4759123141, 5);
else bases = {2, 7, 41, 61, 325};
#undef DEFINE
__uint128_t y(x - 1);
const unsigned iend = std::__countr_zero(y); // __builtin_ctzl
y >>= iend;
for (mint base : bases) {
mint z(base.pow(y));
if (z.val() == 1) continue;
unsigned i = 0;
for (; i ^ iend; i++) {
if (z == none) break;
z *= z;
}
if (i == iend) return false;
}
return true;
}
#ifdef __USE_MY_RAND__
unsigned long long seed = 5489ull;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wstrict-aliasing"
unsigned long long u64_hasher(unsigned long long x) {
double tmp = x;
return *(unsigned long long*)&tmp;
}
#pragma GCC diagnostic pop
unsigned long long rng() { return seed += u64_hasher(seed); }
#else
static std::mt19937_64 rng;
#endif
#ifdef __USE_MY_GCD__
inline __uint128_t gcd(__uint128_t a, __uint128_t b) {
if (!a || !b) return a | b;
int i = std::__countr_zero(a), j = std::__countr_zero(b), k = std::min(i, j);
a >>= i;
b >>= j;
while (true) {
if (a < b) std::swap(a, b);
if (!(a -= b)) break;
a >>= std::__countr_zero(a);
}
return b << k;
}
#else
inline __uint128_t gcd(__uint128_t a, __uint128_t b) { return std::__gcd(a, b); }
#endif
inline __uint128_t rho(__uint128_t value) {
mint::set_mod(value, false);
mint x, y, z, c;
uint64_t i, j;
while (true) {
y = x = rng();
z = 1;
c = rng();
i = 0;
j = 1;
while (++i) {
y = y * y + c;
z *= (x.val() > y.val()) ? (x - y) : (y - x);
if (x == y || !z.val()) break;
if (!(i & 127) || i == j) {
if (__uint128_t g = gcd(z.val(), value); g > 1) return g;
if (i == j) {
x = y;
j <<= 1;
}
}
}
}
}
__uint128_t squfof(__uint128_t value) {
__uint128_t s = sqrtl(1.0l * value) + 1e-6;
while (s * s > value) s--;
if (s * s == value) return s;
__int128 d, po, p, p_prev, q, q_prev, q_, b, r;
long long l, b_, i;
int k = 0;
l = 2 * sqrtl(2.0l * s);
b_ = 3 * l;
if (b_ > 20000000) return rho(value);
while (++k) {
d = k * value;
po = p_prev = p = sqrtl(1.0l * d);
q_prev = 1;
q = d - po * po;
while (q < 0) {
po--;
p_prev--;
p--;
q = d - po * po;
}
if (!q) return gcd(value, k);
for (i = 2; i < b_; i++) {
b = (po + p) / q;
p = b * q - p;
q_ = q;
q = q_prev + b * (p_prev - p);
r = sqrtl(1.0 * q) + 1e-6;
if (!(i & 1) && r * r == q) break;
q_prev = q_;
p_prev = p;
}
if (i >= b_) continue;
b = (po - p) / r;
p_prev = p = b * r + p;
q_prev = r;
q = (d - p_prev * p_prev) / q_prev;
i = 0;
do {
b = (po + p) / q;
p_prev = p;
p = b * q - p;
q_ = q;
q = q_prev + b * (p_prev - p);
q_prev = q_;
i++;
} while (p != p_prev && i < b_);
if (i >= b_) continue;
r = gcd(value, q_prev);
if (r != 1) return r;
}
return 0;
}
template <typename _ModType>
class barrett {
_ModType _M_mod;
int _M_l;
__uint128_t _M_inv;
static constexpr auto _Nd = std::numeric_limits<_ModType>::digits;
static constexpr auto _Nd_ull = std::numeric_limits<unsigned long long>::digits;
static constexpr auto _Nd_ul = std::numeric_limits<unsigned long>::digits;
static constexpr auto _Nd_u = std::numeric_limits<unsigned>::digits;
public:
static_assert(std::is_integral<_ModType>::value && std::is_unsigned<_ModType>::value, "_ModType must be an unsigned integral type");
constexpr barrett() = default;
constexpr barrett(const _ModType& mod)
: _M_mod(mod),
_M_l(_Nd_ull + _Nd - __builtin_clzll((unsigned long long)mod) - ((mod & (mod - 1)) != 0)),
_M_inv(((unsigned __int128)1 << _M_l) / mod + ((mod & (mod - 1)) != 0)) {}
constexpr void set_mod(const _ModType& mod) {
_M_mod = mod;
_M_l = _Nd_ull + _Nd - __builtin_clzll((unsigned long long)mod) - ((mod & (mod - 1)) != 0);
_M_inv = ((unsigned __int128)1 << _M_l) / mod + ((mod & (mod - 1)) != 0);
}
constexpr _ModType mod() const { return _M_mod; }
constexpr _ModType mod(const _ModType& x) const {
__uint128_t tmp = _M_inv * x;
_ModType ret = x - (tmp >> _M_l) * _M_mod;
//while (ret > _M_mod) ret -= _M_mod;
return ret;
}
};
namespace qs {
const size_t L = 510;
uint32_t pcnt, root[L];
mint primes[L];
double log_primes[L];
struct word {
std::bitset<L> mask;
mint lhs, rhs;
size_t bit;
word() : mask(), lhs(), rhs(), bit(L) {}
void merge(const word& rOther) {
lhs *= rOther.lhs;
rhs *= rOther.rhs;
const std::bitset<L> cons(mask & rOther.mask);
for (size_t pos = cons._Find_first(); pos < L; pos = cons._Find_next(pos)) rhs *= primes[pos];
mask ^= rOther.mask;
bit = mask._Find_first();
}
};
std::vector<word> smooth;
std::unordered_map<long long, word> data;
__uint128_t insert(word& o) {
size_t x;
for (x = 0; x < smooth.size(); x++) {
if (smooth[x].bit > o.bit) break;
if (smooth[x].bit == o.bit) o.merge(smooth[x]);
}
if (o.bit == L) {
__uint128_t g(gcd((o.lhs + o.rhs).val(), mint::mod()));
if (g != 1 && g != mint::mod()) return g;
}
smooth.insert(smooth.begin() + x, o);
return 0;
}
__int128 try_insert(mint x, __int128 y) {
word ins;
ins.lhs = x;
ins.rhs = 1;
for (size_t k = 1; k <= pcnt; k++) {
size_t cnt = 0;
while (y % (__int128)primes[k].val() == 0) y /= (__int128)primes[k].val(), ++cnt;
if (cnt) {
ins.mask[k] = cnt & 1;
ins.rhs *= primes[k].pow(cnt >> 1);
}
}
ins.bit = ins.mask._Find_first();
if (y != 1) {
if (data.count(y)) {
ins.merge(data[y]);
ins.rhs *= mint(y);
y = 1;
}
else data[y] = ins;
}
if (y == 1) return insert(ins);
return 0;
}
void append(uint32_t p) {
++pcnt;
primes[pcnt] = p;
log_primes[pcnt] = log(1.0 * p);
}
uint32_t fastpow(__uint128_t base, uint32_t power, uint32_t mod) {
barrett<uint32_t> brt(mod);
base %= mod;
uint32_t ans = 1;
for (; power; base = brt.mod(base * base), power >>= 1) if (power & 1) ans = brt.mod(ans * base);
return ans;
}
void init() {
uint32_t i, j;
uint32_t B = pow(log(1.0 * mint::mod()), 2) * 0.6;
std::vector<char> mark((B >> 1) + 5);
for (i = 3; i * i <= B; i += 2) if (!mark[i >> 1]) for (j = i * i; j <= B; j += i << 1) mark[j >> 1] = true;
for (append(2u), i = 3; i <= B; i += 2) if (!mark[i >> 1]) if (fastpow(mint::mod(), i >> 1, i) == 1) append(i);
barrett<uint64_t> brt;
for (i = 1; i <= pcnt; i++) for (brt.set_mod(primes[i].val()), j = mint::mod() % primes[i].val(); brt.mod(1ull * root[i] * root[i]) != j; ++root[i]);
}
__uint128_t next_prime(__uint128_t x) {
x += ~x & 1;
while (!prime(x += 2));
return x;
}
__uint128_t quadratic_sieve(__int128 value) {
if (__uint128_t tmp = sqrtl(1.0 * value) + 1e-6; tmp * tmp == value) return tmp;
if (value <= 1024) return squfof(value);
__check_time_helper<std::chrono::steady_clock> _Helper;
_Helper.start();
pcnt = 0;
data.clear();
smooth.clear();
mint::set_mod(value, false);
init();
uint32_t M = pcnt << 10;
uint64_t D = sqrtl(sqrtl(2.0L * value) / M);
double bound = log(M * sqrtl(0.5L * value)) * 0.75;
barrett<uint64_t> brt;
while (true) {
do D = next_prime(D);
while ((D & 3) == 1);
mint::set_mod(D, true);
uint64_t y0;
mint y1;
{
mint tmp_value(value);
if (!tmp_value.val()) return D;
mint tmp_y0 = tmp_value.pow((D + 1) >> 2);
if (tmp_y0 * tmp_y0 != tmp_value) continue;
y0 = tmp_y0.val();
y1 = (value - y0 * y0) / D;
y1 = y1 * (D / 2 + 1) / tmp_y0;
}
__int128 A(D);
A *= A;
__int128 B(y1.val() * D + y0);
__int128 C((B * B - value) / A);
mint::set_mod(value, false);
__int128 E = (mint(B) / D).val();
std::vector<double> pos(M), neg(M);
for (uint32_t x = 1; x <= pcnt; x++) {
uint32_t p = primes[x].val();
brt.set_mod(p);
uint32_t s = brt.mod(A);
uint32_t a = s;
uint32_t inv_a = 1;
if (!a) continue;
while (a > 1) inv_a = brt.mod(1ull * (p - inv_a) * (p / a)), a = p % a;
uint32_t u = brt.mod(1ull * (p - brt.mod(B)) * inv_a);
uint32_t v = brt.mod(1ull * root[x] * inv_a);
uint32_t r1 = brt.mod(u + v), r2 = brt.mod(u + p - v);
for (uint32_t su = 0; su < M; su += p) {
if (su + r1 < M) pos[su + r1] += log_primes[x];
if (su + r2 < M) pos[su + r2] += log_primes[x];
if (su >= r1) neg[su - r1] += log_primes[x];
if (su >= r2) neg[su - r2] += log_primes[x];
}
}
for (uint32_t x = 0; x < M; ++x) {
if (pos[x] > bound) if (__uint128_t tmp(try_insert(E + D * x, A * x * x + 2 * B * x + C)); tmp) return _Helper.stop(), tmp;
if (neg[x] > bound) if (__uint128_t tmp(try_insert(E - D * x, A * x * x - 2 * B * x + C)); tmp) return _Helper.stop(), tmp;
}
}
return value;
}
}
inline auto crack(__uint128_t value) {
struct _RetType {
__uint128_t base;
mutable uint64_t power;
bool operator<(const _RetType& rhs) const { return base < rhs.base; }
};
std::set<_RetType> ret;
if (value <= 1) return ret;
if (value == 2) {
ret.insert({2, 1});
return ret;
}
_RetType node{2, 0};
while (value % 2 == 0) ++(node.power), value >>= 1;
if (node.power) ret.insert(node);
if (value <= 1) return ret;
if (prime(value)) {
ret.insert({value, 1ul});
return ret;
}
__uint128_t factor = qs::quadratic_sieve(value);
auto dfs = [&](auto self, __uint128_t v) {
if (v <= 1) return;
if (!prime(v)) {
__uint128_t tmp = rho(v);
self(self, tmp);
self(self, v / tmp);
}
else {
auto pos = ret.find({v, 0});
if (pos == ret.end()) ret.insert({v, 1ul});
else ++(pos->power);
}
};
dfs(dfs, factor);
dfs(dfs, value / factor);
return ret;
}
int main() {
int t;
scanf("%d", &t);
#ifdef __USE_MY_RAND__
seed = std::random_device{}();
seed *= std::random_device{}();
#else
rng.seed(std::random_device{}());
#endif
__check_time_helper<std::chrono::steady_clock> _Helper;
while (t--) {
__uint128_t x = readu128();
_Helper.start();
auto map = crack(x);
_Helper.stop();
auto it = --map.end();
if (map.size() == 1 && it->base == x) puts("Prime");
else put128(it->base), putchar('\n');
}
}
```
by Ruiqun2009 @ 2023-03-11 23:17:55
ruiqun 老师 /bx /bx
by 一扶苏一 @ 2023-03-12 00:39:12
ruiqun 老师 /bx /bx
by poly @ 2023-03-12 07:14:34
此贴结
原因：求 `root` 的时候没有使用二次剩余
by Ruiqun2009 @ 2023-03-12 21:02:31