View difference between Paste ID: PKCcfRLt and p1NAeAiL

 #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;
 }