#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cassert>
#include <cstring>
#include <cinttypes>
#include <iostream>
#include <map>
#include <unordered_map>
#include <algorithm>
#include <bitset>
#include <random>
using namespace std;
#define forn(i, n) for (int i = 0; i < (int)(n); i++)
typedef long double dbl;
typedef __int128 ll;
ll Add(ll a, ll b, ll mod) {
return (a += b) >= mod ? a-mod : a;
}
ll Mul(ll a, ll b, ll mod) {
if (!b)
return 0;
ll a2 = Add(a, a, mod);
a2 = Mul(a2, b >> 1, mod);
return (b & 1) ? Add(a2, a, mod) : a2;
}
int BinPow(int a, int n, int mod) {
if (n == 0)
return 1;
int r = BinPow(a * a % mod, n / 2, mod);
return n % 2 == 0 ? r : (r * a % mod);
}
const int PN = 4500; // number of primes to use in Gauss
int pn = 0, p[PN], v[PN];
const int maxVN = PN+30;
int vn;
ll a[maxVN];
char m0[maxVN][PN];
bitset<PN> m[maxVN];
bitset<maxVN> linear[maxVN];
int ind[maxVN];
bool IsPrime(int x) {
for (int d = 0; d < pn; d++)
if (x % p[d] == 0)
return 0;
return 1;
}
void InitP() {
for (int x = 2; pn < PN; x++)
if (IsPrime(x))
p[pn++] = x;
}
bool NewVector(ll _a) {
a[vn] = _a;
forn(i, PN)
m[vn][i] = (m0[vn][i] = v[i]) & 1;
linear[vn].reset();
linear[vn][vn] = 1;
forn(i, vn)
if (m[vn][ind[i]])
m[vn] ^= m[i], linear[vn] ^= linear[i];
int i = 0;
while (i < PN && !m[vn][i])
i++;
if (i == PN)
return 0;
assert(vn < maxVN);
ind[vn++] = i;
return 1;
}
string str(ll n) {
char s[64];
int k = 0;
while (n)
s[k++] = (int)(n % 10) + '0', n /= 10;
if (!k && !n)
s[k++] = '0';
s[k] = 0;
reverse(s, s + k);
return string(s);
}
// min x : x*x > n
ll sqroot(ll n) {
ll l = 1, r = 1;
while (r * r <= n)
r *= 2;
while (r - l > 1) {
ll m = (l + r) / 2;
(m * m <= n ? l : r) = m;
}
return r;
}
ll Solve(ll n) {
ll x = sqroot(n);
vector<int> root(PN, -1);
forn(j, PN) if (p[j] >= 3)
if (BinPow(n % p[j], (p[j]-1)/2, p[j]) == 1) {
int x0 = (x * x - n) % p[j];
int dx = (2 * x) % p[j];
for (int i = 0; i < p[j]; i++)
if ((x0 + i * (i + dx)) % p[j] == 0) {
root[j] = i;
break;
}
}
const int M = 1e4; // step
vector<ll> rest(M);
vn = 0;
for (int start = 0;; start += M) {
forn(k, M) {
auto i = k + start;
rest[k] = (x + i) * (x + i) - n;
while ((rest[k] & 1) == 0)
rest[k] >>= 1;
}
forn(j, PN) if (p[j] >= 3 && root[j] != -1) {
auto UseRoot = [&](ll i) {
i %= p[j];
if (i < 0) i += p[j];
i += start - start % p[j];
if (i < start) i += p[j];
for (int pos = i; pos < start + M; pos += p[j]) {
int k = pos - start;
rest[k] /= p[j];
while (rest[k] % p[j] == 0)
rest[k] /= p[j];
}
};
UseRoot(root[j]);
UseRoot(-2 * x - root[j]);
}
forn(pos, M) if (rest[pos] == 1) {
memset(v, 0, sizeof(v));
ll a = x + pos + start;
ll x = a * a - n;
forn(j, pn)
while (x % p[j] == 0)
x /= p[j], v[j]++;
if (NewVector(a))
continue;
ll A = 1, B = 1;
memset(v, 0, sizeof(v));
forn(i, vn + 1)
if (linear[vn][i]) {
A = Mul(A, ::a[i], n);
forn(j, pn)
v[j] += m0[i][j];
}
forn(j, PN) {
v[j] /= 2;
while (v[j]--)
B = Mul(B, p[j], n);
}
A %= n, B %= n;
x = gcd(n, A + B);
auto y = gcd(n, max(A, B) - min(A, B));
auto good = [&](int64_t x) { return 1 < x && x < n && n % x == 0; };
fprintf(stderr, "after %d look ups there are %d vectors [pn=%d] : there is candidate\n", pos+start+1, vn, PN);
if (good(x)) return x;
if (good(y)) return y;
}
}
}
int main() {
InitP();
string s;
cin >> s;
ll n = 0;
for (char c : s)
n = n * 10 + c - '0';
cout << str(Solve(n)) << endl;
}
|