#define M_PI 3.1415926

double fn_gaussian(double mean, double sd, double x) {

    return (1.0 / (sd * sqrt(2.0 * M_PI) ) *
	    exp(- pow((x-mean),2.0) / (2.0*sd*sd)));
}

#define GAUSS_END 3.0

/*************************************************************/
int all_fill_array_gauss(int M, int A[PROCS]::[M], 
			  int N[PROCS]::, int total_n)
/*************************************************************/
{
    register int
	i, j;

    double
	p,
	x_seg,
	x[PROCS],
	g[PROCS],
	tot;

    barrier();

    tot   = 0.0;
    p     = (double)PROCS;
    x_seg = 2.0 * GAUSS_END / p;

    for (i=0 ; i<PROCS ; i++) {
	x[i] = - GAUSS_END + (x_seg/2.0) + ((double)i * x_seg);
        g[i] = fn_gaussian(0.0, 1.0, x[i]);
	tot  += g[i];
    }

    j = (int)ceil((((double)total_n * g[MYPROC]) / tot));
    	
    N[MYPROC] = j;

    if (j > M)
	fprintf(stderr,"ERROR: j: %d  M: %d\n",j,M);
    
    for (i=0 ; i<j ; i++)
	A[MYPROC][i] = create_data(); /* random() */

    return all_reduce_to_all_add(N[MYPROC]);
}