import java.util.*;
import java.lang.*;
import java.text.*;

public class Util{
    private static final double Log10e = 0.4342944819032518276511289;
//  private static final double Log10e = 0.434294481903;
    boolean getout=false;
    private double V1[] = new double[20];
    private int nswaps;
    public boolean GAUSS_SINGULAR;
    public double GAUSS_DETERM;

    Util(){}
    public double sign(double number){
        double zero=0.0e0, one=1.0e0;
        if(number >= zero){
            return(one);
        }else{
            return(-one);
        }
    }
    public int Ilog10(double number){
        double val;
        val = Math.abs(number);
        return( (int)log10(val) );
    }
    public double mod(double f1, double f2){
        return(Math.IEEEremainder(f1, f2));
    }
    public double log10(double number){
        return(Log10e * Math.log(number));
    }
    public void AX_Bgauss(int DIMA, double a[][], double b[], 
                          double x[])
    {
        double determ, one=1.0e0;
        boolean stat; int n, k;

        double y[][] = new double[DIMA][1];
        for(n=0 ; n<DIMA; n++)y[n][0] = b[n];
        stat=trianmat(DIMA, a, y, 1);
        GAUSS_SINGULAR = stat;
        if(GAUSS_SINGULAR)return;
        for(determ=one, k=0 ; k<DIMA; k++){
	      x[k] = y[k][0];
	      determ *= a[k][k];
        }
        if(nswaps%2 !=  0)determ = -determ;
        GAUSS_DETERM = determ;
        return;
    }

/* invgauss: Inverser la matrice A par la methode de Gauss.
             L'inverse de A correspond a la matrice X.
             Le determinant de A est retourne.  */
    public void invgauss(int DIMA, double a[][], double x[][])
    {
        double determ, zero=0.0e0, one=1.0e0;
        boolean stat; int n, k;

        for(n=0 ; n<DIMA; n++){
	      for(k=0 ; k<=n; k++){
	          x[n][k] = zero;
	          x[k][n] = x[n][k];
	          if(k == n)x[k][n] = one ;
	      }
        }
        stat=trianmat(DIMA, a, x, DIMA);
        GAUSS_SINGULAR = stat;
        if(GAUSS_SINGULAR)return;
        for(determ=one, k=0; k<DIMA; k++)determ *= a[k][k];
        if(nswaps%2 !=  0)determ = -determ;
        GAUSS_DETERM = determ;
        return;
    }

    /* trianmat: Triangulariser le haut de la matrice A en 
       operant conjointement sur la matrice a NCOLB colonnes B */
    public boolean trianmat(int DIMA, double a[][], double b[][],
	           int NCOLB)
    {
        int m, n, k, l, i, j, index;   
        double delta = 1.0e-20;
        double maxelem, rapport;

        double xcol[] = new double[DIMA];
        nswaps = 0;
        for(m=0; m<DIMA; m++){
	      maxelem = Math.abs(a[m][m]);
	      index = m;
	      for(n=m+1; n<DIMA; n++){
	          if( maxelem < Math.abs(a[n][m]) ){
		        maxelem = Math.abs(a[n][m]);
		        index = n;
	          }
            }
            if(maxelem < delta)return(true);
	      if(index > m){
	          nswaps++; 
	          swaplines(a, index, m, DIMA);
	          swaplines(b, index, m, NCOLB);
            }
	      for(n=m+1; n<DIMA; n++){
	          rapport = a[n][m] / a[m][m];
	          for(k=m; k<DIMA; k++){
		        a[n][k] -= a[m][k] * rapport;
	          }
	          for(k=0; k<NCOLB; k++){
		        b[n][k] -= b[m][k] * rapport;
	          }
            }
        }

    /* La matrice A est sous forme triangulaire et la matrice B a subi
       les memes transformations que A. Pour chaque colonne i de B,
       le vecteur x(i) de Ax(i)=b(i) est determine. */
        for(i=0; i<NCOLB; i++){
	      for(m=DIMA-1; m>=0; m--){
	          xcol[m] = b[m][i];
	          for(n=m+1; n<DIMA; n++)xcol[m] -= a[m][n] * xcol[n];
	          xcol[m] /= a[m][m];
	      }
	      for(m=0; m<DIMA; m++)b[m][i] = xcol[m];
        }
        return(false);
    }
    /* swaplines: Echanger deux lignes d'un tableau */
    public void swaplines(double a[][], int i, int j, int n)
    {
        int m;
        double temp;

        for(m=0; m<n ; m++){
	      temp = a[i][m];
	      a[i][m] = a[j][m];
	      a[j][m] = temp;
        }
        return;
    }
//  Compute A' = A.B
    public void Matprod1(double A[][], double B[][], int n)
    {
        int i, j, k; double zero=0.0e0;

        for(i=0; i<n; i++){
	      for(j=0; j<n; j++){
	          for(V1[j]=zero, k=0; k<n; k++){
		        V1[j] += A[i][k] * B[k][j];
	          }
	      }
	      for(j=0; j<n; j++)A[i][j] = V1[j];
        }
        return;
    }
    public void transpose(double a[][], int n)
    {
        int i, j; double temp;
    
        for(i=0; i<n; i++){
	      for(j=0; j<=i; j++){
	          temp = a[i][j];
	          a[i][j] = a[j][i];
	          a[j][i] = temp;
	      }
        }
        return;
    }
    public void zero_mat(double A[][], int n, int m)
    {
        int i, j; double zero=0.0e0;

        for(i=0; i<n; i++){
	      for(j=0; j<m; j++){
	          A[i][j] = zero;
	      }
        }
        return;
    }
    public void sortAscending(int n, double V[])
    {
        int i, j; double temp;
        for(i=1; i<n; i++){
            for(j=0; j<n-1; j++){
                if(V[j] > V[i]){
                    temp = V[i];              
                    V[i] = V[j];
                    V[j] = temp;
                }
            }
        }
    }
    public void sortDescending(int n, double V[])
    {
        int i, j; double temp;
        for(i=1; i<n; i++){
            for(j=0; j<n-1; j++){
                if(V[j] < V[i]){
                    temp = V[i];              
                    V[i] = V[j];
                    V[j] = temp;
                }
            }
        }
    }
}