#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define USE_GMP
#define USE_MPFR
#include "bnc.h"

#define DIM 512

void get_mpfproblem(MPFMatrix a, MPFVector b, MPFVector ans, long dim)
{
	long int i, j, k;
	mpf_t tmp, sqr2;

	mpf_init2(tmp, prec_mpfvector(ans));
	mpf_init2(sqr2, prec_mpfvector(ans));

	mpf_set_ui(sqr2, 2UL); mpf_sqrt(sqr2, sqr2);

	/* Frank Matrix */
	for(i = 0; i < dim; i++)
	{
		for(j = 0; j < dim; j++)
		{
			if(i < j)
			{
				mpf_set_si(tmp, dim - j);
				set_mpfmatrix_ij(a, i, j, tmp);
			}
			else
			{
				mpf_set_si(tmp, dim - i);
				set_mpfmatrix_ij(a, i, j, tmp);
			}
			mpf_mul(tmp, tmp, sqr2);
			set_mpfmatrix_ij(a, i, j, tmp);
		}
	}

	/* Answer */
	for(i = 0; i < dim; i++)
	{
		mpf_set_si(tmp, i);
		set_mpfvector_i(ans, i, tmp);
	}

	/* Make constant vector */
	mul_mpfmatrix_mpfvec(b, a, ans);

	mpf_clear(tmp);
	mpf_clear(sqr2);
}

int main(int argc, char *argv[])
{
	double start, dtime, startwtime[2], endwtime[2];

	MPFMatrix mpfa;
	MPFVector mpfb, mpfx, mpfans;
	mpf_t reps, aeps;
	long int itimes_mpf;
	double mpftime;

#define MPF_PREC 128

	/* initialize */
	mpf_init(reps);
	mpf_init(aeps);

	mpfa  = init_mpfmatrix(DIM, DIM);
	mpfb  = init_mpfvector(DIM);
	mpfx  = init_mpfvector(DIM);
	mpfans  = init_mpfvector(DIM);

	/* get problem */
	get_mpfproblem(mpfa, mpfb, mpfans, DIM);

	/* run MPFFCG */
	mpf_set_d(reps, 1.0e-20);
	mpf_set_d(aeps, 1.0e-50);

	start = get_secv();
	itimes_mpf = MPFCG(mpfx, mpfa, mpfb, reps, aeps, DIM * 5);
	mpftime = get_secv() - start;

	print_mpfvector(mpfx);

	free_mpfmatrix(mpfa);
	free_mpfvector(mpfb);
	free_mpfvector(mpfx);
	free_mpfvector(mpfans);

	/* end */
	mpf_clear(reps); mpf_clear(aeps);

	/* print itimes */
	printf("mpf_t(%d)  : %ld(%f)\n", MPF_PREC, itimes_mpf, mpftime);

	return EXIT_SUCCESS;
}