#include <cstdlib>
#include <string>
#include <iostream>
#include <boost/cstdint.hpp>
#include <boost/integer.hpp>
#include <gmp.h>

using namespace std;
using namespace boost;

typedef struct sol {
    mpz_t p;
    mpz_t q;
} sol_t;

// Left classifier
void lc(mpz_t * l, const mpz_t * x1, const mpz_t * x2) {
    mpz_t sub1, div1;
    /* Apply l = (((x2-x1)*0.25) + x1 )^2 */

    /* Substract low bound from high bound */
    mpz_sub(sub1, *x2, *x1);
    /* Divide by 4 */
    mpz_tdiv_q_2exp(div1, sub1, 2);
    /* Add low bound */
    mpz_add(*l, div1, *x1);
    /* Powerize to 2 */
    mpz_pow_ui(*l, *l, 2);
}

void rc(mpz_t * r, const mpz_t * x1, const mpz_t * x2) {
    mpz_t sub1, tmp1, tmp2, sum1, div1;
    /* Apply: ((0.75 * (x2-x1)) + x1)^2 +
              ((2*(0.5  * (x2-x1)) + x1)^2) - lc(x1, x2) -
              ((0.75 * (x2-x1)) + x1)^2 */
    /* Substract low bound from high bound */
    mpz_sub(sub1, *x2, *x1);
    /* Divide by 2 */
    mpz_tdiv_q_2exp(div1, sub1, 1);
    /* Add low bound */
    mpz_add(sum1, div1, *x1);
    /* Powerize to 2 */
    mpz_pow_ui(tmp1, sum1, 2);
    /* Multiply by 2 */
    mpz_mul_2exp(tmp1, tmp1, 1);
    /* Left classifier */
    lc(&tmp2, x1, x2);
    /* Add low bound */
    mpz_sub(*r, tmp1, tmp2);
}

void factImpl(const mpz_t * k, mpz_t * l, mpz_t * h) {
    // Continue until granularity is not unity
    mpz_t sub1, div1, tmp1, tmp2;
    while ( mpz_cmp(*h, *l) >= 0 ) {
        /* h - l */
        mpz_sub(sub1, *h, *l);
        /* Divide by 2 */
        mpz_tdiv_q_2exp(div1, sub1, 1);
        /* euclidian left part */
        lc(&tmp1, l, h);
        mpz_sub(tmp1, *k, tmp1);
        mpz_abs(tmp1, tmp1);
        rc(&tmp2, l, h);
        mpz_sub(tmp2, *k, tmp2);
        mpz_abs(tmp2, tmp2);
        if ( mpz_cmp(tmp2, tmp1) > 0 ) {
            mpz_sub(*h, *h, div1);
            if ( mpz_cmp(*l, *h) > 0 ) {
                mpz_add_ui(*h, *l, 1);
            }
        } else {
            mpz_add_ui(*l, *l, 1);
        }
    }
}

void follow_impl(const mpz_t * k,
                 mpz_t * x, mpz_t * y,
                 const int dx, const int dy,
                 const mpz_t * w, const mpz_t * h) {
    mpz_t c, lh, hh, hl, ll, xi, yi, xy, xiy, xyi, xiyi, iix, iiy;
    /* Initial warmup */
    lc(&c, x, y);
    mpz_set_si(iix, dx);
    mpz_set_si(iiy, dy);
    mpz_add(xi, *x, iix);
    mpz_add(yi, *y, iiy);
    mpz_mul(xy, *x, *y);
    mpz_mul(xiy, xi, *y);
    mpz_mul(xyi, *x, yi);
    mpz_mul(xiyi, xi, yi);

    do {
        /* compute ll */
        mpz_sub(ll, xy, *k);
        mpz_abs(ll, ll);
        /* compute hl */
        mpz_sub(hl, xyi, c);
        mpz_abs(hl, hl);
        /* compute lh */
        mpz_sub(lh, xiy, c);
        mpz_abs(lh, lh);
        /* compute hh */
        mpz_sub(hh, xiyi, c);
        mpz_abs(hh, hh);

        /* Check if found */
        if ( mpz_sgn(ll) == 0 ) {
            break;
        } else if ( mpz_cmp(hl, hh) > 0 && mpz_cmp(lh, hh) > 0) {
            mpz_add(*x, *x, iix);
            mpz_add(*y, *y, iiy);
        } else if ( mpz_cmp(hl, lh) > 0 ) {
            mpz_add(*x, *x, iix);
        } else {
            mpz_add(*y, *y, iiy);
        }
        mpz_add(xi, *x, iix);
        mpz_add(yi, *y, iiy);
        mpz_mul(xy, *x, *y);
        mpz_mul(xiy, xi, *y);
        mpz_mul(xyi, *x, yi);
        mpz_mul(xiyi, xi, yi);
    } while( mpz_cmp(*w, xi) > 0 && mpz_cmp_ui(xi, 1) > 0 &&
            mpz_cmp(*h, yi) > 0 && mpz_cmp_ui(yi, 1) > 0);
}
void
resolve(const mpz_t * product, sol *prims) {
    mpz_t hl, hh, x, y, pq;   // Hyperbol bounds (low bound and high bound)

    mpz_set_ui(prims->p, 0);
    mpz_set_ui(prims->q, 0);

    mpz_set(hh, *product);
    factImpl(product, &hl, &hh);
    mpz_set(x, hl);
    mpz_set(y, hh);
    std::cerr << "Analysis hyperbol direction: up-right" << std::endl;
    follow_impl( product, &x, &y, 1, -1, product, product );
    mpz_mul(pq, x, y);
    mpz_set(x, hl);
    mpz_set(y, hh);
    if ( mpz_cmp(pq, *product) != 0) {
        std::cerr << "Analysis hyperbol direction: down-left" << std::endl;
        follow_impl( product, &x, &y, -1, 1, product, product );
        mpz_mul(pq, x, y);
    }
    if ( mpz_cmp(pq, *product) != 0) {
        mpz_set(prims->p, x);
        mpz_set(prims->p, y);
    }
}

int main(int argc, char *argv[]) {
    typedef mpz_t bigint_t;
    if(2 == argc) {
        bigint_t pq;
        sol_t res;
        mpz_set_ui(res.p, 0);
        mpz_set_ui(res.q, 0);
        mpz_init_set_str(pq, argv[1], 10);
        resolve(&pq, &res);
        printf("solution: pq = (p*q) =\n");
        mpz_out_str(stdout, 10, res.p);
        printf("\n*\n");
        mpz_out_str(stdout, 10, res.q);
    } else {
        printf("Usage: %s <p*q>\n", argv[0]);
        return -1;
    }
    return 0;
}