/* remesa.c */

#include "remes.h"
#define PI 3.14159265358979323846

/* Subroutine to get operator instructions. */

remesa()
{
char s[40];
double a, b, c, q, r;
int i, k, ncheb;
double func(), cos(), sqrt();

/* Ask for error criterion */
printf( "Relative error (y or n) ? " );
gets( s ); /* Read in operator's response. */
relerr = 0;
if( s[0] == 'y' ) /* "y" means yes. */
	relerr = 1;

/* Enable Automatic search. */
search = 1;

getconf:

printf( "Configuration word = " );
pconfig();
printf( " ? " );
gets(s);
if( s[0] != '\0' )
	{
/* Display a menu of configuration bits */
	if( s[0] == '?' )
		{
		k = config;
		config = 0xffff;
		pconfig();
		config = k;
		goto getconf;
		}
	else
		{
		sscanf( s, "%x", &config );
		pconfig();
		}
	printf( "\n" );
	}

printf( "Degree of numerator polynomial? " );
gets( s );		/* read line */
sscanf( s, "%d", &n );	/* decode characters */

printf( "Degree of denominator polynomial? " );
gets( s );
sscanf( s, "%d", &d );

printf ( "Start of approximation interval ? " );
gets( s );
sscanf( s, "%lf", &apstrt );

getwid:
printf ( "Width of approximation interval ? " );
gets( s );
sscanf( s, "%lf", &apwidt );
if( apwidt <= 0.0 )
	{
	puts( "? must be > 0" );
	goto getwid;
	}
apend = apstrt + apwidt;
nd1 = n + d + 1;
spread = 1.0e37;
iter = 0;
askitr = 1;

/* Supply initial guesses for solution points. */

if( (d == 0) && ((config & ZER) == 0) )
	{ /* there is no denominator polynomial */
	neq = n + 2; /* The number of equations to solve */
	mm[neq] = apend;
	}
else
	{
	neq = nd1;
	}
ncheb = nd1; /* Degree of Chebyshev error estimate */

/* Find ncheb+1 extrema of Chebyshev polynomial */
a = ncheb;
mm[0] = apstrt;
for( i=1; i<ncheb; i++ )
	{
	r = 0.5 * (1.0 - cos( (PI * i) / a ) );
	if( config & SQL )
		r = r * sqrt(r);
	else if( config & SQH )
		r = sqrt(r);
	mm[i] = apstrt + r * apwidt;
	}
mm[ncheb] = apend;

if( (d == 0) && ((config & ZER) == 0) )
	{ /* The operative locations are maxima. */
	for( i=0; i<=neq; i++ )
		xx[i] = mm[i];
	}
else
	{ /* Zeros of Chebyshev polynomial */
	a = 2.0 * ncheb;
	for( i=0; i<=ncheb; i++ )
		{
		r = 0.5 * (1.0 - cos( PI * (2*i+1) / a ) );
/* Squeeze toward low end of approximation interval. */
		if( config & SQL )
			r = r * sqrt(r);
/* Squeeze toward high end. */
		else if( config & SQH )
			r = sqrt(r);
		xx[i] = apstrt + r * apwidt;
		}
/* Deviation at solution points */
	dev = 0.0;
	}

/* If form is xR(x), avoid x = 0 */
if( config & (XPX | X2PX) )
	{
	if( config & CW )
		qx = 1.0e-15 + xx[0]/2.0;
	else
		qx =  1.0e-15 + apstrt * (1.0 + 1.0e-15);
	mm[0] = qx;
	printf( "mm[0] = %.15E\n", qx );
	if( (d == 0) && ((config & ZER) == 0) )
		xx[0] = qx;
	}
/* Initialize step sizes */
stpini();
}



/* Subroutine to initialize step sizes. */
stpini()
{
int i;

if( search )
	{
	xx[neq+1] = apend;
	delta = 0.25;
	if( (d > 0) || ((config & ZER) != 0) )
		{
		step[0] = xx[0] - apstrt;
		for( i=1; i<neq; i++ )
			step[i] = xx[i] - xx[i-1];
		}
	else
		{
		step[0] = 0.0625 * (xx[1] - xx[0]);
		for( i=1; i<neq; i++ )
			step[i] = 0.0625 * (xx[i+1] - xx[i]);
		}
	step[neq] = step[neq-1];
	}
}


/* Bits of configuration word */
#define NCNFG 11

static char *cnfmsg[NCNFG] = {
"PXSQ",
"XPX",
"PADE",
"CW",
"SQL",
"ZER",
"X2PX",
"SQH",
"SPECIAL",
"TRUNC",
"PXCU"
};

/* Subroutine to Display the configuration word */

pconfig()
{
int i, k;

k = 1;
for( i=0; i<NCNFG; i++ )
	{
	if( k & config )
		printf( "%x:%s ", k, cnfmsg[i] );
	k <<= 1;
	}
printf( "hex value = %x ", config );
}