/* Alignment macros */
#define SPE_ALIGNMENT 16
#define SPE_ALIGNMENT_FULL 128
#define SPE_ALIGN __attribute__((aligned(16)))
#define SPE_ALIGN_FULL __attribute__((aligned(128)))

/* Rounding macros */

/* Go _up_ to the next address based on an alignment */
#define ROUND_UP_ALIGN(value, alignment) (((value) + \
 ((alignment) - 1))&(~((alignment)-1)))

/* Go _down_ to the next address based on an alignment */
#define ROUND_DOWN_ALIGN(value, alignment) ((value)&(~((alignment) - 1)))

/* Divide and round up */
#define DIVIDE_ROUND_UP(value, divisor) (((value) + (divisor) - 1) / (divisor))

/* In case these aren't already defined */
#ifndef MIN
#define MIN(X,Y) (((X) < (Y)) ? (X) : (Y))
#endif
#ifndef MAX
#define MAX(X,Y) (((X) > (Y)) ? (X) : (Y))
#endif  

/* Simplify Dereferencing on the SPE (to deref other types, change the typecase at the */
/* beginning) */
#define SPE_DEREF_UINT32(base, offset) *((unsigned int *)(((char *)(base)) + (offset))) 

/* Give a 16-byte structure for "status" so that it can be used in an array without */
/* causing alignment errors */
typedef struct {
	int status SPE_ALIGN;
} spe_remote_function_status_t;

                
float spe_calculate_standard_deviation(int num_values, float *values_ea) {
	/* Previous Data Definitions Go Here */
	/* ..... */

	/* Calculate the number of iterations we are going to need */
	int num_values_remaining = num_values;
	float *current_values_ea = values_ea;

	while(num_values_remaining > 0) {
		int values_this_cycle = MIN(num_values_remaining, MAX_VALUES);

		/* Previous code goes here */
		/* Replace values_ea with current_values_ea */
		/* Replace num_values with values_this_cycle */
		/* ..... */

		num_values_remaining -= values_this_cycle;
		current_values_ea += values_this_cycle;

	}

	/* Rest of function from original */

	/* Finalize Computation */
	avg = sum / (float)num_values;
	variance = (sum_squares - (sum * avg)) / (float)num_values;
	std_dev = sqrt(variance);

	return std_dev;
}

                
while(num_values_remaining > 0) {
		int values_this_cycle = MIN(num_values_remaining, MAX_VALUES);

		/* Check to see if we the leading values are unaligned.  If so, go one */
            /* value at a time until they are aligned. */
		if(((int)current_values_ea % DMA_TRANSFER_UNIT) != 0) {
			values_this_cycle = 1;
		} else {
			/* Check to see if this is a size that is transferable in one block */
			if((values_this_cycle % DMA_TRANSFER_UNIT_FLOAT) != 0) {
				/* Check to see if we can transfer a larger sub-block in a */
                        /* single transfer */
				int tmp_num_values = values_this_cycle - (values_this_cycle %
                         DMA_TRANSFER_UNIT_FLOAT);
				if(tmp_num_values > 0) {
					/* If so, use this value */
					values_this_cycle = tmp_num_values;
				} else {
					/* Otherwise, just go one value at a time */
					values_this_cycle = 1;
				}
			}
		}

		/* Load values in from main memory pointer */

		/* If we are doing one at a time, then we need to check for */
            /* alignment issues */
		if(values_this_cycle == 1) {
			/* Unaligned addresses require similarly unaligned destinations - */
                  /* calculate the offset */
			int natural_offset = (((int)current_values_ea / sizeof(float)) %
                   DMA_TRANSFER_UNIT_FLOAT);

			/* Perform the Transfer */
			spu_mfcdma32(&ls_values[natural_offset],
                   (unsigned int)current_values_ea, values_this_cycle * sizeof(float),
                   DEFAULT_DMA_TAG, MFC_GET_CMD);
			spu_mfcstat(MFC_TAG_UPDATE_ALL);

			/* Move it back into position for processing */
			ls_values[0] = ls_values[natural_offset];
		} else {
			/* Regular properly-aligned, properly-sized transfer */
			spu_mfcdma32(ls_values, (unsigned int)current_values_ea,
                   values_this_cycle * sizeof(float), DEFAULT_DMA_TAG, MFC_GET_CMD);
			spu_mfcstat(MFC_TAG_UPDATE_ALL);
		}


		for(i = 0; i < values_this_cycle; i++) {
			sum += ls_values[i];
			sum_squares += ls_values[i]*ls_values[i];
		}

		num_values_remaining -= values_this_cycle;
		current_values_ea += values_this_cycle;
	}

