/* Parallel Factoring program written in C, using MPI

   By Doug Hoyte, April 2003
 */

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <malloc.h>
#include <mpi.h>
#include <mp.h>


MINT *factorloop(MINT *lrg, MINT *bot, MINT *top) {

  MINT *ind1, *ind2, *step, *mod, *quo, *zero, *tp;

  ind1 = itom(0);
  ind2 = itom(0);

  move(bot, ind1);

  step = itom(2);
  mod = itom(0);
  quo = itom(0);
  zero = itom(0);

  tp = ind1;

  for(;;) {
    mdiv(lrg, tp, quo, mod);

    if (mcmp(mod, zero) == 0) break;

    if (tp == ind1) {
      madd(tp, step, ind2);
      tp = ind2;
      if (mcmp(tp, top) > 0) return NULL;
    } else {
      madd(tp, step, ind1);
      tp = ind1;
      if (mcmp(tp, top) > 0) return NULL;
    }

  }

  mfree(step);
  mfree(mod);
  mfree(quo);
  mfree(zero);

  return tp;

}


MINT *gettopbound(MINT *tp) {

  MINT *t2, *junk, *one, *two, *mod;

  t2 = itom(0);
  junk = itom(0);
  one = itom(1);
  two = itom(2);
  mod = itom(0);

  msqrt(tp, t2, junk);

  madd(t2, one, t2);

  mdiv(t2, two, junk, mod);

  if (mcmp(mod, one) != 0) madd(t2, one, t2);

  mfree(one);
  mfree(two);
  mfree(mod);
  mfree(junk);

  return t2;

}



MINT *atomint(char *num) {

  MINT *tp, *cval, *radix;

  tp = itom(0);
  radix = itom(10);

  if (isdigit(*num)) {
    while(1) {

      cval = itom(*num - '0');
      madd(tp, cval, tp);
      mfree(cval);

      num++;

      if (!isdigit(*num)) break;

      mult(tp, radix, tp);

    }
  }

  mfree(radix);

  return tp;

}



void minttoa(MINT *tp, char *buf, int size) {

  MINT *count, *radix, *mod, *zero;
  int len=0, tpswap, iup, idown;
  char *tpstr;

  count = itom(0);
  zero = itom(0);
  mod = itom(0);
  radix = itom(10);

  memset(buf, '\0', size);

  /* Cheap hack to copy MINT: */
  rpow(tp, 1, count);

  while(1) {

    mdiv(count, radix, count, mod);

    tpstr = mtox(mod);
    buf[len++] = *tpstr;
    free(tpstr);

    if (mcmp(zero, count) >= 0) break;

    if (len>=size) return;
  }


  /* Reverse the number: */

  iup = strlen(buf) - 1;
  idown = 0;
  for(; idown < iup; idown++, iup--) {
    tpswap = buf[idown];
    buf[idown] = buf[iup];
    buf[iup] = tpswap;
  }


}



MINT *dofactor(char *num, int rank, int numprocs) {

  MINT *mnum, *mrank, *mnumprocs;
  MINT *root, *part, *top, *bot;
  MINT *junk, *res;
  MINT *one, *two, *mod;

  mnum = atomint(num);
  mrank = itom(rank); 
  mnumprocs = itom(numprocs); 

  root = gettopbound(mnum);

  part = itom(0);
  junk = itom(0);
  mod = itom(0);
  one = itom(1);
  two = itom(2);

  mdiv(root, mnumprocs, part, junk);

  if (rank == 0) {
    bot = itom(3);
  } else {
    bot = itom(0);
    mult(part, mrank, bot);

    mdiv(bot, two, junk, mod);
    if (mcmp(mod, one) != 0) madd(bot, one, bot);
  }

  if (rank == numprocs-1) {
    top = itom(0);
    madd(root,one,top);
  } else {
    top = itom(0);
    madd(bot, part, top);
  }

  res = factorloop(mnum, bot, top);

  mfree(junk);
  mfree(one);
  mfree(two);
  mfree(mod);
  mfree(mnum);
  mfree(mrank);
  mfree(mnumprocs);
  mfree(root);
  mfree(part);
  mfree(bot);
  mfree(top);

  return res;

}



int main (int argc, char *argv[]) {

  MINT *tp;
  int i, myid, numprocs;
  char buf[1024];
  MPI_Status status;
  double startwtime, endwtime;

  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD,&myid);
  MPI_Comm_size(MPI_COMM_WORLD,&numprocs);

  if (myid == 0) {

    // This process here is the control process

    printf("Please enter a number you wish to factor:\n");
    if (fgets(buf, sizeof(buf), stdin) == NULL) exit(0);
    if (!isdigit(*buf)) {
      printf("Invalid number!\n");
      exit(0);
    }

    for(i=0; isdigit(buf[i]); i++)   ;
    buf[i] = '\0';

    for (i=1; i<numprocs; i++)
      MPI_Send(buf, strlen(buf)+1, MPI_CHAR, i, 0, MPI_COMM_WORLD);

    memset(buf, '\0', sizeof(buf));

    

    /* Start the timer */
    startwtime = MPI_Wtime();

    for(i=1; i<numprocs; i++) {

      MPI_Recv(buf, sizeof(buf), MPI_CHAR, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, &status);
      if (*buf != 'P') break;

    }

    /* Stop the timer */
    endwtime = MPI_Wtime();


    if (*buf == 'P') {
      printf("Your number was found to be prime...\n");
    } else {
      printf("A factor was found on node %d: %s\n", status.MPI_SOURCE, buf);
    }
    printf("It took %f seconds to compute this.\n", endwtime-startwtime);


    MPI_Finalize();

    exit(0);
  }


  // All processes here are work-horse processes:

  memset(buf, '\0', sizeof(buf));

  MPI_Recv(buf, sizeof(buf), MPI_CHAR, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, &status);

  tp = dofactor(buf, myid-1, numprocs-1);

  if (tp == NULL) {
    MPI_Send("P", 2, MPI_CHAR, 0, 0, MPI_COMM_WORLD);
  } else {
    minttoa(tp, buf, sizeof(buf));
    MPI_Send(buf, strlen(buf)+1, MPI_CHAR, 0, 0, MPI_COMM_WORLD);
  }

  MPI_Finalize();

  return 0;

}