// test_Big_simeq.java  tough simultaneous equations to solv accurately 

import java.math.BigDecimal;

public class test_Big_simeq
{
  final int nbit = 300; // bits precision 

  public test_Big_simeq()
  {
    int n;
    BigDecimal err = new BigDecimal(0.0);

    System.out.println("test_Big_simeq.java on tough equations to solve");
    n = 1;
    while(n < 128)
    {
      n=n*2;
      System.out.println("initializing big matrix, n="+n+", n*n="+(n*n));
      BigDecimal AA[][] = new BigDecimal[n][n];
      BigDecimal Y[] = new BigDecimal[n];
      BigDecimal X[] = new BigDecimal[n];
      BigDecimal one = new BigDecimal(1.0);
      BigDecimal zero = new BigDecimal(0.0);
      BigDecimal tmp = new BigDecimal(0.0);
      BigDecimal sol = new BigDecimal(0.0);

      one = one.setScale(nbit, BigDecimal.ROUND_DOWN);  
      zero = zero.setScale(nbit, BigDecimal.ROUND_DOWN);
      tmp = tmp.setScale(nbit, BigDecimal.ROUND_DOWN);
      for(int i=0; i<n; i++)
      {
        for(int j=0; j<n; j++)
        {
          AA[i][j] = one.divide(new BigDecimal((i+1)+(j+1)),
                                BigDecimal.ROUND_DOWN);
          AA[i][j] = AA[i][j].setScale(nbit, BigDecimal.ROUND_DOWN);
        }
      }
      for(int i=0; i<n; i++)
      {
	tmp = zero;
        for(int j=0; j<n; j++)
        {
	  sol = new BigDecimal(j+1);
          tmp = tmp.add(sol.multiply(AA[i][j]));
        }
        Y[i] = tmp;
      }

      new Big_simeq(AA, Y, X);
      if(n<=8)
      {
        System.out.println("simeq");
        for(int i=0; i<n; i++)
        {
            System.out.println("X["+i+"]="+X[i].toString());
        }
      } 
      err = zero;
      for(int i=0; i<n; i++)
      {
	  err = err.add((X[i].subtract(new BigDecimal(i+1))).abs());
      }
      err = err.setScale(nbit, BigDecimal.ROUND_DOWN);
      System.out.println("sum of error="+err);
      err = (err.divide(new BigDecimal(n), BigDecimal.ROUND_DOWN));
      err = err.setScale(nbit, BigDecimal.ROUND_DOWN);
      System.out.println("avg error="+err);
    }
    System.out.println("test double to Big to double");
    double A[][] = new double[3][4];
    double U[] = new double[3];
    A[0][0] = 1.0; A[0][1] = 1.0; A[0][2] = 1.0; A[0][3] = 6.0;
    A[1][0] = 1.0; A[1][1] = 2.0; A[1][2] = 3.0; A[1][3] = 14.0;
    A[2][0] = 2.0; A[2][1] = 1.0; A[2][2] = 3.0; A[2][3] = 13.0;
    new Big_simeq(3, A, U);
    System.out.println("U012="+U[0]+","+U[1]+","+U[2]);

    System.out.println("test Hilbert double to Big to double");
    n = 10;
    System.out.println("n="+n+", smallest="+(1.0/(double)(n*(n-1))));
    double B[][] = new double[n][n+1];
    double UB[] = new double[n];
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
	  B[i][j] = 1.0/(double)(i+j+1);
      B[i][n] = 0.0;
    }
    for(int i=0; i<n; i++)
      for(int j=0; j<n; j++)
	  B[i][n] += B[i][j]*(double)j ;
    System.out.println("B[0]0123="+B[0][0]+", "+B[0][1]+", "+B[0][2]+", "+B[0][3]);
    System.out.println("B[1]0123="+B[1][0]+", "+B[1][1]+", "+B[1][2]+", "+B[1][3]);
    System.out.println("B[2]0123="+B[2][0]+", "+B[2][1]+", "+B[2][2]+", "+B[2][3]);
    System.out.println("B[3]0123="+B[3][0]+", "+B[3][1]+", "+B[3][2]+", "+B[3][3]);

    new Big_simeq(n, B, UB);
    System.out.println("UB012="+UB[0]+", "+UB[1]+", "+UB[2]);
    System.out.println("UBn-1="+UB[n-1]+", "+UB[n-2]+", "+UB[n-3]);

    n = 16;
    System.out.println("n="+n+", smallest="+(1.0/(double)(n*(n-1))));
    double Bb[][] = new double[n][n+1];
    double UBb[] = new double[n];
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
	  Bb[i][j] = 1.0/(double)(i+j+1);
      Bb[i][n] = 0.0;
    }
    for(int i=0; i<n; i++)
      for(int j=0; j<n; j++)
	  Bb[i][n] += Bb[i][j]*(double)j ;
    System.out.println("Bb[0]0123="+Bb[0][0]+", "+Bb[0][1]+", "+Bb[0][2]+", "+Bb[0][3]);
    System.out.println("Bb[1]0123="+Bb[1][0]+", "+Bb[1][1]+", "+Bb[1][2]+", "+Bb[1][3]);
    System.out.println("Bb[2]0123="+Bb[2][0]+", "+Bb[2][1]+", "+Bb[2][2]+", "+Bb[2][3]);
    System.out.println("Bb[3]0123="+Bb[3][0]+", "+Bb[3][1]+", "+Bb[3][2]+", "+Bb[3][3]);

    new Big_simeq(n, Bb, UBb);
    System.out.println("UBb012="+UBb[0]+", "+UBb[1]+", "+UBb[2]);
    System.out.println("UBbn-1="+UBb[n-1]+", "+UBb[n-2]+", "+UBb[n-3]);

    System.out.println("finished test_Big_simeq.java");
  } // end test_Big_simeq

  public static void main (String[] args)
  {
    new test_Big_simeq();
  } // end main of test_Big_simeq.java
} // end class test_Big_simeq