#define DMA_TRANSFER_UNIT 16
#define DMA_TRANSFER_UNIT_FLOAT (DMA_TRANSFER_UNIT / sizeof(float))

                
#ifndef MY_MATH_SPE_H
#define MY_MATH_SPE_H

#include "speport.h"

typedef struct {
	int num_values SPE_ALIGN;
	float *values SPE_ALIGN;
	float result_sum SPE_ALIGN;
	float result_sum_squares SPE_ALIGN;
} spe_std_dev_params_t;

#endif   

 
void spe_calculate_sums(int num_values, float *values_ea, float *result_sum,
 float *result_sum_squares);
                  

 
void spe_calculate_sums(int num_values, float *values_ea, float *result_sum, float
 *result_sum_squares) { 
                  

 
/* Perform Task */
spe_calculate_sums(params.num_values, params.values, &params.result_sum,
&params.result_sum_squares);
                  

 
*result_sum = sum;
*result_sum_squares = sum_squares;
return;
                  

                
                        /* ... header files ... */
#ifndef USE_SPE
/* Non-SPE Version */
/* ... */
#else
/* SPE-specific includes */
#include "speport.h"
#include "my_math_spe.h"

#define NUM_SPE_THREADS 4
float calculate_standard_deviation(int num_values, float *values) {
	/* Keep SPE Thread Contexts Here */
	static spe_remote_function_ptr_t std_dev_func[NUM_SPE_THREADS] =
       {NULL, NULL, NULL, NULL};

	/* Parameters to pass to contexts */
	spe_std_dev_params_t params[NUM_SPE_THREADS] SPE_ALIGN;
	spe_remote_function_status_t status[NUM_SPE_THREADS] SPE_ALIGN;

	/* Working variables */
	int thread_idx;
	int workunit_size = ROUND_UP_ALIGN(DIVIDE_ROUND_UP(num_values, NUM_SPE_THREADS),
       (SPE_ALIGNMENT_FULL/sizeof(float)));
	int current_workunit_size;
	int remaining_values = num_values;
	float *current_value_position = values;
	float sum = 0.0, sum_squares = 0.0, avg, variance, std_dev;

	/* Start up SPE processes if this is our first run */
	if(std_dev_func[0] == NULL) {
		for(thread_idx = 0; thread_idx < NUM_SPE_THREADS; thread_idx++) {
			std_dev_func[thread_idx] = spe_remote_function_start("./spe_std_dev",
                   NULL);
			if(std_dev_func[thread_idx] == NULL) {
				fprintf(stderr, "Error starting thread!\n");
				exit(1);
			}
		}
	}

	/* Split up the Data and Call the Function */
	for(thread_idx = 0; thread_idx < NUM_SPE_THREADS; thread_idx++) {
		current_workunit_size = MIN(workunit_size, remaining_values);

		/* Make parameters */
		params[thread_idx].num_values = current_workunit_size;
		params[thread_idx].values = current_value_position;

		/* Call the function */
		if(spe_remote_call(std_dev_func[thread_idx], &params[thread_idx],
             SPE_RUN_NONBLOCK, &status[thread_idx]) < 0) {
			fprintf(stderr, "Error running function\n");
			exit(1);
		}

		/* Move counter and pointer for next function call */
		remaining_values -= current_workunit_size;
		current_value_position += current_workunit_size;
	}

	/* Gather the results */
	for(thread_idx = 0; thread_idx < NUM_SPE_THREADS; thread_idx++) {
		/* Wait for thread to complete */
		spe_wait_completion(&status[thread_idx], 0);

		/* Gather data from parameters */
		sum += params[thread_idx].result_sum;
		sum_squares += params[thread_idx].result_sum_squares;
	}
 
	/* Finalize calculation */
	avg = sum / (float)num_values;
	variance = (sum_squares - (sum * avg)) / (float)num_values;
	std_dev = sqrt(variance);

	/* Return result */
	return std_dev;
}
#endif