gmpfactor.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gmp.h>
#include <signal.h>
#include <time.h>

//f(x)=13231/x
//f'(x)=-13231/(x^2)
//g'(x)=f'(x)+1=1-13231/(x^2)
//g(x)=x+13231/x
static void g(mpf_t res, mpf_t x, mpf_t n) {
  // res = n/x + x
  mpf_div(res, n, x);
  mpf_add(res, res, x);
}

// g(res1|res2)=x
// thanks for hints on how to solve this, wolfram alpha!
static mpf_t inv_g_halfx, inv_g_root;
static void inv_g(mpf_t res1, mpf_t res2, mpf_t x, mpf_t n) {
  // res^2 - res*x = -n
  // res^2 - 2*res*(x/2) + (x/2)^2 = (x/2)^2 - n
  // (res-(x/2))^2 = (x/2)^2 - n
  // res-(x/2) = +/- sqrt((x/2)^2 - n)
  // res = (x/2) +/- sqrt((x/2)^2 - n)
  
  // halfx = x/2
  mpf_div_ui(inv_g_halfx, x, 2);
  
  // root = sqrt(halfx*halfx - n)
  mpf_mul(inv_g_root, inv_g_halfx, inv_g_halfx);
  mpf_sub(inv_g_root, inv_g_root, n);
  mpf_sqrt(inv_g_root, inv_g_root);
  
  // res1 = halfx-root
  mpf_sub(res1, inv_g_halfx, inv_g_root);
  
  // res2 = halfx+root
  mpf_add(res2, inv_g_halfx, inv_g_root);
}

static int stop = 0;

void int_handler() {
  if (stop == 1) {
    exit(0);
  }
  stop = 1;
}

int main(int argc, char *argv[]) {
  // I only want to work with 1024-bit stuff and smaller, so
  // 4096 bits should be sufficient for the floats
  mpf_set_default_prec(4096);
  
  // Initialize temporary static variables for other functions - yes, this is
  // totally non-reentrant, but would you want to pass them as parameters?
  mpf_init(inv_g_halfx);
  mpf_init(inv_g_root);
  
  // Initialize our own temporary variables
  mpf_t max_x, min_g, cur_g, x1, x2, n, lastx1, stepsize;
  mpf_inits(max_x, min_g, cur_g, x1, x2, n, lastx1, stepsize, NULL);
  mpz_t x1_rounded, x2_rounded, x1_x2_product, stepsize_int, n_root_int;
  mpz_inits(x1_rounded, x2_rounded, x1_x2_product, stepsize_int, n_root_int, NULL);

  if (argc != 2) { puts("please invoke with the product of two primes"); return 1; }
  
  mpz_t int_n;
  if (mpz_init_set_str(int_n, argv[1], 10)) {
    puts("libgmp says that you didn't supply a valid base10 number; exiting\n");
    return 1;
  }
  mpf_set_z(n, int_n);
  
  printf("rounded square root of the product: ");
  mpz_sqrt(n_root_int, int_n);
  mpz_out_str(stdout, 10, n_root_int);
  printf("\n");
  
  // max_x = floor(sqrt(n))
  mpf_sqrt(max_x, n);
  mpf_floor(max_x, max_x);

  // min_g = ceil(g(max_x, n))
  g(min_g, max_x, n);
  mpf_ceil(min_g, min_g);
  
  signal(SIGINT, int_handler);
  
  for (mpf_set(cur_g, min_g); 1; mpf_add_ui(cur_g, cur_g, 1)) {
    mpf_set(lastx1, x1);

    inv_g(x1, x2, cur_g, n);

    printf("running, stepsize=");
    if (mpf_cmp_ui(lastx1, 0) == 0) {
      printf("0");
    } else {
      mpf_sub(stepsize, lastx1, x1);
      mpz_set_f(stepsize_int, stepsize);
      mpz_out_str(stdout, 10, stepsize_int);
    }
    printf(" p1=");
    mpz_out_str(stdout, 10, x1_rounded);
    printf(" p2=");
    mpz_out_str(stdout, 10, x2_rounded);
    printf("\n");

    // check whether res1%1=0. because floating-point calculations aren't precise, we have to
    // round res1 and res2, cast them to integers, multiply them with each other and
    // check whether the result is n.
    mpz_set_f(x1_rounded, x1);
    mpz_set_f(x2_rounded, x2);
    mpz_mul(x1_x2_product, x1_rounded, x2_rounded);
    if (mpz_cmp(x1_x2_product, int_n) == 0 || stop) {
      if (stop) {
        printf("stopped, p1=");
      } else {
        printf("solved, p1=");
      }
      mpz_out_str(stdout, 10, x1_rounded);
      printf(" p2=");
      mpz_out_str(stdout, 10, x2_rounded);
      printf("\n");
      return 0;
    }
    if (mpf_cmp_ui(x1, 1) < 0) {
      printf("p1 must be <1, but that's impossible. terminating.\n");
      return 1;
    }
  }
}