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

extern dbl f1(dbl x, dbl y)
{
	return x*y;
}

extern dbl f2(dbl x, dbl y)
{
	return x - y;
}

extern dbl f3(dbl x, dbl y, int n, dbl a[])
{
	return a[0]*x + a[1]*y;
}

void test1()
{
	/* integral of x*y */
	fprintf(stderr, "comp          = %30.26lf\n", GaussIntegral2D(f1, 0.0, 1.0, 0.0, 1.0));
	fprintf(stderr, "analytic sol. = %30.26lf\n", 1.0/4);
}

void test1marginal()
{
	dbl	x, y;
	
	/* int_{0.0}^{1.0} x*y dx = y  (y:constant) */
	for (y=0.00; y<=1.00; y+=0.20) {
		fprintf(stderr, "y = %lf   comp = %26.18f\n",
			y, GaussIntegral2Dmarginal1(f1, 0.0, 1.0, y));
	}
	
	/* int_{0.0}^{1.0} x*y dy = x  (x:constant) */
	for (x=0.00; x<=1.00; x+=0.20) {
		fprintf(stderr, "x = %lf   comp = %26.18f\n",
			x, GaussIntegral2Dmarginal2(f1, x, 0.0, 1.0));
	}
}

void test2()
{
	/* integral of x - y */
	fprintf(stderr, "comp          = %30.26lf\n", GaussIntegral2D(f2, 0.0, 1.0, 1.0, 2.0));
	fprintf(stderr, "analytic sol. = %30.26lf\n", -1.0);
}

void test3()
{
	static dbl	a[2] = {1.00, -1.00};
	
	/* integral of 1.00 x + (-1.00) y */
	vectorfprint(stderr, 2, a);
	fprintf(stderr, "comp          = %30.26lf\n",
		GaussIntegral2DParams(f3, 0.0, 1.0, 1.0, 2.0, 2, a));
}

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