#line 1 "verify/verify-unit-test/enumerate-convex.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"
//
#line 2 "template/template.hpp"
using namespace std;
// intrinstic
#include <immintrin.h>
#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <cctype>
#include <cfenv>
#include <cfloat>
#include <chrono>
#include <cinttypes>
#include <climits>
#include <cmath>
#include <complex>
#include <cstdarg>
#include <cstddef>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <deque>
#include <fstream>
#include <functional>
#include <initializer_list>
#include <iomanip>
#include <ios>
#include <iostream>
#include <istream>
#include <iterator>
#include <limits>
#include <list>
#include <map>
#include <memory>
#include <new>
#include <numeric>
#include <ostream>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <stack>
#include <streambuf>
#include <string>
#include <tuple>
#include <type_traits>
#include <typeinfo>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
// utility
#line 1 "template/util.hpp"
namespace Nyaan {
using ll = long long;
using i64 = long long;
using u64 = unsigned long long;
using i128 = __int128_t;
using u128 = __uint128_t;
template <typename T>
using V = vector<T>;
template <typename T>
using VV = vector<vector<T>>;
using vi = vector<int>;
using vl = vector<long long>;
using vd = V<double>;
using vs = V<string>;
using vvi = vector<vector<int>>;
using vvl = vector<vector<long long>>;
template <typename T>
using minpq = priority_queue<T, vector<T>, greater<T>>;
template <typename T, typename U>
struct P : pair<T, U> {
template <typename... Args>
P(Args... args) : pair<T, U>(args...) {}
using pair<T, U>::first;
using pair<T, U>::second;
P &operator+=(const P &r) {
first += r.first;
second += r.second;
return *this;
}
P &operator-=(const P &r) {
first -= r.first;
second -= r.second;
return *this;
}
P &operator*=(const P &r) {
first *= r.first;
second *= r.second;
return *this;
}
template <typename S>
P &operator*=(const S &r) {
first *= r, second *= r;
return *this;
}
P operator+(const P &r) const { return P(*this) += r; }
P operator-(const P &r) const { return P(*this) -= r; }
P operator*(const P &r) const { return P(*this) *= r; }
template <typename S>
P operator*(const S &r) const {
return P(*this) *= r;
}
P operator-() const { return P{-first, -second}; }
};
using pl = P<ll, ll>;
using pi = P<int, int>;
using vp = V<pl>;
constexpr int inf = 1001001001;
constexpr long long infLL = 4004004004004004004LL;
template <typename T>
int sz(const T &t) {
return t.size();
}
template <typename T, typename U>
inline bool amin(T &x, U y) {
return (y < x) ? (x = y, true) : false;
}
template <typename T, typename U>
inline bool amax(T &x, U y) {
return (x < y) ? (x = y, true) : false;
}
template <typename T>
inline T Max(const vector<T> &v) {
return *max_element(begin(v), end(v));
}
template <typename T>
inline T Min(const vector<T> &v) {
return *min_element(begin(v), end(v));
}
template <typename T>
inline long long Sum(const vector<T> &v) {
return accumulate(begin(v), end(v), 0LL);
}
template <typename T>
int lb(const vector<T> &v, const T &a) {
return lower_bound(begin(v), end(v), a) - begin(v);
}
template <typename T>
int ub(const vector<T> &v, const T &a) {
return upper_bound(begin(v), end(v), a) - begin(v);
}
constexpr long long TEN(int n) {
long long ret = 1, x = 10;
for (; n; x *= x, n >>= 1) ret *= (n & 1 ? x : 1);
return ret;
}
template <typename T, typename U>
pair<T, U> mkp(const T &t, const U &u) {
return make_pair(t, u);
}
template <typename T>
vector<T> mkrui(const vector<T> &v, bool rev = false) {
vector<T> ret(v.size() + 1);
if (rev) {
for (int i = int(v.size()) - 1; i >= 0; i--) ret[i] = v[i] + ret[i + 1];
} else {
for (int i = 0; i < int(v.size()); i++) ret[i + 1] = ret[i] + v[i];
}
return ret;
};
template <typename T>
vector<T> mkuni(const vector<T> &v) {
vector<T> ret(v);
sort(ret.begin(), ret.end());
ret.erase(unique(ret.begin(), ret.end()), ret.end());
return ret;
}
template <typename F>
vector<int> mkord(int N, F f) {
vector<int> ord(N);
iota(begin(ord), end(ord), 0);
sort(begin(ord), end(ord), f);
return ord;
}
template <typename T>
vector<int> mkinv(vector<T> &v) {
int max_val = *max_element(begin(v), end(v));
vector<int> inv(max_val + 1, -1);
for (int i = 0; i < (int)v.size(); i++) inv[v[i]] = i;
return inv;
}
vector<int> mkiota(int n) {
vector<int> ret(n);
iota(begin(ret), end(ret), 0);
return ret;
}
template <typename T>
T mkrev(const T &v) {
T w{v};
reverse(begin(w), end(w));
return w;
}
template <typename T>
bool nxp(T &v) {
return next_permutation(begin(v), end(v));
}
// 返り値の型は入力の T に依存
// i 要素目 : [0, a[i])
template <typename T>
vector<vector<T>> product(const vector<T> &a) {
vector<vector<T>> ret;
vector<T> v;
auto dfs = [&](auto rc, int i) -> void {
if (i == (int)a.size()) {
ret.push_back(v);
return;
}
for (int j = 0; j < a[i]; j++) v.push_back(j), rc(rc, i + 1), v.pop_back();
};
dfs(dfs, 0);
return ret;
}
// F : function(void(T&)), mod を取る操作
// T : 整数型のときはオーバーフローに注意する
template <typename T>
T Power(T a, long long n, const T &I, const function<void(T &)> &f) {
T res = I;
for (; n; f(a = a * a), n >>= 1) {
if (n & 1) f(res = res * a);
}
return res;
}
// T : 整数型のときはオーバーフローに注意する
template <typename T>
T Power(T a, long long n, const T &I = T{1}) {
return Power(a, n, I, function<void(T &)>{[](T &) -> void {}});
}
template <typename T>
T Rev(const T &v) {
T res = v;
reverse(begin(res), end(res));
return res;
}
template <typename T>
vector<T> Transpose(const vector<T> &v) {
using U = typename T::value_type;
if(v.empty()) return {};
int H = v.size(), W = v[0].size();
vector res(W, T(H, U{}));
for (int i = 0; i < H; i++) {
for (int j = 0; j < W; j++) {
res[j][i] = v[i][j];
}
}
return res;
}
template <typename T>
vector<T> Rotate(const vector<T> &v, int clockwise = true) {
using U = typename T::value_type;
int H = v.size(), W = v[0].size();
vector res(W, T(H, U{}));
for (int i = 0; i < H; i++) {
for (int j = 0; j < W; j++) {
if (clockwise) {
res[W - 1 - j][i] = v[i][j];
} else {
res[j][H - 1 - i] = v[i][j];
}
}
}
return res;
}
} // namespace Nyaan
#line 58 "template/template.hpp"
// bit operation
#line 1 "template/bitop.hpp"
namespace Nyaan {
__attribute__((target("popcnt"))) inline int popcnt(const u64 &a) {
return __builtin_popcountll(a);
}
inline int lsb(const u64 &a) { return a ? __builtin_ctzll(a) : 64; }
inline int ctz(const u64 &a) { return a ? __builtin_ctzll(a) : 64; }
inline int msb(const u64 &a) { return a ? 63 - __builtin_clzll(a) : -1; }
template <typename T>
inline int gbit(const T &a, int i) {
return (a >> i) & 1;
}
template <typename T>
inline void sbit(T &a, int i, bool b) {
if (gbit(a, i) != b) a ^= T(1) << i;
}
constexpr long long PW(int n) { return 1LL << n; }
constexpr long long MSK(int n) { return (1LL << n) - 1; }
} // namespace Nyaan
#line 61 "template/template.hpp"
// inout
#line 1 "template/inout.hpp"
namespace Nyaan {
template <typename T, typename U>
ostream &operator<<(ostream &os, const pair<T, U> &p) {
os << p.first << " " << p.second;
return os;
}
template <typename T, typename U>
istream &operator>>(istream &is, pair<T, U> &p) {
is >> p.first >> p.second;
return is;
}
template <typename T>
ostream &operator<<(ostream &os, const vector<T> &v) {
int s = (int)v.size();
for (int i = 0; i < s; i++) os << (i ? " " : "") << v[i];
return os;
}
template <typename T>
istream &operator>>(istream &is, vector<T> &v) {
for (auto &x : v) is >> x;
return is;
}
istream &operator>>(istream &is, __int128_t &x) {
string S;
is >> S;
x = 0;
int flag = 0;
for (auto &c : S) {
if (c == '-') {
flag = true;
continue;
}
x *= 10;
x += c - '0';
}
if (flag) x = -x;
return is;
}
istream &operator>>(istream &is, __uint128_t &x) {
string S;
is >> S;
x = 0;
for (auto &c : S) {
x *= 10;
x += c - '0';
}
return is;
}
ostream &operator<<(ostream &os, __int128_t x) {
if (x == 0) return os << 0;
if (x < 0) os << '-', x = -x;
string S;
while (x) S.push_back('0' + x % 10), x /= 10;
reverse(begin(S), end(S));
return os << S;
}
ostream &operator<<(ostream &os, __uint128_t x) {
if (x == 0) return os << 0;
string S;
while (x) S.push_back('0' + x % 10), x /= 10;
reverse(begin(S), end(S));
return os << S;
}
void in() {}
template <typename T, class... U>
void in(T &t, U &...u) {
cin >> t;
in(u...);
}
void out() { cout << "\n"; }
template <typename T, class... U, char sep = ' '>
void out(const T &t, const U &...u) {
cout << t;
if (sizeof...(u)) cout << sep;
out(u...);
}
struct IoSetupNya {
IoSetupNya() {
cin.tie(nullptr);
ios::sync_with_stdio(false);
cout << fixed << setprecision(15);
cerr << fixed << setprecision(7);
}
} iosetupnya;
} // namespace Nyaan
#line 64 "template/template.hpp"
// debug
#line 1 "template/debug.hpp"
namespace DebugImpl {
template <typename U, typename = void>
struct is_specialize : false_type {};
template <typename U>
struct is_specialize<
U, typename conditional<false, typename U::iterator, void>::type>
: true_type {};
template <typename U>
struct is_specialize<
U, typename conditional<false, decltype(U::first), void>::type>
: true_type {};
template <typename U>
struct is_specialize<U, enable_if_t<is_integral<U>::value, void>> : true_type {
};
void dump(const char& t) { cerr << t; }
void dump(const string& t) { cerr << t; }
void dump(const bool& t) { cerr << (t ? "true" : "false"); }
void dump(__int128_t t) {
if (t == 0) cerr << 0;
if (t < 0) cerr << '-', t = -t;
string S;
while (t) S.push_back('0' + t % 10), t /= 10;
reverse(begin(S), end(S));
cerr << S;
}
void dump(__uint128_t t) {
if (t == 0) cerr << 0;
string S;
while (t) S.push_back('0' + t % 10), t /= 10;
reverse(begin(S), end(S));
cerr << S;
}
template <typename U,
enable_if_t<!is_specialize<U>::value, nullptr_t> = nullptr>
void dump(const U& t) {
cerr << t;
}
template <typename T>
void dump(const T& t, enable_if_t<is_integral<T>::value>* = nullptr) {
string res;
if (t == Nyaan::inf) res = "inf";
if constexpr (is_signed<T>::value) {
if (t == -Nyaan::inf) res = "-inf";
}
if constexpr (sizeof(T) == 8) {
if (t == Nyaan::infLL) res = "inf";
if constexpr (is_signed<T>::value) {
if (t == -Nyaan::infLL) res = "-inf";
}
}
if (res.empty()) res = to_string(t);
cerr << res;
}
template <typename T, typename U>
void dump(const pair<T, U>&);
template <typename T>
void dump(const pair<T*, int>&);
template <typename T>
void dump(const T& t,
enable_if_t<!is_void<typename T::iterator>::value>* = nullptr) {
cerr << "[ ";
for (auto it = t.begin(); it != t.end();) {
dump(*it);
cerr << (++it == t.end() ? "" : ", ");
}
cerr << " ]";
}
template <typename T, typename U>
void dump(const pair<T, U>& t) {
cerr << "( ";
dump(t.first);
cerr << ", ";
dump(t.second);
cerr << " )";
}
template <typename T>
void dump(const pair<T*, int>& t) {
cerr << "[ ";
for (int i = 0; i < t.second; i++) {
dump(t.first[i]);
cerr << (i == t.second - 1 ? "" : ", ");
}
cerr << " ]";
}
void trace() { cerr << endl; }
template <typename Head, typename... Tail>
void trace(Head&& head, Tail&&... tail) {
cerr << " ";
dump(head);
if (sizeof...(tail) != 0) cerr << ",";
trace(std::forward<Tail>(tail)...);
}
} // namespace DebugImpl
#ifdef NyaanDebug
#define trc(...) \
do { \
cerr << "## " << #__VA_ARGS__ << " = "; \
DebugImpl::trace(__VA_ARGS__); \
} while (0)
#else
#define trc(...) (void(0))
#endif
#ifdef NyaanLocal
#define trc2(...) \
do { \
cerr << "## " << #__VA_ARGS__ << " = "; \
DebugImpl::trace(__VA_ARGS__); \
} while (0)
#else
#define trc2(...) (void(0))
#endif
#line 67 "template/template.hpp"
// macro
#line 1 "template/macro.hpp"
#define each(x, v) for (auto&& x : v)
#define each2(x, y, v) for (auto&& [x, y] : v)
#define all(v) (v).begin(), (v).end()
#define rep(i, N) for (long long i = 0; i < (long long)(N); i++)
#define repr(i, N) for (long long i = (long long)(N)-1; i >= 0; i--)
#define rep1(i, N) for (long long i = 1; i <= (long long)(N); i++)
#define repr1(i, N) for (long long i = (N); (long long)(i) > 0; i--)
#define reg(i, a, b) for (long long i = (a); i < (b); i++)
#define regr(i, a, b) for (long long i = (b)-1; i >= (a); i--)
#define fi first
#define se second
#define ini(...) \
int __VA_ARGS__; \
in(__VA_ARGS__)
#define inl(...) \
long long __VA_ARGS__; \
in(__VA_ARGS__)
#define ins(...) \
string __VA_ARGS__; \
in(__VA_ARGS__)
#define in2(s, t) \
for (int i = 0; i < (int)s.size(); i++) { \
in(s[i], t[i]); \
}
#define in3(s, t, u) \
for (int i = 0; i < (int)s.size(); i++) { \
in(s[i], t[i], u[i]); \
}
#define in4(s, t, u, v) \
for (int i = 0; i < (int)s.size(); i++) { \
in(s[i], t[i], u[i], v[i]); \
}
#define die(...) \
do { \
Nyaan::out(__VA_ARGS__); \
return; \
} while (0)
#line 70 "template/template.hpp"
namespace Nyaan {
void solve();
}
int main() { Nyaan::solve(); }
#line 4 "verify/verify-unit-test/enumerate-convex.test.cpp"
//
#line 2 "math/enumerate-convex.hpp"
#line 6 "math/enumerate-convex.hpp"
using namespace std;
#line 2 "math/stern-brocot-tree.hpp"
#line 6 "math/stern-brocot-tree.hpp"
using namespace std;
// x / y (x > 0, y > 0) を管理、デフォルトで 1 / 1
// 入力が互いに素でない場合は gcd を取って格納
// seq : (1, 1) から (x, y) へのパス。右の子が正/左の子が負
template <typename Int>
struct SternBrocotTreeNode {
using Node = SternBrocotTreeNode;
Int lx, ly, x, y, rx, ry;
vector<Int> seq;
SternBrocotTreeNode() : lx(0), ly(1), x(1), y(1), rx(1), ry(0) {}
SternBrocotTreeNode(Int X, Int Y) : SternBrocotTreeNode() {
assert(1 <= X && 1 <= Y);
Int g = gcd(X, Y);
X /= g, Y /= g;
while (min(X, Y) > 0) {
if (X > Y) {
Int d = X / Y;
X -= d * Y;
go_right(d - (X == 0 ? 1 : 0));
} else {
Int d = Y / X;
Y -= d * X;
go_left(d - (Y == 0 ? 1 : 0));
}
}
}
SternBrocotTreeNode(const pair<Int, Int> &xy)
: SternBrocotTreeNode(xy.first, xy.second) {}
SternBrocotTreeNode(const vector<Int> &_seq) : SternBrocotTreeNode() {
for (const Int &d : _seq) {
assert(d != 0);
if (d > 0) go_right(d);
if (d < 0) go_left(d);
}
assert(seq == _seq);
}
// pair<Int, Int> 型で分数を get
pair<Int, Int> get() const { return make_pair(x, y); }
// 区間の下限
pair<Int, Int> lower_bound() const { return make_pair(lx, ly); }
// 区間の上限
pair<Int, Int> upper_bound() const { return make_pair(rx, ry); }
// 根からの深さ
Int depth() const {
Int res = 0;
for (auto &s : seq) res += abs(s);
return res;
}
// 左の子に d 進む
void go_left(Int d = 1) {
if (d <= 0) return;
if (seq.empty() or seq.back() > 0) seq.push_back(0);
seq.back() -= d;
rx += lx * d, ry += ly * d;
x = rx + lx, y = ry + ly;
}
// 右の子に d 進む
void go_right(Int d = 1) {
if (d <= 0) return;
if (seq.empty() or seq.back() < 0) seq.push_back(0);
seq.back() += d;
lx += rx * d, ly += ry * d;
x = rx + lx, y = ry + ly;
}
// 親の方向に d 進む
// d 進めたら true, 進めなかったら false を返す
bool go_parent(Int d = 1) {
if (d <= 0) return true;
while (d != 0) {
if (seq.empty()) return false;
Int d2 = min(d, abs(seq.back()));
if (seq.back() > 0) {
x -= rx * d2, y -= ry * d2;
lx = x - rx, ly = y - ry;
seq.back() -= d2;
} else {
x -= lx * d2, y -= ly * d2;
rx = x - lx, ry = y - ly;
seq.back() += d2;
}
d -= d2;
if (seq.back() == 0) seq.pop_back();
if (d2 == Int{0}) break;
}
return true;
}
// SBT 上の LCA
static Node lca(const Node &lhs, const Node &rhs) {
Node n;
for (int i = 0; i < min<int>(lhs.seq.size(), rhs.seq.size()); i++) {
Int val1 = lhs.seq[i], val2 = rhs.seq[i];
if ((val1 < 0) != (val2 < 0)) break;
if (val1 < 0) n.go_left(min(-val1, -val2));
if (val1 > 0) n.go_right(min(val1, val2));
if (val1 != val2) break;
}
return n;
}
friend ostream &operator<<(ostream &os, const Node &rhs) {
os << "\n";
os << "L : ( " << rhs.lx << ", " << rhs.ly << " )\n";
os << "M : ( " << rhs.x << ", " << rhs.y << " )\n";
os << "R : ( " << rhs.rx << ", " << rhs.ry << " )\n";
os << "seq : {";
for (auto &x : rhs.seq) os << x << ", ";
os << "} \n";
return os;
}
friend bool operator<(const Node &lhs, const Node &rhs) {
return lhs.x * rhs.y < rhs.x * lhs.y;
}
friend bool operator==(const Node &lhs, const Node &rhs) {
return lhs.x == rhs.x and lhs.y == rhs.y;
}
};
/**
* @brief Stern-Brocot Tree
*/
#line 9 "math/enumerate-convex.hpp"
// 下向き凸包の頂点列挙
// (xl, yl) 始点, x in [xl, xr]
// inside(x, y) : (x, y) が凸包内部か？
// candicate(x, y, c, d) : (x, y) が凸包外部にあるとする。
// 凸包内部の点 (x + sc, y + sd) が存在すればそのような s を返す
// 存在しなければ任意の値 (-1 でもよい) を返す
template <typename Int>
vector<pair<Int, Int>> enumerate_convex(
Int xl, Int yl, Int xr, const function<bool(Int, Int)>& inside,
const function<Int(Int, Int, Int, Int)>& candicate) {
assert(xl <= xr);
// inside かつ x <= xr
auto f = [&](Int x, Int y) { return x <= xr && inside(x, y); };
// (a, b) から (c, d) 方向に進めるだけ進む
auto go = [&](Int a, Int b, Int c, Int d) -> Int {
assert(f(a, b));
Int r = 1, s = 0;
while (f(a + r * c, b + r * d)) r *= 2;
while ((r /= 2) != 0) {
if (f(a + r * c, b + r * d)) s += r, a += r * c, b += r * d;
}
return s;
};
// (a, b) が out, (a + c * k, b + d * k) が in とする
// out の間進めるだけ進む
auto go2 = [&](Int a, Int b, Int c, Int d, Int k) {
assert(!inside(a, b) and inside(a + c * k, b + d * k));
Int ok = 0, ng = k;
while (ok + 1 < ng) {
Int m = (ok + ng) / 2;
(inside(a + c * m, b + d * m) ? ng : ok) = m;
}
return ok;
};
vector<pair<Int, Int>> ps;
Int x = xl, y = yl;
assert(inside(x, y) and go(x, y, 0, -1) == 0);
ps.emplace_back(x, y);
while (x < xr) {
Int a, b;
if (f(x + 1, y)) {
a = 1, b = 0;
} else {
SternBrocotTreeNode<Int> sb;
while (true) {
assert(f(x + sb.lx, y + sb.ly));
assert(!f(x + sb.rx, y + sb.ry));
if (f(x + sb.lx + sb.rx, y + sb.ly + sb.ry)) {
Int s = go(x + sb.lx, y + sb.ly, sb.rx, sb.ry);
assert(s > 0);
sb.go_right(s);
} else {
Int s = candicate(x + sb.rx, y + sb.ry, sb.lx, sb.ly);
if (s <= 0 || !inside(x + sb.lx * s + sb.rx, y + sb.ly * s + sb.ry)) {
a = sb.lx, b = sb.ly;
break;
} else {
Int t = go2(x + sb.rx, y + sb.ry, sb.lx, sb.ly, s);
sb.go_left(t);
}
}
}
}
Int s = go(x, y, a, b);
x += a * s, y += b * s;
ps.emplace_back(x, y);
}
return ps;
}
#line 2 "math/isqrt.hpp"
#line 4 "math/isqrt.hpp"
using namespace std;
// floor(sqrt(n)) を返す (ただし n が負の場合は 0 を返す)
long long isqrt(long long n) {
if (n <= 0) return 0;
long long x = sqrt(n);
while ((x + 1) * (x + 1) <= n) x++;
while (x * x > n) x--;
return x;
}
#line 2 "math/two-square.hpp"
#line 2 "internal/internal-math.hpp"
#line 2 "internal/internal-type-traits.hpp"
#line 4 "internal/internal-type-traits.hpp"
using namespace std;
namespace internal {
template <typename T>
using is_broadly_integral =
typename conditional_t<is_integral_v<T> || is_same_v<T, __int128_t> ||
is_same_v<T, __uint128_t>,
true_type, false_type>::type;
template <typename T>
using is_broadly_signed =
typename conditional_t<is_signed_v<T> || is_same_v<T, __int128_t>,
true_type, false_type>::type;
template <typename T>
using is_broadly_unsigned =
typename conditional_t<is_unsigned_v<T> || is_same_v<T, __uint128_t>,
true_type, false_type>::type;
#define ENABLE_VALUE(x) \
template <typename T> \
constexpr bool x##_v = x<T>::value;
ENABLE_VALUE(is_broadly_integral);
ENABLE_VALUE(is_broadly_signed);
ENABLE_VALUE(is_broadly_unsigned);
#undef ENABLE_VALUE
#define ENABLE_HAS_TYPE(var) \
template <class, class = void> \
struct has_##var : false_type {}; \
template <class T> \
struct has_##var<T, void_t<typename T::var>> : true_type {}; \
template <class T> \
constexpr auto has_##var##_v = has_##var<T>::value;
#define ENABLE_HAS_VAR(var) \
template <class, class = void> \
struct has_##var : false_type {}; \
template <class T> \
struct has_##var<T, void_t<decltype(T::var)>> : true_type {}; \
template <class T> \
constexpr auto has_##var##_v = has_##var<T>::value;
} // namespace internal
#line 4 "internal/internal-math.hpp"
namespace internal {
#line 10 "internal/internal-math.hpp"
using namespace std;
// a mod p
template <typename T>
T safe_mod(T a, T p) {
a %= p;
if constexpr (is_broadly_signed_v<T>) {
if (a < 0) a += p;
}
return a;
}
// 返り値：pair(g, x)
// s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g
template <typename T>
pair<T, T> inv_gcd(T a, T p) {
static_assert(is_broadly_signed_v<T>);
a = safe_mod(a, p);
if (a == 0) return {p, 0};
T b = p, x = 1, y = 0;
while (a != 0) {
T q = b / a;
swap(a, b %= a);
swap(x, y -= q * x);
}
if (y < 0) y += p / b;
return {b, y};
}
// 返り値 : a^{-1} mod p
// gcd(a, p) != 1 が必要
template <typename T>
T inv(T a, T p) {
static_assert(is_broadly_signed_v<T>);
a = safe_mod(a, p);
T b = p, x = 1, y = 0;
while (a != 0) {
T q = b / a;
swap(a, b %= a);
swap(x, y -= q * x);
}
assert(b == 1);
return y < 0 ? y + p : y;
}
// T : 底の型
// U : T*T がオーバーフローしない かつ 指数の型
template <typename T, typename U>
T modpow(T a, U n, T p) {
a = safe_mod(a, p);
T ret = 1 % p;
while (n != 0) {
if (n % 2 == 1) ret = U(ret) * a % p;
a = U(a) * a % p;
n /= 2;
}
return ret;
}
// 返り値 : pair(rem, mod)
// 解なしのときは {0, 0} を返す
template <typename T>
pair<T, T> crt(const vector<T>& r, const vector<T>& m) {
static_assert(is_broadly_signed_v<T>);
assert(r.size() == m.size());
int n = int(r.size());
T r0 = 0, m0 = 1;
for (int i = 0; i < n; i++) {
assert(1 <= m[i]);
T r1 = safe_mod(r[i], m[i]), m1 = m[i];
if (m0 < m1) swap(r0, r1), swap(m0, m1);
if (m0 % m1 == 0) {
if (r0 % m1 != r1) return {0, 0};
continue;
}
auto [g, im] = inv_gcd(m0, m1);
T u1 = m1 / g;
if ((r1 - r0) % g) return {0, 0};
T x = (r1 - r0) / g % u1 * im % u1;
r0 += x * m0;
m0 *= u1;
if (r0 < 0) r0 += m0;
}
return {r0, m0};
}
} // namespace internal
#line 2 "prime/fast-factorize.hpp"
#line 6 "prime/fast-factorize.hpp"
using namespace std;
#line 2 "misc/rng.hpp"
#line 2 "internal/internal-seed.hpp"
#line 4 "internal/internal-seed.hpp"
using namespace std;
namespace internal {
unsigned long long non_deterministic_seed() {
unsigned long long m =
chrono::duration_cast<chrono::nanoseconds>(
chrono::high_resolution_clock::now().time_since_epoch())
.count();
m ^= 9845834732710364265uLL;
m ^= m << 24, m ^= m >> 31, m ^= m << 35;
return m;
}
unsigned long long deterministic_seed() { return 88172645463325252UL; }
// 64 bit の seed 値を生成 (手元では seed 固定)
// 連続で呼び出すと同じ値が何度も返ってくるので注意
// #define RANDOMIZED_SEED するとシードがランダムになる
unsigned long long seed() {
#if defined(NyaanLocal) && !defined(RANDOMIZED_SEED)
return deterministic_seed();
#else
return non_deterministic_seed();
#endif
}
} // namespace internal
#line 4 "misc/rng.hpp"
namespace my_rand {
using i64 = long long;
using u64 = unsigned long long;
// [0, 2^64 - 1)
u64 rng() {
static u64 _x = internal::seed();
return _x ^= _x << 7, _x ^= _x >> 9;
}
// [l, r]
i64 rng(i64 l, i64 r) {
assert(l <= r);
return l + rng() % u64(r - l + 1);
}
// [l, r)
i64 randint(i64 l, i64 r) {
assert(l < r);
return l + rng() % u64(r - l);
}
// choose n numbers from [l, r) without overlapping
vector<i64> randset(i64 l, i64 r, i64 n) {
assert(l <= r && n <= r - l);
unordered_set<i64> s;
for (i64 i = n; i; --i) {
i64 m = randint(l, r + 1 - i);
if (s.find(m) != s.end()) m = r - i;
s.insert(m);
}
vector<i64> ret;
for (auto& x : s) ret.push_back(x);
sort(begin(ret), end(ret));
return ret;
}
// [0.0, 1.0)
double rnd() { return rng() * 5.42101086242752217004e-20; }
// [l, r)
double rnd(double l, double r) {
assert(l < r);
return l + rnd() * (r - l);
}
template <typename T>
void randshf(vector<T>& v) {
int n = v.size();
for (int i = 1; i < n; i++) swap(v[i], v[randint(0, i + 1)]);
}
} // namespace my_rand
using my_rand::randint;
using my_rand::randset;
using my_rand::randshf;
using my_rand::rnd;
using my_rand::rng;
#line 2 "modint/arbitrary-montgomery-modint.hpp"
#line 4 "modint/arbitrary-montgomery-modint.hpp"
using namespace std;
template <typename Int, typename UInt, typename Long, typename ULong, int id>
struct ArbitraryLazyMontgomeryModIntBase {
using mint = ArbitraryLazyMontgomeryModIntBase;
inline static UInt mod;
inline static UInt r;
inline static UInt n2;
static constexpr int bit_length = sizeof(UInt) * 8;
static UInt get_r() {
UInt ret = mod;
while (mod * ret != 1) ret *= UInt(2) - mod * ret;
return ret;
}
static void set_mod(UInt m) {
assert(m < (UInt(1u) << (bit_length - 2)));
assert((m & 1) == 1);
mod = m, n2 = -ULong(m) % m, r = get_r();
}
UInt a;
ArbitraryLazyMontgomeryModIntBase() : a(0) {}
ArbitraryLazyMontgomeryModIntBase(const Long &b)
: a(reduce(ULong(b % mod + mod) * n2)){};
static UInt reduce(const ULong &b) {
return (b + ULong(UInt(b) * UInt(-r)) * mod) >> bit_length;
}
mint &operator+=(const mint &b) {
if (Int(a += b.a - 2 * mod) < 0) a += 2 * mod;
return *this;
}
mint &operator-=(const mint &b) {
if (Int(a -= b.a) < 0) a += 2 * mod;
return *this;
}
mint &operator*=(const mint &b) {
a = reduce(ULong(a) * b.a);
return *this;
}
mint &operator/=(const mint &b) {
*this *= b.inverse();
return *this;
}
mint operator+(const mint &b) const { return mint(*this) += b; }
mint operator-(const mint &b) const { return mint(*this) -= b; }
mint operator*(const mint &b) const { return mint(*this) *= b; }
mint operator/(const mint &b) const { return mint(*this) /= b; }
bool operator==(const mint &b) const {
return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
}
bool operator!=(const mint &b) const {
return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
}
mint operator-() const { return mint(0) - mint(*this); }
mint operator+() const { return mint(*this); }
mint pow(ULong n) const {
mint ret(1), mul(*this);
while (n > 0) {
if (n & 1) ret *= mul;
mul *= mul, n >>= 1;
}
return ret;
}
friend ostream &operator<<(ostream &os, const mint &b) {
return os << b.get();
}
friend istream &operator>>(istream &is, mint &b) {
Long t;
is >> t;
b = ArbitraryLazyMontgomeryModIntBase(t);
return (is);
}
mint inverse() const {
Int x = get(), y = get_mod(), u = 1, v = 0;
while (y > 0) {
Int t = x / y;
swap(x -= t * y, y);
swap(u -= t * v, v);
}
return mint{u};
}
UInt get() const {
UInt ret = reduce(a);
return ret >= mod ? ret - mod : ret;
}
static UInt get_mod() { return mod; }
};
// id に適当な乱数を割り当てて使う
template <int id>
using ArbitraryLazyMontgomeryModInt =
ArbitraryLazyMontgomeryModIntBase<int, unsigned int, long long,
unsigned long long, id>;
template <int id>
using ArbitraryLazyMontgomeryModInt64bit =
ArbitraryLazyMontgomeryModIntBase<long long, unsigned long long, __int128_t,
__uint128_t, id>;
#line 2 "prime/miller-rabin.hpp"
#line 4 "prime/miller-rabin.hpp"
using namespace std;
#line 8 "prime/miller-rabin.hpp"
namespace fast_factorize {
template <typename T, typename U>
bool miller_rabin(const T& n, vector<T> ws) {
if (n <= 2) return n == 2;
if (n % 2 == 0) return false;
T d = n - 1;
while (d % 2 == 0) d /= 2;
U e = 1, rev = n - 1;
for (T w : ws) {
if (w % n == 0) continue;
T t = d;
U y = internal::modpow<T, U>(w, t, n);
while (t != n - 1 && y != e && y != rev) y = y * y % n, t *= 2;
if (y != rev && t % 2 == 0) return false;
}
return true;
}
bool miller_rabin_u64(unsigned long long n) {
return miller_rabin<unsigned long long, __uint128_t>(
n, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
}
template <typename mint>
bool miller_rabin(unsigned long long n, vector<unsigned long long> ws) {
if (n <= 2) return n == 2;
if (n % 2 == 0) return false;
if (mint::get_mod() != n) mint::set_mod(n);
unsigned long long d = n - 1;
while (~d & 1) d >>= 1;
mint e = 1, rev = n - 1;
for (unsigned long long w : ws) {
if (w % n == 0) continue;
unsigned long long t = d;
mint y = mint(w).pow(t);
while (t != n - 1 && y != e && y != rev) y *= y, t *= 2;
if (y != rev && t % 2 == 0) return false;
}
return true;
}
bool is_prime(unsigned long long n) {
using mint32 = ArbitraryLazyMontgomeryModInt<96229631>;
using mint64 = ArbitraryLazyMontgomeryModInt64bit<622196072>;
if (n <= 2) return n == 2;
if (n % 2 == 0) return false;
if (n < (1uLL << 30)) {
return miller_rabin<mint32>(n, {2, 7, 61});
} else if (n < (1uLL << 62)) {
return miller_rabin<mint64>(
n, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
} else {
return miller_rabin_u64(n);
}
}
} // namespace fast_factorize
using fast_factorize::is_prime;
/**
* @brief Miller-Rabin primality test
*/
#line 12 "prime/fast-factorize.hpp"
namespace fast_factorize {
using u64 = uint64_t;
template <typename mint, typename T>
T pollard_rho(T n) {
if (~n & 1) return 2;
if (is_prime(n)) return n;
if (mint::get_mod() != n) mint::set_mod(n);
mint R, one = 1;
auto f = [&](mint x) { return x * x + R; };
auto rnd_ = [&]() { return rng() % (n - 2) + 2; };
while (1) {
mint x, y, ys, q = one;
R = rnd_(), y = rnd_();
T g = 1;
constexpr int m = 128;
for (int r = 1; g == 1; r <<= 1) {
x = y;
for (int i = 0; i < r; ++i) y = f(y);
for (int k = 0; g == 1 && k < r; k += m) {
ys = y;
for (int i = 0; i < m && i < r - k; ++i) q *= x - (y = f(y));
g = gcd(q.get(), n);
}
}
if (g == n) do
g = gcd((x - (ys = f(ys))).get(), n);
while (g == 1);
if (g != n) return g;
}
exit(1);
}
using i64 = long long;
vector<i64> inner_factorize(u64 n) {
using mint32 = ArbitraryLazyMontgomeryModInt<452288976>;
using mint64 = ArbitraryLazyMontgomeryModInt64bit<401243123>;
if (n <= 1) return {};
u64 p;
if (n <= (1LL << 30)) {
p = pollard_rho<mint32, uint32_t>(n);
} else if (n <= (1LL << 62)) {
p = pollard_rho<mint64, uint64_t>(n);
} else {
exit(1);
}
if (p == n) return {i64(p)};
auto l = inner_factorize(p);
auto r = inner_factorize(n / p);
copy(begin(r), end(r), back_inserter(l));
return l;
}
vector<i64> factorize(u64 n) {
auto ret = inner_factorize(n);
sort(begin(ret), end(ret));
return ret;
}
map<i64, i64> factor_count(u64 n) {
map<i64, i64> mp;
for (auto &x : factorize(n)) mp[x]++;
return mp;
}
vector<i64> divisors(u64 n) {
if (n == 0) return {};
vector<pair<i64, i64>> v;
for (auto &p : factorize(n)) {
if (v.empty() || v.back().first != p) {
v.emplace_back(p, 1);
} else {
v.back().second++;
}
}
vector<i64> ret;
auto f = [&](auto rc, int i, i64 x) -> void {
if (i == (int)v.size()) {
ret.push_back(x);
return;
}
rc(rc, i + 1, x);
for (int j = 0; j < v[i].second; j++) rc(rc, i + 1, x *= v[i].first);
};
f(f, 0, 1);
sort(begin(ret), end(ret));
return ret;
}
} // namespace fast_factorize
using fast_factorize::divisors;
using fast_factorize::factor_count;
using fast_factorize::factorize;
/**
* @brief 高速素因数分解(Miller Rabin/Pollard's Rho)
* @docs docs/prime/fast-factorize.md
*/
#line 2 "math/gaussian-integer.hpp"
// x + yi
template <typename T>
struct Gaussian_Integer {
T x, y;
using G = Gaussian_Integer;
Gaussian_Integer(T _x = 0, T _y = 0) : x(_x), y(_y) {}
Gaussian_Integer(const pair<T, T>& p) : x(p.fi), y(p.se) {}
T norm() const { return x * x + y * y; }
G conj() const { return G{x, -y}; }
G operator+(const G& r) const { return G{x + r.x, y + r.y}; }
G operator-(const G& r) const { return G{x - r.x, y - r.y}; }
G operator*(const G& r) const {
return G{x * r.x - y * r.y, x * r.y + y * r.x};
}
G operator/(const G& r) const {
G g = G{*this} * r.conj();
T n = r.norm();
g.x += n / 2, g.y += n / 2;
return G{g.x / n - (g.x % n < 0), g.y / n - (g.y % n < 0)};
}
G operator%(const G& r) const { return G{*this} - G{*this} / r * r; }
G& operator+=(const G& r) { return *this = G{*this} + r; }
G& operator-=(const G& r) { return *this = G{*this} - r; }
G& operator*=(const G& r) { return *this = G{*this} * r; }
G& operator/=(const G& r) { return *this = G{*this} / r; }
G& operator%=(const G& r) { return *this = G{*this} % r; }
G operator-() const { return G{-x, -y}; }
G operator+() const { return G{*this}; }
bool operator==(const G& g) const { return x == g.x && y == g.y; }
bool operator!=(const G& g) const { return x != g.x || y != g.y; }
G pow(__int128_t e) const {
G res{1}, a{*this};
while (e) {
if (e & 1) res *= a;
a *= a, e >>= 1;
}
return res;
}
friend G gcd(G a, G b) {
while (b != G{0, 0}) {
trc(a, b, a / b, a % b);
swap(a %= b, b);
}
return a;
}
friend ostream& operator<<(ostream& os, const G& rhs) {
return os << rhs.x << " " << rhs.y;
}
};
#line 6 "math/two-square.hpp"
// 解が存在しない場合 (-1, -1) を返す
Gaussian_Integer<__int128_t> solve_norm_equation_prime(long long p) {
if (p % 4 == 3) return {-1, -1};
if (p == 2) return {1, 1};
long long x = 1;
while (true) {
x++;
long long z = internal::modpow<long long, __int128_t>(x, (p - 1) / 4, p);
if (__int128_t(z) * z % p == p - 1) {
x = z;
break;
}
}
long long y = 1, k = (__int128_t(x) * x + 1) / p;
while (k > 1) {
long long B = x % k, D = y % k;
if (B < 0) B += k;
if (D < 0) D += k;
if (B * 2 > k) B -= k;
if (D * 2 > k) D -= k;
long long nx = (__int128_t(x) * B + __int128_t(y) * D) / k;
long long ny = (__int128_t(x) * D - __int128_t(y) * B) / k;
x = nx, y = ny;
k = (__int128_t(x) * x + __int128_t(y) * y) / p;
}
return {x, y};
}
// p^e が long long に収まる
vector<Gaussian_Integer<__int128_t>> solve_norm_equation_prime_power(
long long p, long long e) {
using G = Gaussian_Integer<__int128_t>;
if (p % 4 == 3) {
if (e % 2 == 1) return {};
long long x = 1;
for (int i = 0; i < e / 2; i++) x *= p;
return {G{x}};
}
if (p == 2) return {G{1, 1}.pow(e)};
G pi = solve_norm_equation_prime(p);
vector<G> pows(e + 1);
pows[0] = 1;
for (int i = 1; i <= e; i++) pows[i] = pows[i - 1] * pi;
vector<G> res(e + 1);
for (int i = 0; i <= e; i++) res[i] = pows[i] * (pows[e - i].conj());
return res;
}
// 0 <= arg < 90 の範囲の解のみ返す
vector<Gaussian_Integer<__int128_t>> solve_norm_equation(long long N) {
using G = Gaussian_Integer<__int128_t>;
if (N < 0) return {};
if (N == 0) return {G{0}};
auto pes = factor_count(N);
for (auto& [p, e] : pes) {
if (p % 4 == 3 && e % 2 == 1) return {};
}
vector<G> res{G{1}};
for (auto& [p, e] : pes) {
vector<G> cur = solve_norm_equation_prime_power(p, e);
vector<G> nxt;
for (auto& g1 : res) {
for (auto& g2 : cur) nxt.push_back(g1 * g2);
}
res = nxt;
}
for (auto& g : res) {
while (g.x <= 0 || g.y < 0) g = G{-g.y, g.x};
}
return res;
}
// x,y 両方非負のみ, 辞書順で返す
vector<pair<long long, long long>> two_square(long long N) {
if (N < 0) return {};
if (N == 0) return {{0, 0}};
vector<pair<long long, long long>> ans;
for (auto& g : solve_norm_equation(N)) {
ans.emplace_back(g.x, g.y);
if (g.y == 0) ans.emplace_back(g.y, g.x);
}
sort(begin(ans), end(ans));
return ans;
}
#line 9 "verify/verify-unit-test/enumerate-convex.test.cpp"
//
using namespace Nyaan;
vector<pair<long long, long long>> calc(ll N) {
ll m = isqrt(N);
// (0, m) を中心とする半径 sqrt(N) の円
auto inside = [&](ll x, ll y) {
return y >= m or x * x + (y - m) * (y - m) <= N;
};
auto candicate = [&](ll x, ll y, ll c, ll d) {
// (x + sc)^2 + (y - m + sd)^2 <= N
ll A = c * c + d * d;
ll B = 2 * c * x + 2 * d * (y - m);
// A s^2 + B s + const <= 0
ll num = -B, den = 2 * A;
ll quo = num / den, rem = num % den;
if (rem < 0) quo--, rem += den;
if (2 * rem > den) quo++, rem -= den;
return quo;
};
auto ans = enumerate_convex<ll>(0, 0, m, inside, candicate);
vector<pair<long long, long long>> res;
each2(x, y, ans) if (x * x + (y - m) * (y - m) == N) {
res.emplace_back(x, m - y);
}
sort(begin(res), end(res));
return res;
}
void check(long long N) {
auto ac = two_square(N);
auto ad = calc(N);
assert(ac == ad);
}
void q() {
rep1(N, 1000) check(N);
rep(t, 100) check(rng(1001, TEN(9)));
check(TEN(18));
trc2("OK");
inl(a, b);
out(a + b);
}
void Nyaan::solve() {
int t = 1;
// in(t);
while (t--) q();
}
verify/verify-unit-test/enumerate-convex.test.cpp