#include <iostream>
#include <mpi.h>
#include <stdlib.h>
#include <iomanip>
#include <math.h>


int main(int argc, char *argv[]) {
        double k, timeIteration;
        double h, dt;
        const double T = 0.0001; //Last moment of time

        int myrank, size; // MPI Initializings
        MPI_Status Status;

        double begin,end;

        MPI_Init(&argc, &argv);
        begin = MPI_Wtime();

        MPI_Comm_size(MPI_COMM_WORLD, &size);
        MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

        int pointsCount;
        int partSize;
        int lastPart;

        if (myrank == 0) {
                pointsCount = 50000;
        }

        MPI_Bcast(&pointsCount, 1, MPI_INT, 0, MPI_COMM_WORLD);

        k = 1;
        h = 1 / double(pointsCount);
        dt = h * h / (2 * k);
        partSize = pointsCount / size;
        lastPart = pointsCount % size;
        if (lastPart != 0 && myrank < lastPart) {
                ++partSize;
        }

        double temperatureArray[partSize];
        for (int i = 0; i < partSize; ++i) {
                temperatureArray[i] = 1;
        }
        if (myrank == 0) {
                temperatureArray[0] = 0;
        }
        if (myrank == size - 1) {
                temperatureArray[partSize - 1] = 0;
        }

        double leftElement = 0;
        double rightElement = 0;
        timeIteration = dt;

        while (timeIteration <= T) {
        double newTemperatureArray[partSize];
        newTemperatureArray[0] = 0;
        newTemperatureArray[partSize - 1] = 0;
                if (myrank != size - 1) {
                        MPI_Sendrecv(&temperatureArray[partSize - 1], 1, MPI_DOUBLE, myrank + 1, myrank, // передача
                                &rightElement, 1, MPI_DOUBLE, myrank + 1, MPI_ANY_TAG, // прием
                                MPI_COMM_WORLD, &Status);
                }

                if (myrank != 0) {
                        MPI_Sendrecv(&temperatureArray[0], 1, MPI_DOUBLE, myrank - 1, myrank, // передача
                                &leftElement, 1, MPI_DOUBLE, myrank - 1, MPI_ANY_TAG, // прием
                                MPI_COMM_WORLD, &Status);
                }

        if (myrank != 0) {
                    newTemperatureArray[0] = k * (temperatureArray[1] - 2 * temperatureArray[0] + leftElement) / (h * h) * dt + temperatureArray[0];
        }
                for (int i = 1; i < partSize - 1; ++i) {
                        newTemperatureArray[i] = k * (temperatureArray[i + 1] - 2 * temperatureArray[i] + temperatureArray[i - 1]) / (h * h) * dt + temperatureArray[i];
                }
        if (myrank != size - 1) {
                    newTemperatureArray[partSize - 1] = k * (rightElement - 2 * temperatureArray[partSize - 1] + temperatureArray[partSize - 2]) / (h * h) * dt + temperatureArray[partSize - 1];
        }

        for (int i = 0; i < partSize; ++i) {
            temperatureArray[i] = newTemperatureArray[i];
        }
                timeIteration += dt;
        }

        MPI_Barrier (MPI_COMM_WORLD);

        if (myrank == 0) {
                int offset = 0;
                double resultArray[pointsCount];
                for (int i = 0; i < partSize; ++i) {
                        resultArray[i] = temperatureArray[i];
                }
                offset += partSize;
                for (int i = 1; i < size; ++i) {
                        MPI_Probe(i, MPI_ANY_TAG, MPI_COMM_WORLD, &Status);
                        MPI_Get_count(&Status, MPI_DOUBLE, &partSize);

                        MPI_Recv(&resultArray[offset], partSize, MPI_DOUBLE, i, MPI_ANY_TAG, MPI_COMM_WORLD, &Status);
                        offset += partSize;
                }

                end = MPI_Wtime();

                //Counting by formula
                const double precision = 1e-12;
                double prevIteration = 0;
                double currentIteration = 0;
                double difference = 0;
                double formulaResultArray[11];
                formulaResultArray[0] = 0;
                formulaResultArray[10] = 0;

                //Так как sin(0) = 0 и sin(Pi * (2m + 1)) = 0 о считать значение
                //в начальной и конечной точках не имеет смысла
                for (int i = 1; i < 10; ++i) {
                        double point = (pointsCount / 10 * i) * h;
                        currentIteration = 0;
                        prevIteration = 0;
                        difference = 1e-5;
                        double iteration = 0;

                        while (difference > precision) {
                                prevIteration = currentIteration;
                                difference = exp(-1 * k * M_PI * M_PI * (2 * iteration + 1) * (2 * iteration + 1) * T / (1 * 1)) / (2 * iteration + 1) * sin(M_PI * (2 * iteration + 1) * point / 1);

                                currentIteration = prevIteration + difference;
                                ++iteration;
                        }
                        formulaResultArray[i] = 4 * currentIteration / M_PI;
                }


                std::cout <<  "1 point: " << resultArray[pointsCount - 1] << "    formula: 0" << std::endl;
                for (int i = 1; i < 10; ++i) {
                        int number = pointsCount / 10 * i - 1;
                        std::cout << i << " point: " << resultArray[number] << "    formula: " << formulaResultArray[i] << std::endl;
                }
        std::cout <<  "11 point: " << resultArray[pointsCount - 1] << "    formula: 0"<< std::endl;

                double time = end - begin;
                std::cout <<  std::setprecision(3) << "Program worked: " << time << " seconds" << std::endl;


        }
        else {
                MPI_Send(temperatureArray, partSize, MPI_DOUBLE, 0, myrank, MPI_COMM_WORLD);
        }

    MPI_Finalize();
        return 0;
}