/*
      	Title: This is a test of function sumprd.
 
      	The test solves problem 2(a) of Section 3.3.
*/


#include <stdio.h>

 
float sumprd(float *a, float *b,int n);

float a[5001], b[5001];


main()
{
	int n, j;
	float sum, truesum, error;
	
	n = 1;
	while (n>0) 
	{	
		printf("\nHow many terms(max=5000): ");
		scanf("%d", &n);
		if (n == 0) 
			return(0);
     	 
		for(j=1; j <= n; j++)  
		{
			a[j] = 1.0/j;
			b[j] = 1.0/(j+1);
		}
     	 
		sum = sumprd(a,b,n);
		truesum = ((double) n) / (double)(n+1);
		error = truesum-sum; 
     	 
		printf("\n\t true = %8.6e", truesum);
		printf("\n\t approx = %8.6e", sum);
		printf("\n\t error = %8.6e", error);
     				 
		sum = 0.0;
		for (j=1; j <= n; j++)
			sum = sum + a[j]*b[j];
		error = truesum - sum;  
		
		printf("\n\t Single Precision Accumulation");
		printf("\n\t approx = %8.6e", sum);
		printf("\n\t error = %8.6e \n", error);
		return (0);
	}
	return(0);
}
	

float sumprd(float *a, float *b,int n)
{

/*

      	This calculates the inner product
 
                        i=n
             sumprd=  sum   a(i)*b(i)
                        i=1
                     
      	The products and sums are done in double
      	precision, and the final result is converted
      	back to single precision.
*/

	double dsum = 0.0;
	float result;
	int i;

	for(i=1; i <= n; i++)
		dsum = dsum + ((double) a[i])*((double) b[i]);
	result = dsum;
	return(result);
}