/*
Some not so big primes
 117773,524287, 612613,1300051, 2147483647, 2305843009213693951,618970019642690137449562111,
162259276829213363391578010288127, 170141183460469231731687303715884105727
The largest unsigned long (32 bits) 4,294,967,295
Assumptions: 
The number "n" that we are factoring is such that sqroot(sqroot(n)) is an unsigned long. This 
implies that n is limited to somewhere under 40 decimal places. Here is one "small" test product 61746852851
and another 321185031931, 796428143263, and another 1125897758834689
*/


#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <signal.h>
#include <string.h>
#include "mpi.h"
#include "gmp.h"

int main(int argc,char *argv[]){
    int done = 0,master=0,boot_builder=1,true=1,first_time=1,start_factor_loop=1,ack=1; 
    int myid, numprocs, i,j,mpi_error;
    char *block_buffer,*boot_buffer,n_buffer[1000],block_floor[1000],block_roof[1000],s[100];
    mpz_t n,sqroot,sqsqroot, block_base, isprime,aprime, test_value,boot_base,processor_block;
    unsigned long boot_top,block_top,long_i,long_j; 
    double startwtime, endwtime;
    FILE *file_pointer;
    char file_name[50];
    int  namelen;
    char *zero_string="0";
    char processor_name[MPI_MAX_PROCESSOR_NAME];
    MPI_Status mpistatus;
    MPI_Init(&argc,&argv);
    MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
    MPI_Comm_rank(MPI_COMM_WORLD,&myid);
    MPI_Get_processor_name(processor_name,&namelen);
    mpz_init_set_ui(boot_base,(unsigned long)2);
    sprintf(file_name,"iofile%d",myid);
    file_pointer=fopen(file_name,"w");
    fprintf(file_pointer,"iofile%d\n",myid);
    fflush(file_pointer);
    if (myid == master) {
        printf("\nDo you want to factor another?\n 1= no \n 0=yes\n ");
        fflush(stdout);
        scanf("%d",&done);
        if (!done) {
            printf("\nPlease input the number to factor\n");
            fflush(stdout);
            scanf("%s",n_buffer);
            if (first_time) {
                mpz_init(n);        //the number I know is the product of two primes
                mpz_init(sqroot);   // I only need to know the primes l.e. sqroot n.
                mpz_init(sqsqroot); // To get those primes I do a bootstrap sieve using primes
                                    // l.e. sqroot(sqroot(n))
                mpz_init(processor_block);   //each processor is asked to compute primes in this block
                mpz_init(test_value);      // just a gmp variable .
                mpz_init(isprime);         //   same
                mpz_init(aprime);           //   same
                first_time=0;
            } else {
                free(boot_buffer);
            }

            /*get the boot strap processor to work*/
            start_factor_loop=1; 
            mpz_set_str (n,n_buffer, 10);
            mpz_sqrt (sqroot, n); //greatest int in square root
            mpz_sqrt (sqsqroot, sqroot); //greatest int in square root of square root
            boot_top=mpz_get_ui(sqsqroot);
            MPI_Send(&boot_top, 1, MPI_LONG,boot_builder,0, MPI_COMM_WORLD);
            /**********************************************************/
            /*Get the processors computing blocks of primes to work*/
            mpz_sub(processor_block,sqroot,sqsqroot);
            mpz_get_str(s,10,processor_block);
            mpz_fdiv_q_ui(processor_block,processor_block,(unsigned long)(numprocs-2));
            mpz_get_str(s,10,processor_block);
            printf("processor_block size=%s\n",s);
            fflush(stdout);
            for (i=2;i<numprocs;i++) {
                mpz_mul_ui(test_value,processor_block,(unsigned long)(i-2));//test_value=(i-2)*processor_block
                mpz_add(block_base, test_value,sqsqroot);// add the base to find where block starts.
                mpz_get_str(block_floor,10, block_base);
                printf("block floor %s\n",block_floor);
                fflush(stdout);
                if (i==(numprocs-1))
                    mpz_get_str(block_roof,10, sqroot);
                else {

                    mpz_add(block_base, block_base,processor_block); //add one more.
                    mpz_get_str(block_roof,10, block_base); //"block_base" now points to "block_roof" 
                }
                printf("block roof %s\n",block_roof);
                fflush(stdout);
                MPI_Send(block_floor, 1000, MPI_CHAR,i,0, MPI_COMM_WORLD);
                MPI_Send(block_roof, 1000, MPI_CHAR,i,0, MPI_COMM_WORLD);
            }
            /********************************************/
            /*Receive primes and test from the various blocks...beginning with the boot block*/
            MPI_Recv(&boot_top,1,MPI_LONG,boot_builder,0, MPI_COMM_WORLD,&mpistatus);
            boot_buffer=malloc(boot_top*sizeof(char));
            MPI_Recv(boot_buffer, (int)boot_top, MPI_CHAR,boot_builder,0 ,MPI_COMM_WORLD,&mpistatus);
            for (j=0;j<(int)boot_top;j++) {
                if (boot_buffer[j]==1) {
                    mpz_add_ui(aprime,boot_base,(unsigned long)j);
                    mpz_fdiv_r(test_value,n,aprime);
                    if (mpz_cmp_ui(test_value,(unsigned long)0)==0) {
                        mpz_out_str(stdout,10,n);printf(" equals  ");
                        mpz_out_str(stdout,10,aprime);printf(" times ");
                        mpz_fdiv_q(test_value,n,aprime);
                        mpz_out_str(stdout,10,test_value);printf(" \n ");

                    }
                }
            }


            /*Now the other blocks*/
            for (i=2;i<numprocs;i++) {
                mpz_mul_ui(test_value,processor_block,(unsigned long)(i-2));
                mpz_add(block_base, test_value,sqsqroot);
                if (i==(numprocs-1))
                    mpz_set(test_value,sqroot);
                else
                    mpz_add(test_value, block_base,processor_block);
                mpz_sub(test_value,test_value,block_base);
                block_top=mpz_get_ui(test_value);
                if (start_factor_loop) {
                    start_factor_loop=0;
                } else {
                    free(block_buffer);
                }   
                block_buffer=malloc((int)block_top*sizeof(char));              
                MPI_Recv(block_buffer,(int)block_top, MPI_CHAR,i,0 ,MPI_COMM_WORLD,&mpistatus);//(int)block_top
                //     MPI_Send(&ack, 1, MPI_INT,i,0 ,MPI_COMM_WORLD);
                for (j=0;j<(int)block_top;j++) {
                    if (block_buffer[j]==1) {
                        mpz_add_ui(aprime,block_base,(unsigned long)j);
                        mpz_fdiv_r(test_value,n,aprime);
                        if (mpz_cmp_ui(test_value,(unsigned long)0)==0) {
                            mpz_out_str(stdout,10,n);printf(" equals  ");
                            mpz_out_str(stdout,10,aprime);printf(" times ");
                            mpz_fdiv_q(test_value,n,aprime);
                            mpz_out_str(stdout,10,test_value);printf(" \n ");

                        }
                    }
                }
            }

        }  //done
    }  //master_block

    if (myid==boot_builder) {
        if (first_time) {
            mpz_init(test_value);
            mpz_init(isprime);
            first_time=0;
        } else {
            free(boot_buffer);
        }
        MPI_Recv(&boot_top, 1, MPI_LONG,master,0 ,MPI_COMM_WORLD,&mpistatus);
        boot_buffer=malloc((int)boot_top*sizeof(char));
        for (long_i=0;long_i<boot_top;long_i++)boot_buffer[long_i]=1;
        for (long_i=0;long_i<boot_top;long_i++) {
            if (boot_buffer[long_i]==1) {
                mpz_add_ui(isprime,boot_base,long_i);
                for (long_j=long_i+1;long_j<boot_top;long_j++) {
                    if (boot_buffer[long_j]==1) {
                        mpz_add_ui(test_value,boot_base,long_j);
                        mpz_fdiv_r(test_value,test_value,isprime);
                        if (mpz_get_ui(test_value)==(unsigned long)0)
                            boot_buffer[long_j]=0;
                    } //if
                } //for

            }//if
        }//for
        MPI_Send( &boot_top,1,MPI_LONG,master,0, MPI_COMM_WORLD);
        MPI_Send(boot_buffer, boot_top, MPI_CHAR,master,0, MPI_COMM_WORLD);
        for (i=2;i<numprocs;i++) {
            MPI_Send( &boot_top,1,MPI_LONG,i,0, MPI_COMM_WORLD);
            MPI_Send(boot_buffer, boot_top, MPI_CHAR,i,0, MPI_COMM_WORLD);
            // MPI_Recv(&ack, 1, MPI_INT,i,0 ,MPI_COMM_WORLD,&mpistatus);
        }

    }
    if ((myid!=boot_builder)&&(myid!=master)) {
        if (first_time) {
            mpz_init (isprime);
            mpz_init (test_value);
            mpz_init (aprime);
            first_time=0;
        } else {
            free(block_buffer);
        }

        MPI_Recv(&boot_top,1,MPI_LONG,boot_builder,0, MPI_COMM_WORLD,&mpistatus);
        boot_buffer=malloc(boot_top*sizeof(char));
        MPI_Recv(boot_buffer, (int)boot_top, MPI_CHAR,boot_builder,0 ,MPI_COMM_WORLD,&mpistatus);
        //  MPI_Send(&ack,1,MPI_INT,boot_builder,0, MPI_COMM_WORLD);
        MPI_Recv(block_floor, 1000, MPI_CHAR,master,0, MPI_COMM_WORLD,&mpistatus);
        MPI_Recv(block_roof, 1000, MPI_CHAR,master,0, MPI_COMM_WORLD,&mpistatus);
        mpz_init_set_str(block_base,block_floor,10);
        mpz_init_set_str(test_value,block_roof,10);
        mpz_sub(test_value,test_value,block_base);
        block_top=mpz_get_ui(test_value);
        block_buffer=malloc(block_top*sizeof(char));
        for (long_i=0;long_i<block_top;long_i++)block_buffer[long_i]=1;
        for (i=0;i<boot_top;i++) {
            if (boot_buffer[i]==1) {
                mpz_add_ui(aprime,boot_base,(unsigned long)i);
                for (long_j=0;long_j<block_top;long_j++) {
                    if (block_buffer[long_j]==1) {
                        mpz_add_ui(test_value,block_base,long_i);
                        mpz_fdiv_r(test_value,test_value,aprime);
                        if (mpz_get_ui(test_value)==(unsigned long)0)
                            block_buffer[long_i]=0;
                    } //if
                } //for
            }//if
        }//for 
        MPI_Send(block_buffer, (int)block_top, MPI_CHAR,master,0, MPI_COMM_WORLD);
        //  MPI_Recv(&ack, 1, MPI_INT,master,0 ,MPI_COMM_WORLD,&mpistatus);
    }
    MPI_Bcast(&ack,1, MPI_INT,master,MPI_COMM_WORLD);
    mpi_error=MPI_Finalize();
    if (MPI_SUCCESS!=mpi_error)
        return mpi_error;
    else
        return 0;

}