/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:08 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_daccum s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_daccum.h"
/*     program DRDACCUM
 *>> 1996-06-18 DRDACCUM Krogh Special code for C conversion.
 *>> 1996-05-28 DRDACCUM Krogh Added external state. & moved up formats
 *>> 1995-09-15 DRDACCUM Krogh Remove '0' in format (again?)
 *>> 1994-10-19 DRDACCUM Krogh  Changes to use M77CON
 *>> 1994-08-09 DRDACCUM WVS Removed '0' in formats
 *>> 1991-11-20 DRDACCUM CLL Edited for Fortran 90.
 *>> 1987-12-09 DRDACCUM Lawson  Initial Code.
 *     Demonstration driver for DACCUM.
 *--D replaces "?": DR?ACCUM, ?ACCUM, ?HFTI, ?COPY, ?MPVAL
 *     ------------------------------------------------------------------ */
		/* PARAMETER translations */
#define	LDIM	(NMAX + 2)
#define	LPMAX	(NMAX + 2)
#define	MDATA	12
#define	NB	1
#define	NMAX	8
#define	TAU	1.0e-5
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int i, ip[NMAX], ir1, irow, j, krank, n, ncount, ndeg, nrows,
	 _i, _r;
	double a[NMAX][LDIM], b[LDIM], dof, r, rnorm[NB], sigfac, u, work[NMAX],
	 yfit;
	static double p[LPMAX];
	static double x[MDATA]={2.0e0,4.0e0,6.0e0,8.0e0,10.0e0,12.0e0,
	 14.0e0,16.0e0,18.0e0,20.0e0,22.0e0,24.0e0};
	static double y[MDATA]={2.2e0,4.0e0,5.0e0,4.6e0,2.8e0,2.7e0,3.8e0,
	 5.1e0,6.1e0,6.3e0,5.0e0,2.0e0};
	static int _aini = 1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const B = &b[0] - 1;
	long *const Ip = &ip[0] - 1;
	double *const P = &p[0] - 1;
	double *const Rnorm = &rnorm[0] - 1;
	double *const Work = &work[0] - 1;
	double *const X = &x[0] - 1;
	double *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		P[1] = 13.0e0;
		P[2] = 11.0e0;
		_aini = 0;
	}
 
	/*     ------------------------------------------------------------------ */
 
	/*     ------------------------------------------------------------------ */
	n = NMAX;
	ndeg = n - 1;
	ir1 = 1;
	nrows = 1;
 
	for (irow = 1; irow <= MDATA; irow++)
	{
		u = (X[irow] - P[1])/P[2];
		i = ir1;
		a[0][i - 1] = 1.;
		for (j = 2; j <= (ndeg + 1); j++)
		{
			a[j - 1][i - 1] = a[j - 2][i - 1]*u;
		}
		B[i] = Y[irow];
		daccum( (double*)a, LDIM, n, b, LDIM, NB, &ir1, nrows, &ncount );
	}
 
	printf("DRDACCUM..  Demo driver for DACCUM.\n");
	printf(" MDATA = %4d, NCOUNT = %4ld\n", MDATA, ncount);
	dhfti( (double*)a, LDIM, ir1 - 1, n, b, LDIM, NB, TAU, &krank,
	 rnorm, work, ip );
	printf(" KRANK = %4ld\n", krank);
 
	/*                      The following stmt does a type conversion. */
	dof = ncount - n;
	sigfac = Rnorm[1]/sqrt( dof );
	dcopy( n, b, 1, &P[3], 1 );
	/*++ Code for .C. is active */
   printf(
      "\n NDEG =%2ld          RNORM =%8.4f          SIGFAC =%8.4f",
       ndeg, Rnorm[1], sigfac );
   printf(
   "\n\n P(1),P(2) =         %15.5f%15.5f\n\n P(3),...,P(NDEG+3) =",
   p[0], p[1]);
   for (i = 2; i < (n + 2); i+=3){
      for (j = i; j <= (i < n ? i+2 : n+1); j++)
           printf("%15.5f", p[j]);
           if (i < n-1) printf("\n                     ");}
   printf("\n");
	printf(" \n  I    X       Y       YFIT  R=Y-YFIT\n \n");
	/*++ Code for ~.C. is inactive
	 *      print
	 *     *'(/'' NDEG ='',I2,10X,''RNORM ='',F8.4,10X,''SIGFAC ='',F8.4//
	 *     *'' P(1),P(2) ='',9X,2F15.5//'' P(3),...,P(NDEG+3) ='',3F15.5/
	 *     *(21X,3F15.5))', NDEG,RNORM(1),SIGFAC,(P(I),I=1,N+2)
	 *++ End */
	for (i = 1; i <= MDATA; i++)
	{
		yfit = dmpval( p, ndeg, X[i] );
		r = Y[i] - yfit;
		printf(" %2ld%6.0f%9.3f%9.3f%10.3f\n", i, X[i], Y[i], yfit, r);
	}
	exit(0);
 
} /* end of function */