#include "stdio.h"
#include "math.h"
#include "mpi.h"

int main(int argc, char *argv[]) 
{
/* mpi control declarations */
   int numprocs, numslaves, master, slave, myrank, tag;
   MPI_Status status;  
/* other declarations */   
   int indices[2], n, termsperproc, numterms, lastslaveterm;
   double nn, zetasum, partialsum, s;     
/* universal parameter and data initalizations */
   master = 0;
   numterms = 100000; 
   s = 2.0;    
/* mpi driver */
   tag =1;
   MPI_Init(&argc,&argv);
   MPI_Comm_size( MPI_COMM_WORLD, &numprocs );
   MPI_Comm_rank( MPI_COMM_WORLD, &myrank );
   numslaves = numprocs-1;
   if (myrank == master)
   /* master node code */
   { termsperproc = numterms/numslaves;
     lastslaveterm = termsperproc*numslaves;
     for (slave = 1; slave <= numslaves; slave=slave+1)
     {  indices[0] = termsperproc*(slave-1)+1;
        indices[1] = termsperproc*slave;
        MPI_Send(&indices,2,MPI_INT,slave,tag,MPI_COMM_WORLD);
     }
     zetasum = 0.0;   
     for (n = lastslaveterm+1; n <= numterms; n=n+1)
     {  nn = n;
        zetasum = zetasum + exp(-s*log(nn)); 
     }
     printf("tailsum = %.16f \n", zetasum);    
     for (slave = 1; slave <= numslaves; slave=slave+1)
     {  MPI_Recv(&partialsum,1,MPI_DOUBLE,slave,tag,MPI_COMM_WORLD,&status);
        zetasum = zetasum + partialsum; 
        printf("slave = %d, partial sum = %.16f, zetasum = %.16f \n", slave, partialsum, zetasum);
     }           
   }
   else
   /* slave node code */     
   {  MPI_Recv(&indices,2,MPI_INT,master,tag,MPI_COMM_WORLD,&status);
      partialsum = 0;
      for (n = indices[0]; n <= indices[1];n = n+1)
      {  nn = n;
         partialsum = partialsum + exp(-s*log(nn));
      } 
      MPI_Send(&partialsum,1,MPI_DOUBLE,master,tag,MPI_COMM_WORLD);  
   }  
/* mpi cleanup */
   if (myrank == master)
      printf("zetasum = %.16f \n", zetasum);
   MPI_Finalize();  
   return 0;
}