#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "gmp.h"
#include "ecm.h"

#define USE_ECM /* use ECM when no more factors are given */

/* usage: aliq < format.compact */

unsigned int number_of_digits(n) mpz_t n;
{
  return ((unsigned int) (1.0+log(mpz_get_d(n))/log(10.0)));
}

void
find_prime_factor_by_ecm (mpz_t p, mpz_t n)
{
  double B1;
  char *s;
  int res;

  s = mpz_get_str (NULL, 10, n);
  printf ("[c%u = ", strlen (s));
  mpz_out_str (stdout, 10, n);
  printf ("]");
  fflush (stdout);
  free (s);

  if (mpz_divisible_ui_p (n, 2))
    mpz_set_ui (p, 2);
  else if (mpz_divisible_ui_p (n, 3))
    mpz_set_ui (p, 3);
  else if (mpz_divisible_ui_p (n, 5))
    mpz_set_ui (p, 5);
  else if (mpz_divisible_ui_p (n, 7))
    mpz_set_ui (p, 7);
  else if (mpz_perfect_power_p (n))
    {
      unsigned long k;
      
      for (k = 2; ; k++)
	if (mpz_root (p, n, k))
	  return;
    }
  else
    {
      /* n is prime to 2 and 3 */
      B1 = 1.0;
      while ((res = ecm_factor (p, n, B1, NULL)) == 0 || mpz_probab_prime_p (p, 25) == 0)
	{
	  if (res == 0)
	    B1 += sqrt (B1);
#ifdef DEBUG
	  printf ("n="); mpz_out_str (stdout, 10, n);
	  printf (" B1=%f\n", B1);
#endif
	}
    }

  s = mpz_get_str (NULL, 10, p);
  printf ("[p%u]", strlen (s));
  free (s);
  fflush (stdout);
}

int
main (argc,argv) int argc; char *argv[];
{
  unsigned long int n0,i,j,k,ind,maxind;
  mpz_t n,p,s,m,r,maxn;
  char c;
#ifndef NOBOUND
  unsigned long int B,*pr; double b; char *prime;
#endif

  mpz_init (r);
  mpz_inp_str (r, stdin, 10);
  n0 = mpz_get_ui (r);
  mpz_init_set_ui(n,n0);
#ifndef NOBOUND
  mpz_inp_str (r, stdin, 10);
  B=mpz_get_ui(r);
  /* construct table of small primes:
     p(k) >= k*(ln(p(k))-2) for k>=5 [Dusart] ==> the number of primes
      less or equal to B is <= B/(ln(B)-2) */
  if (B<=7) b=4.1;
  else { b=(double)B; b=b/(log(b)-2.0); }
  pr=(unsigned long int*)malloc((int)b*sizeof(int));
  prime = (char*) malloc(B);
  for (i=2;i<B;i++) prime[i]='1';
  j=2; do {
      for (i=j*j;i<B;i+=j) prime[i]='0';
      while (prime[++j]=='0');
    } while (j*j<B);
  for (k=0,i=2;i<B;i++)
    if (prime[i]=='1') pr[k++]=i;
#endif
  mpz_init(s); ind=0; mpz_init(p); mpz_init(m);
  mpz_init_set(maxn,n); maxind=0;
  while (mpz_cmp_ui (n, 1) != 0
#ifndef USE_ECM
	 && !feof (stdin)
#endif
	 )
    {
      printf("%lu:%lu = ",n0,ind);
      mpz_out_str(stdout,10,n);
      printf(" = ");
      ind++;
      mpz_set_ui(s,1); i=0; mpz_set(m,n);
#ifndef NOBOUND
    /* divide by small primes */
    while (mpz_cmp_ui(n,1)!=0 && mpz_probab_prime_p(n,25)==0 && i<k) {
      j=0; while (mpz_fdiv_ui(n,pr[i])==0) {
	j++;
	mpz_fdiv_q_ui(n,n,pr[i]);
      }
      if (j>0) {
	printf("%d",pr[i]);
	if (j>1) printf("^%lu",j);
	if (mpz_cmp_ui(n,1)!=0)
	  putchar('*');
	mpz_set_ui(p,pr[i]);
	mpz_pow_ui(p,p,j+1);
	mpz_sub_ui(p,p,1);
	mpz_fdiv_q_ui(p,p,pr[i]-1);
	mpz_mul(s,s,p);
      }
      i++;
    }
#endif
    /* divide by given primes */
    while (mpz_cmp_ui (n, 1) != 0 && mpz_probab_prime_p (n, 25) == 0) {
#ifdef USE_ECM
      if (feof (stdin) || mpz_inp_str (p, stdin, 10) == 0)
	find_prime_factor_by_ecm (p, n);
#else
      mpz_inp_str (p, stdin, 10);
#endif
      j=0; while (1) {
	mpz_mod(r,n,p);
	if (mpz_cmp_ui(r,0)==0) { j++; mpz_divexact(n,n,p); }
	else break;
      }
      if (j>0) {
	mpz_out_str(stdout,10,p);
	if (j>1) printf("^%lu",j);
	if (mpz_cmp_ui(n,1)!=0)
	  putchar('*');
	fflush (stdout);
	mpz_pow_ui(r,p,j+1);
	mpz_sub_ui(r,r,1);
	mpz_sub_ui(p,p,1);
	mpz_divexact(r,r,p);
	mpz_mul(s,s,r);
      }
      else if (!feof(stdin)) {
	printf("error: given prime p="); mpz_out_str(stdout,10,p);
	printf(" does not divide n="); mpz_out_str(stdout,10,n);
	putchar('\n');
	exit(1);
      }
      else {
	printf("?\n");
	fflush (stdout);
	break;
      }
    }
#ifndef USE_ECM
    if (!feof (stdin))
#endif
      {
	if (mpz_cmp_ui (n, 1) !=0 && mpz_probab_prime_p (n, 25) == 0)
	  {
	    printf ("error: could not find factor of ");
	    mpz_out_str (stdout, 10, n);
	    putchar ('\n');
	    exit (1);
	  }
	/* rest is prime */
	if (mpz_cmp_ui (n,1) != 0)
	  {
	    mpz_out_str (stdout, 10, n); 
	    mpz_add_ui (n, n, 1);
	    mpz_mul (s, s, n);
	  }
	putchar ('\n');
	mpz_sub (n, s, m);
	if (mpz_cmp (n, maxn) > 0)
	  {
	    mpz_set (maxn, n);
	    maxind = ind;
	  }
      }
    /* parse optional comment # ... $ */
    while (isblank(c = getchar ()));
    if (c == '#')
      while ((c = getchar ()) != '\n');
    }
  printf("largest term has %u digits at index %d\n",
	 number_of_digits(maxn),maxind);
}