// 去除广播
#include "mpi.h"
#include <math.h>
#include <stdio.h>
#define MIN(a,b)  ((a)<(b)?(a):(b))
#define VALUE_TO_INDEX(x) ((x - 3) / 2)
#define INDEX_TO_VALUE(x) (x * 2 + 3)
typedef long long LL;

int main (int argc, char *argv[])
{
    LL     count;        /* Local prime count */
    double elapsed_time; /* Parallel execution time */
    LL     first;        /* Index of first multiple */
    LL     global_count; /* Global prime count */
    LL     high_value;   /* Highest value on this proc */
    LL     i;
    int    id;           /* Process ID number */
    LL     low_value;    /* Lowest value on this proc */
    char  *marked;       /* Portion of 2,...,'n' */
    LL     n;            /* Sieving from 2, ..., 'n' */
    int    p;            /* Number of processes */
    int    proc0_size;   /* Size of proc 0's subarray */
    LL     prime;        /* Current prime */
    LL     size;         /* Elements in 'marked' */
    LL     low_index;
    LL     high_index;
    char  *cand;         /* Get the prime */

    MPI_Init (&argc, &argv);

    /* Start the timer */

    MPI_Comm_rank (MPI_COMM_WORLD, &id);
    MPI_Comm_size (MPI_COMM_WORLD, &p);
    MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time = -MPI_Wtime();

    if (argc != 2) {
        if (!id) printf ("Command line: %s <m>\n", argv[0]);
        MPI_Finalize();
        exit (1);
    }

    n = atoll(argv[1]);

    /* Figure out this process's share of the array, as
        well as the integers represented by the first and
        last array elements */

    low_index = 0 + (LL) id * (n - 1) / p / 2;            // index为3开始的奇数映射的下标，3-->0，5-->1，依此类推
    high_index = -1 + (LL) (id + 1) * (n - 1) / p / 2;
    low_value = INDEX_TO_VALUE(low_index);
    high_value = INDEX_TO_VALUE(high_index);
    size = high_index - low_index + 1;

    /* Allocate this process's share of the array. */

    marked = (char *) malloc (size);

    if (marked == NULL) {
        printf ("Cannot allocate enough memory for marked\n");
        MPI_Finalize();
        exit (1);
    }

    for (i = 0; i < size; i++) marked[i] = 0;

    LL N = (LL) sqrt(n);
    cand = (char *) malloc (N + 1);

    if (cand == NULL) {
        printf ("Cannot allocate enough memory for cand\n");
        MPI_Finalize();
        exit (1);
    }

    for (i = 0; i < N + 1; i++) cand[i] = 0;

    prime = 3;
    do {
        if (prime * prime > low_value) {
            first = VALUE_TO_INDEX(prime * prime) - low_index;
        } else {
            if (!(low_value % prime)) first = 0;
            else {
                LL t = (prime - low_value % prime) + low_value;
                if (t % 2) first = VALUE_TO_INDEX(t) - low_index;
                else first = VALUE_TO_INDEX(t + prime) - low_index;
            }
        }

        for (i = first; i < size; i += prime) marked[i] = 1;

        // 进程单独计算下一个质数，消除广播
        for (i = prime; i <= N; i += prime) cand[i] = 1;
        do prime++; while (prime % 2 == 0 || cand[prime]);

    } while (prime * prime <= n);

    count = 0;
    for (i = 0; i < size; i++)
        if (!marked[i]) count++;
    
    if (p > 1) MPI_Reduce (&count, &global_count, 1, MPI_INT, MPI_SUM,
        0, MPI_COMM_WORLD);
    else global_count = count;

    /* Stop the timer */

    elapsed_time += MPI_Wtime();

    /* Print the results */

    if (!id) {
        printf ("There are %lld primes less than or equal to %lld\n",
            global_count + 1, n);
        printf ("SIEVE (%d) %10.6f\n", p, elapsed_time);
    }
    MPI_Finalize ();
    return 0;
}