/* remess.c */
#include "remes.h"

/* Search subroutine */

remess()
{
double a, b, q, xm, ym, xn, yn, xx0, xx1;
int i, j, meq, emsign, ensign, steps;
double approx(), func(), geterr();

/* Search for maxima of error */
eclose = 1e30;
farther = 0.0;
j = 0;
meq = neq;

if( (d > 0) || ((config & ZER) != 0) )
	{
	j = 1;
	meq += 1;
	}
xx0 = apstrt;

for( i=0; i<meq; i++ )
	{
	steps = 0;
	if( (d > 0) || ((config & ZER) != 0) )
		xx1 = xx[i]; /* Next zero */
	else
		xx1 = mm[i+1]; /* Next maximum */
	if( i == meq-1 )
		xx1 = apend;
	xm = mm[i];
	ym = geterr(xm);
	emsign = esign; /* Sign of error */
	q = step[i];
	xn = xm + q;
/* Cannot skip over adjacent boundaries */
	if( xn < xx0 )
		goto revers;
	if( xn >= xx1 )
		goto revers;
	yn = geterr(xn);
	ensign = esign;
	if( yn < ym )
		{
revers:
		q = -q;
		xn = xm;
		yn = ym;
		ensign = emsign;
		}

/* March until error becomes smaller. */

	while( yn >= ym )
		{
		if( ++steps > 10 )
			goto tsear;
		ym = yn;
		xm = xn;
		emsign = ensign;
		a = xm + q;
		if( a == xm )
			goto tsear;
/* Must not skip over the zeros either side. */
		if( a <= xx0 )
			goto tsear;
		if( a >= xx1 )
			goto tsear;
		xn = a;
		yn = geterr(xn);
		ensign = esign;
		}

tsear:
	mm[i] = xm; /* Position of maximum */
	yy[i] = ym; /* Value of maximum */

	if( ym == 0.0 )
		printf( "yy[%d] = 0\n", i );
	if( eclose > ym )
		eclose = ym;

	if( farther < ym )
		farther = ym;

 /* No denominator */
	if( (d == 0) && ((config & ZER) == 0) )
		xx[i] = xm;
/* Walk to next zero. */
	if( (d > 0) || ((config & ZER) != 0) )
		xx0 = xx1;
	else
		xx0 = 0.5*(xm + xx1);
	} /* end of search loop */


/* Decrease step size if error spread increased. */
q = (farther - eclose);
/* Relative error spread */
if( eclose != 0.0 )
	q /= eclose;
if( q >= spread )
	{ /* Spread is increasing; decrease step size. */
	if( config & TRUNC )
		delta *= 0.875;
	else
		delta *= 0.5;
	printf( "delta = %.4E\n", delta );
	}
 spread = q;

printf( "peak error = %.4E, relative error spread = %.4E\n",
		farther, spread );

/* Calculate new step sizes */

if( (d == 0) && ((config & ZER) == 0) )
	{
	if( spread < 0.25 )
		q = 0.0625;
	else
		q = 0.5;

	for( i=0; i<neq; i++ )
		step[i] *= q;
	}
else
	{

	for( i=0; i<neq; i++ )
		{
		q = yy[i+1];
		if( q != 0.0 )
			q = yy[i]/q  - 1.0;
		else
			q = 0.0625;
		if( q > 0.25 )
			q = 0.25;
		q *= mm[i+1] - mm[i];
		step[i] = q * delta;
		}
	step[neq] = step[neq-1];


/* Insert new locations for the zeros. */
	for( i=0; i<neq; i++ )
		{
		xm = xx[i] - step[i];
		if( xm <= apstrt )
			continue;
		if( xm >= apend )
			continue;
		if( xm <= mm[i] )
			{
			printf( "xx[%d] < mm\n", i );
			xm = 0.5 * (mm[i] + xx[i]);
			}
		if( xm >= mm[i+1] )
			{
			printf( "xx[%d] > mm\n", i );
			xm = 0.5 * (mm[i+1] + xx[i]);
			}
		xx[i] = xm;
		}
	}
sdone:	;
}