/* Factoring by determining how much over the mark a guess
 * at phi(n) is. This method is only valid for an RSA composite number,
 * if gives junk answers for others.
 */
#include <stdio.h>
#include "isqrt.c"

/* Modular arithmetic routines */
#include "xpymn.c"
#include "xbyymn.c"
#include "xtoymn.c"

main () { long n,s,p,q,test,guess,x,y,oversht;

  /* get number to work on */
  printf ("Product of two odd primes: N=\n");
  scanf("%ld",&n); printf ("N=%ld\n",n);
  if (n<0) {
     printf ("N must be greater than zero\n"); return 1;
  } else if (0 == (n & 1)) {
     printf ("N must be odd, this number has a factor of 2\n"); return 2;
  }

  /* Form a good guess for phi (which could have been supplied) */
  s = isqrt(n); guess = n-s-s+1; /* known to not undershoot */

  /* Compute the overshoot of the guess for phi. */
  test = xtoymn((long)2,guess,n); oversht = 0;
  while (test != 1) {
    if ( test & 1 ) { test += n; }
    test >>= 1; oversht++;
  }
  /* Actually, if test is not an explicit power of two, you know
   * you have overshot by at least log2(n), but we'll ignore that speedup here.
   */

  /* derive actual phi */
  guess -= oversht;

  /* and try to factor by it. */
  y = (n-guess+1)/2; x = isqrt(y*y-n);
  p = y-x; q = y+x;
  if (p*q==n) {
     printf ("Factors: %ld and %ld\n",p,q);
  } else {
     printf ("N is not the product of TWO prime numbers?\n");
     printf (" (Or has very very very small order)\n");
     printf (" but we know lambda(N) (The actual order)\n");
     printf (" is a multiple of %ld.  *IF* N is the product of two\n", guess);
     printf (" primes, then phi(N) is also a multiple\n");
  }

  return p;
}