/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:16 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_silup s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_silup.h"
		/* PARAMETER translations */
#define	H	.5e0
#define	NDEG1	2
#define	NDEG2	10
#define	NTAB	40
#define	NX	3
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int i, ndeg;
	float ans, eopt[1], xt[2], y, yt[NTAB];
	static float x[NX]={-.1e0,.1e0,8.3e0};
	static long iopt[3]={0,1,0};
	static long lup = 3;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Eopt = &eopt[0] - 1;
	long *const Iopt = &iopt[0] - 1;
	float *const X = &x[0] - 1;
	float *const Xt = &xt[0] - 1;
	float *const Yt = &yt[0] - 1;
		/* end of OFFSET VECTORS */
 
	/*>> 2001-05-22 DRSILU Krogh Minor change for making .f90 version.
	 *>> 1994-12-21 DRSILU  Krogh  Latest version.
	 *--S replaces "?": DR?ILU, ?ILUP
	 * Demonstration driver for SILUP.
	 * Given table of sin(x), x = 0, .5, 1., ..., interpolates for x =
	 * -.1, .1, and, 8.3 using polynomial degrees 2, 3, 4, ..., 10.  Errors
	 * are estimated in all cases.
	 * */
	/*                       Set IOPT to get an error estimate. */
 
	/* Compute the XT and YT tables */
	Xt[1] = 0.e0;
	Xt[2] = H;
	for (i = 1; i <= NTAB; i++)
	{
		Yt[i] = sinf( (float)( i - 1 )*H );
	}
	printf("IOP(1) NDEG    X         Y     Est. Error True Error\n");
	for (i = 1; i <= NX; i++)
	{
		ans = sinf( X[i] );
		for (ndeg = NDEG1; ndeg <= NDEG2; ndeg++)
		{
			silup( X[i], &y, NTAB, xt, yt, ndeg, &lup, iopt, eopt );
			printf("%5ld%5ld%9.4f%12.8f%10.2e%11.2e\n", Iopt[1], ndeg, X[i], y, Eopt[1], y - ans);
		}
		fputc( '\n', stdout );
	}
	exit(0);
} /* end of function */