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

extern dbl poly_1(dbl x)
{
	return x*x*x - 3.0*x*x + 7.0*x - 0.25;
}

extern dbl int_poly_1(dbl x) /* integral of poly_1(x) */
{
	return x*x*x*x/4.00 - 3.0*x*x*x/3.00 + 7.0*x*x/2.00 - 0.25*x;
}

extern dbl poly (dbl x, int n, dbl *a)
/*
  poly(x) = x^n + a[n-1]x^(n-1) + a[n-2]x^(n-2) + ... + a[1]^x + a[0]
  */
{
	int	i, k;
	dbl	mul;
	
	mul = 1;
	for (i=0, k=n-1; i<n; i++, k--) {
		mul = a[k] + mul*x;
	}
	return(mul);
}

void test1()
{
	/* integral of cos(x) */
	fprintf(stderr, "comp = %l30.26f\n", TrapezoidIntegral(cos, -1.00, 1.00, 100));
	fprintf(stderr, "comp = %l30.26f\n", GaussIntegral(cos, -1.00, 1.00));
	fprintf(stderr, "real = %l30.26f\n", sin(1.00)-sin(-1.00));
	fprintf(stderr, "\n");
	
	/* integral of e(x) */
	fprintf(stderr, "comp = %l30.26f\n", TrapezoidIntegral(exp, -2.00, 5.00, 100));
	fprintf(stderr, "comp = %l30.26f\n", GaussIntegral(exp, -2.00, 5.00));
	fprintf(stderr, "real = %l30.26f\n", exp(5.00)-exp(-2.00));
	fprintf(stderr, "\n");
	
	/* integral of x^3 - 3.0 x^2 + 7.0 x - 0.25 */
	fprintf(stderr, "comp = %lf\n", TrapezoidIntegral(poly_1, -1.00, 10.00, 10000));
	fprintf(stderr, "comp = %lf\n", GaussIntegral(poly_1, -1.00, 10.00));
	fprintf(stderr, "real = %lf\n", int_poly_1(10.00)-int_poly_1(-1.00));
	fprintf(stderr, "\n");
}

void test2()
{
	int dim, i;
	dbl a[BUFSIZ], x;
	
	/* integral of x^3 - 3.0 x^2 + 7.0 x - 0.25 */
	dim = 3;
	a[0]=(-0.25);
	a[1]=7.0;
	a[2]=(-3.0);
	fprintf(stderr, "comp = %l30.26f\n", GaussIntegralParams(poly, -1.00, 10.00, dim, a));
	fprintf(stderr, "\n");
}

main()
{
	/* GaussIntegralInitialize("gauss10"); */
	test1();
	/*test2();*/
}