/*
	Title: Find a near-minimax approximation
	
	This program calculates a near-minimax approximation
	to a given function 'fcn' on the interval [-1,1].
	The approximation is based on interpolation at the
	Chebyshev zeros on [-1,1].  The program then
	estimates the maximum error by calculating the error
	at a large sample of evenly spaced points on [-1,1].
	The program can be easily modified to print the
	divided differences and the errors. This is taken
	from Section 6.3.
*/

#include <stdio.h>
#include <math.h>

#define	MAX_N	20

/*  The following statement defines the function to be
	approximated. It can also be specified using a
	function subprogram.
*/
     
#define	fcn(x)	exp(x)


void divdif(int n, double x[MAX_N+1], double f[MAX_N+1]);

      
main()
{
	const double zero = 0.0, one = 1.0, p02 = .02;
	double fx[MAX_N+1], x[MAX_N+1];
	double pi, h, z, p, error, errmax;
	int i, j, n;

	
	/*  Initialize variables  */
      
	pi = 4*atan(one);
	printf("\n\n What is the degree n? : ");
	scanf("%d", &n);
	h = pi/(2*(n+1));

	/*  Set up interpolation nodes and corresponding
		function values */

	for (i=0; i <= n; i++)  
	{
		x[i] = cos((2*i+1)*h);
		fx[i] = fcn(x[i]);
	}

	/*  Calculate the divided differences for the
		interpolation polynomial  */
      
	divdif(n,x,fx);

	/* Calculate the maximum interpolation error at 101
		evenly spaced points on [-1,1].  
	*/
      
	errmax = zero;
	for (i=0; i <= 100; i++)  
	{
		z = -one + p02*i;
		
		/*  Use nested multiplication to evaluate the
			interpolation polynomial at z. 
		*/
        	
		p = fx[n];
		for (j=n-1; j >= 0; j--) 
			p = fx[j] + (z - x[j])*p;
		error = fcn(z) - p;
		if (fabs(error) > errmax)
			errmax = fabs(error);
	}

	
	/*  Print maximum error  */

	printf("\n For n = %d", n);
	printf("\n Maximum error = %e\n\n", errmax);
}



void divdif(int n, double x[MAX_N+1], double f[MAX_N+1])
{

/*
	Input: (1) 'n' denotes a positive integer.
	(2) 'x' denotes a vector of points x(0),...,x(n).
	(3) 'f' denotes a vector of values of some function
		evaluated at the nodes in x.
	Output: The vector f is converted to a vector of divided
	differences: f(i)=f[x(0),...,x(i)], i=0,...,n
*/

	int i,j;
 
	for(i=1; i <= n; i++)
		for(j=n; j >= i; j--)
			f[j] = (f[j] - f[j-1])/(x[j] - x[j-i]);

}