2020.7.6 -- Miller_RabinとPollard_Rhoアルゴリズム

#include "bits/stdc++.h"
using namespace std;
#define ll long long
#define ull unsigned ll
#define pb push_back
const int RhoLimit = 10000; //    ,    
const int Rhoc = 12;        //         
vector<ll> Fac;

//gcd
ll gcd(ll x, ll y)
{
    if (y == 0)
        return x;
    return gcd(y, x % y);
}
//mod   
ll qmul(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 qpow(ll a, ll n, ll mod)
{
    ll res = 1;
    while (n)
    {
        if (n & 1)
            res = qmul(res, a, mod);
        a = qmul(a, a, mod);
        n >>= 1;
    }
    return res;
}
//  Pollard-Rho              
bool Miller_Rabin(ll n)
{
    if (n == 2)
        return 1;
    if (n < 2 || !(n & 1))
        return 0;
    ll m = n - 1, k = 0; //n-1 = 2^k + m
    while (!(m & 1))
    {
        k++;
        m >>= 1;
    }
    int prime[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 53};
    for (int i = 0; i < 14; i++)
    {
        ll a = prime[i] % (n - 1) + 1;
        ll x = qpow(a, m, n);
        ll y;
        for (int j = 1; j <= k; j++)
        {
            y = qmul(x, x, n);
            if (y == 1 && x != 1 && x != n - 1)
                return false;
            x = y;
        }
        if (y != 1)
            return false;
    }
    return true;
}
//     
ll rand_(ll x)
{
    ull p = rand() * rand();
    p *= rand() * rand();
    p += rand();
    return p % x;
}
//
ll Pollard_Rho(ll n)
{
    ll x = rand_(n - 2) + 2;
    ll y = x;
    int step = 2;
    for (int i = 1;; i++)
    {
        x = qmul(x, x, n) + Rhoc;
        if (x >= n)
            x -= n;
        if (x == y)
            return -1;
        ll d = gcd(abs(x - y), n);
        if (d > 1)
            return d;
        if (i == step)
        {
            i = 0;
            y = x;
            step <<= 1;
        }
    }
}

void solve(ll n)
{
    if (Miller_Rabin(n))
    {
        Fac.pb(n);
        return;
    }
    ll p;
    for (int i = 0; i != RhoLimit; ++i)
    {
        p = Pollard_Rho(n);
        if (p != -1)
            break;
    }
    if (p == -1)
        return;
    solve(p);
    solve(n / p);
}

int main()
{
    ll t;
    while (cin >> t)
    {
        solve(t);
        sort(Fac.begin(), Fac.end());
        for (auto i : Fac)
            cout << i << " ";
        cout << endl;
        Fac.clear();
    }
    return 0;
}