/* test_mpf_simeq.c */
/* gcc -o test_mpf_simeq -I. mpf_simeq.c test_mpf_simeq.c -L. -lgmp */
#include <stdio.h>
#include "gmp.h"
#include "mpf_simeq.h"
#define n 4
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

int main()
{
  mpf_set_default_prec(50*3.32); /* 50 digits */
  mpf_t one, zero, sum, temp;
  mpf_t A[n*n];
  mpf_t X[n];
  mpf_t Y[n];
  int i, j;

  printf("test_mpf_simeq, n=%d  \n", n);
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++) mpf_init(A[i*n+j]);
    mpf_init(X[i]);
    mpf_init(Y[i]);
  }
  mpf_init(one); /* must init every variable */
  mpf_init(sum);
  mpf_init(temp);
  mpf_init(zero);

  mpf_set_str(one, "1.0", 10);
  mpf_set_si(zero, 0);
  
  for(i=0; i<n; i++)
  {
    mpf_set_str(sum, "0", 10);  
    for(j=0; j<n; j++)
    {
      mpf_set_si(temp, max(i+1,j+1));
      mpf_div(temp, one, temp);
      mpf_set(A[i*n+j], temp);
      gmp_printf("A[%d,%d] = %50.47Ff \n", i, j, A[i*n+j]);
      mpf_add(sum, sum, temp);
    }
    mpf_set(Y[i], sum);
    gmp_printf("Y[%d] =   %50.47Ff \n", i, sum);
  }

  mpf_simeq(n, A, Y, X);
  for(i=0; i<n; i++)
  {
    gmp_printf("X[%d]=%50.47Ff \n", i, X[i]);
  }
  for(i=0; i<n; i++)
  {
    mpf_sub(temp, zero, Y[i]);
    mpf_set(sum, temp);
    for(j=0; j<n; j++)
    {
      mpf_mul(temp, A[i*n+j], X[j]);
      mpf_add(sum, sum, temp);
    }
    gmp_printf("row %d, sum err = %Fg \n", i, sum);
  }

  printf("\n test_mpf_simeq finished \n");
  return 0;
} /* end main of test_mpf_simeq.c */