NJUST 1722 所有的平方差 Pollard_Rho大数分解

原题链接：所有的平方差

代码如下：

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <set>

using namespace std;

typedef long long int ll;
const int MAX_SIZE = 3000;
const ll mod = 10000000019;
ll fac[MAX_SIZE], ct;
ll N;

ll GCD(ll x, ll y)
{
    return !y ? x : GCD(y, x % y);
}

ll mod_mul(ll a, ll b, ll mod)
{
    ll res = 0;
    while (b)
    {
        if (b & 1)
            res = (res + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return res;
}

ll mod_pow(ll a, ll n, ll mod)
{
    ll res = 1;
    while (n)
    {
        if (n & 1)
            res = mod_mul(res, a, mod);
        a = mod_mul(a, a, mod);
        n >>= 1;
    }
    return res;
}

bool Miller_Rabin(ll n)
{
    if (n == 2)
        return true;
    if (n < 2 || !(n & 1))
        return false;
    ll m = n - 1, k = 0;
    while (!(m & 1))
    {
        k++;
        m >>= 1;
    }
    for (int i = 1; i <= 20; i++)
    {
        ll a = rand() % (n - 1) + 1;
        ll x = mod_pow(a, m, n);
        ll y;
        for (int j = 1; j <= k; j++)
        {
            y = mod_mul(x, x, n);
            if (y == 1 && x != 1 && x != n - 1)
                return false;
            x = y;
        }
        if (y != 1)
            return false;
    }
    return true;
}

ll Pollard_Rho(ll n, ll c)
{
    ll i = 1, k = 2;
    ll x = rand() % (n - 1) + 1, y = x;
    while (1)
    {
        i++;
        x = (mod_mul(x, x, n) + c) % n;
        ll d = GCD((y - x + n) % n, n);
        if (d != 1 && d != n)
            return d;
        if (y == x)
            return n;
        if (i == k)
        {
            y = x;
            k <<= 1;
        }
    }
}

void find(ll n, ll c)
{
    if (n == 1)
        return;
    if (Miller_Rabin(n))
    {
        fac[ct++] = n;
        return;
    }
    ll p = n;
    ll k = c;
    while (p >= n)
        p = Pollard_Rho(p, c--);
    find(p, k);
    find(n / p, k);
}

ll ans, s;
set<ll> S;
void dfs(ll cur, ll val)
{
    if (cur == ct)
    {
        ll q = val, p = N / val;
        if (((p + q) & 1) == 0 && !S.count(q))
        {
            S.insert(q);
            ll a = (p + q) / 2, b = (p - q) / 2;
            a = mod_mul(a, a, mod - 1);
            b = mod_mul(b, b, mod - 1);
            ans += (a + b);
            ans %= (mod - 1);
        }
        return;
    }
    if (val * fac[cur] <= s)
        dfs(cur + 1, val * fac[cur]);
    dfs(cur + 1, val);
}

int main()
{
    //freopen("test.txt", "r", stdin);

    int T, Case = 1;
    scanf("%lld", &T);
    while (T--)
    {
        scanf("%lld", &N);
        printf("Case %d: ", Case++);
        if (N == 1)
        {
            printf("1\n");
            continue;
        }
        ct = 0;
        find(N, 113);
        sort(fac, fac + ct);
        ans = 0;
        s = (ll)sqrt(N * 1.0);
        S.clear();
        dfs(0, 1);
        if (S.size() > 0)
            printf("%lld\n", mod_pow(N, ans, mod));
        else
            printf("-1\n");
    }
    return 0;
}