//******************************************************************************
// Lab.java:	Applet
// CS 20c lab 2 
// Brian Frazier
// frazier@cs.caltech.edu
// This applet implements both fast polynomial multiplication and
// LU decomposition of matrices
//
//******************************************************************************
import java.applet.*;
import java.awt.*;
import java.util.*;
import java.lang.Math;


/*
 *
 * Complex
 *
 */
//The class Complex encapsulates all the necessary 
//functions of a Complex variable

class Complex 
{
        public double Re;        //Real part of complex variable
        public double Im;        //Imaginery part of the complex variable

        // Complex() method constructs a complex and initializes its values 
        public Complex()
        {
                Re=0.0;
                Im=0.0;
        }

        // add method adds two complexs together
        public void add(Complex a, Complex b)
        {
                Re = (a.Re + b.Re);
                Im = (b.Im + a.Im);
        }

        // subtract method subtracts two complexs
        public void subtract(Complex a,Complex b)
        {
                Re = a.Re - b.Re;
                Im = a.Im - b.Im;
        }

        // multiply method multiplies two complexs together
        public void multiply(Complex a,Complex b)
        {
                Re = (a.Re * b.Re) - (a.Im * b.Im);
                Im = (a.Re * b.Im) + (b.Re * a.Im);

        }

        // divide method creates an inverse of the functions in the denominator
        // then multiplies the numerator by the inverse of the denominator
        public void divide(Complex a,Complex b)
        {
                Complex c = new Complex();
                Re = (b.Re / ((b.Re * b.Re)+(b.Im * b.Im)));
                Im = (((-1.0)*b.Im) / ((b.Re * b.Re)+(b.Im * b.Im)));
                c.Re = Re;
                c.Im = Im;
                multiply(a,c);
        }                                                              
        
		// rootun method returns the root of unity of Omega(k,n)
		public void rootun(double k,double n)
		{
				Re = Math.cos((2*k*Math.PI)/n);
				Im = Math.sin((2*k*Math.PI)/n);
		}
}


//
//
// Matrix
//
//
// Matrix is a class that stores a matrix in a 2 dimensional
// array with the first index representing the rows and the 
// second index representing the columns
class Matrix 
{
	int rows,columns,flag;
	double mat[][];

	public Matrix(int r, int c)
	{
		rows = r;
		columns = c;
		flag=0;
		mat = new double [rows][columns];
	}

}


//
//
// MatrixOp
//
//
// MatrixOp is a class that contains functions to operate
// on matrices.  All the function take input of the class Matrix
// and return matrices of the type Matrix
class MatrixOp
{

	// Matrix add is a method that simply add two matrices
	// a and b together by adding each component a(x,y) to 
	// b(x,y).  A matrix of the proper size is returned
	// with the result
	public Matrix add(Matrix a, Matrix b)
	{
		Matrix c = new Matrix(a.rows,a.columns);
		for(int x=0;x<c.rows;x++)
		{   for(int y=0;y<c.columns;y++)
			{
			// add individual elements together
			c.mat[x][y] = a.mat[x][y] + b.mat[x][y];
		}}
		return c;
	}


	// scalmul is a method to multiply a given matrix
	// by a scaler.  This is done through the multiplication
	// of each element by the respective scaler.  A matrix
	// of the same size and the input matrix is returned
	public Matrix scalmul(Matrix a, double scale)
	{
		Matrix c = new Matrix(a.rows,a.columns);
		for(int x=0;x<a.rows;x++)
		{
			for(int y=0;y<a.columns;y++)
			{
				// multiply individual element by scaler
				c.mat[x][y] = scale *  a.mat[x][y];
			}
		}
		return c;
	
	}

	// sub is a method that subtract Matrix b from Matrix a
	// this is done through the subtraction of individual elements
	// of the respective matrices
	public Matrix sub(Matrix a, Matrix b)
	{
		Matrix c = new Matrix(a.rows,a.columns);
		for(int x=0;x<c.rows;x++)
		{
			for(int y=0;y<c.columns;y++)
			{
			// subtract element in b(x,y) from a(x,y)
			c.mat[x][y] = a.mat[x][y] - b.mat[x][y];
			}
		}
		return c;
	}

	
	// mul is a method to multiply two matrices a and b
	// together. It does this by multiplying each
	// element in a row in "a" with the respective element of 
	// each column in "b" and then summing the product of a 
	// row times a column
	public Matrix mul(Matrix a, Matrix b)
	{
		Matrix c = new Matrix(a.rows,b.columns);
		double temp = 0.0;
		for(int z=0;z<c.rows;z++)
		{
			for(int x=0;x<c.columns;x++)
			{
				// sum up products of individual
				// row times a column
				for(int y=0;y<a.columns;y++)
				{ temp= temp + a.mat[z][y] * b.mat[y][x]; }
				// place row x column sum in respective
				// resultant matrix position
				c.mat[z][x] = temp;
				temp = 0.0;
			}
		}
		return c;
	}


	// LU performs LU decomposition of a matrix "a" passed
	// to it. It returns a matrix twice the size of the 
	// original passed to it with L in the left side
	// and U in the right side.  LU assumes that a is square
	public Matrix LU(Matrix a)
	{
		// create a new Matrix to store L and U in
		Matrix b = new Matrix(a.rows,2*(a.rows));
		// create a matrix to get the recursive LU results
		Matrix ret = new Matrix((a.rows)-1,2*((a.rows)-1));
		Matrix u,v,r,r2,t;
		
		double alpha = 0.0;

		// If the matrix is a single element
		// return L = 1 and U = A
		if(a.rows==1)
		{
			b.mat[0][0] = 1.0;
			b.mat[0][1] = a.mat[0][0];
			return b;
		}
		// If an alpha of 0 encountered
		// return error
		if(a.mat[0][0]==0)
		{
			b.flag=1;
			return b;
		}
		alpha = 1/(a.mat[0][0]);
		// create u and v to be n-1 in size
		u = new Matrix((a.rows)-1,(a.columns)-1);
		v = new Matrix((a.rows)-1,(a.columns)-1);
		// copy the first column values in a to u
		for(int y=0;y<u.columns;y++)
		{	
			u.mat[y][0] = a.mat[y+1][0];
		}
		// copy the first row values in a to v
		for(int x=0;x<v.rows;x++)
		{
			v.mat[0][x] = a.mat[0][x+1];
		}
		r = new Matrix((a.rows)-1,(a.columns)-1);
		r2 = new Matrix((a.rows)-1,(a.columns)-1);
		t = new Matrix((a.rows)-1,(a.columns)-1);
		// r = u * v
		r = mul(u,v);
		// r2 = (alpha*u*v)
		r2 = scalmul(r,alpha);
		for(int x=0;x<v.rows;x++)
		{
			for(int y=0;y<v.columns;y++)
			{
				r.mat[x][y] = a.mat[x+1][y+1];
			}
		}
		// resurse on the matrix a(n-1)
		ret = LU(sub(r,r2));
		b.flag=ret.flag;
		
		// copy the components of L(n-1) and U(n-1)
		// returned in ret to the matrix to return
		b.mat[0][0]=1;
		b.mat[0][(ret.columns/2)+1]=1/alpha;
		for(int x=0;x<ret.rows;x++)
		{
			for(int y=0;y<ret.rows;y++)
			{
				b.mat[x+1][y+1] = ret.mat[x][y];
				b.mat[x+1][(y+1)+(ret.columns/2)+1] = ret.mat[x][y+(ret.columns/2)];
			}
		}
		// let L(2 to n) elements = u*alpha;
		t = scalmul(u,alpha);

		for(int x=0;x<ret.rows;x++)
		{
			//zero out the to components of L from 2 to n
			b.mat[0][x+1]=0.0;
			//let U top row elements from 2 to n be elements of v
			b.mat[0][(x+1)+(ret.columns/2)+1]=v.mat[0][x];
			b.mat[x+1][0]=t.mat[x][0];
			// let first columns of u from 2 to n be 0
			b.mat[x+1][(ret.columns/2)+1]=0.0;
		}
		// return LU
		return b;
	}


	// Trid performs tridiagonalization on a symmetric matrix.  
	// it assumes that a is symmetric
	public Matrix Trid(Matrix a)
	{
		double scaler,h,k, au, f;
		int l;
		Matrix ret = new Matrix(2,a.rows);
		
		for(int i=(a.rows-1);i>=1;i--)
		{
			l = i - 1;
			scaler=0.0;
			h=0.0;
			if(l > 0)
			{
				for(int x=0;x<=l;x++)
				{
					scaler = scaler + Math.abs(a.mat[i][x]);
				}
				if(scaler == 0.0)
					ret.mat[1][i] = a.mat[i][l];
				else
				{
					for(int y=0;y<=l;y++)
					{
						a.mat[i][y] = a.mat[i][y] / scaler;
						h = h + a.mat[i][y]*a.mat[i][y];
					}
					f=a.mat[i][l];
					if(f >= 0.0)
						au = Math.sqrt(h);
					else
						au = Math.sqrt(h * -1);
					ret.mat[1][i]=scaler * au;
					h = h - f*au;
					a.mat[i][l]=f-au;
					f=0.0;
					for(int z=0;z<=l;z++)
					{
						au=0.0;
						for(int c=0;c<=z;c++)
							au = au + a.mat[z][c]*a.mat[i][c];
						for(int c=z+1;c<=l;c++)
							au = au + a.mat[c][z]*a.mat[i][c];
						ret.mat[1][z]=au/h;
						f = f + ret.mat[1][z]*a.mat[i][z];
					}
					k=f/(h+h);
					for(int z=0;z<=l;z++)
					{
						f=a.mat[i][z];
						ret.mat[1][z]=au=ret.mat[1][z] - k*f;
						for(int c=0;c<=z;c++)
							a.mat[z][c] = a.mat[z][c] - (f * ret.mat[1][c] + au*a.mat[i][c]);
					}
				}
			}
			else
				ret.mat[1][i]=a.mat[i][l];
			ret.mat[0][i]=h;
		}
		ret.mat[1][0]=0.0;
		for(int i=0;i<a.rows;i++)
			ret.mat[0][i]=a.mat[i][i];
		return ret;
	}


	// FFT is a recursive FFT method to transform a 
	// matrix a into its fourier components
	public Matrix FFT(Matrix a)
	{
		Matrix y = new Matrix(a.rows,a.columns);
		Matrix b,c,p,q;
		Complex roo = new Complex();
		Complex roo1 = new Complex();
		Complex p1 = new Complex();
		Complex q1 = new Complex();
		Complex y1 = new Complex();
		Complex y2 = new Complex();

		// If a is one element, return it
		if(a.columns==1)
		{
			y.mat[0][0]=a.mat[0][0];
			y.mat[1][0]=a.mat[1][0];
			return y;
		}

		// divide "a" up into two vectors half the size
		// b represents the odd coefficients and 
		// c represents the even coefficients
		b = new Matrix(a.rows,(a.columns)/2);
		c = new Matrix(a.rows,(a.columns)/2);
		for(int x=0; x<(a.columns/2);x++)
		{
			b.mat[0][x]=a.mat[0][2*x];
			b.mat[1][x]=a.mat[1][2*x];
			c.mat[0][x]=a.mat[0][2*x+1];
			c.mat[1][x]=a.mat[1][2*x+1];
		}

		// recurse on the even and odd elements
		p = FFT(b);
		q = FFT(c);

		// compute first half of y to return
		// y = p + w*q (where w is ith root of unity)
		for(int i=0;i<=((a.columns)/2)-1;i++)
		{
			roo.rootun(i,a.columns);
			// convert vector to complex class
			p1.Re=p.mat[0][i];
			p1.Im=p.mat[1][i];
			q1.Re=q.mat[0][i];
			q1.Im=q.mat[1][i];

			y1.multiply(roo,q1);
			y2.add(y1,p1);
			y.mat[0][i]=y2.Re;
			y.mat[1][i]=y2.Im;			
		}

		// compute second half of y to return
		for(int i=(a.columns/2);i<a.columns;i++)
		{
			// generate root of unity of i
			roo.rootun(i,a.columns);
			p1.Re=p.mat[0][i-((a.columns)/2)];
			p1.Im=p.mat[1][i-((a.columns)/2)];
			q1.Re=q.mat[0][i-((a.columns)/2)];
			q1.Im=q.mat[1][i-((a.columns)/2)];
			
			y1.multiply(roo,q1);
			y2.add(y1,p1);
			
			y.mat[0][i]=y2.Re;
			y.mat[1][i]=y2.Im;
		}
		return y;
	}

	// invFFT performs the inverse fast fourier
	// transform on a matrix a.  It assumes that
	// a has number of elements on order of 2^n
	public Matrix invFFT(Matrix a)
	{
		Matrix y = new Matrix(a.rows,a.columns);
		Matrix b,c,p,q;
		Complex roo = new Complex();
		Complex roo1 = new Complex();
		Complex p1 = new Complex();
		Complex q1 = new Complex();
		Complex y1 = new Complex();
		Complex y2 = new Complex();

		// if "a" is a single element, return it
		if(a.columns==1)
		{
			y.mat[0][0]=a.mat[0][0];
			y.mat[1][0]=a.mat[1][0];
			return y;
		}

		// divide up elements of "a" into even and odd
		// b represents odd elements of a
		// c represents even elements of a
		b = new Matrix(a.rows,(a.columns)/2);
		c = new Matrix(a.rows,(a.columns)/2);
		for(int x=0; x<(a.columns/2);x++)
		{
			b.mat[0][x]=a.mat[0][2*x];
			b.mat[1][x]=a.mat[1][2*x];
			c.mat[0][x]=a.mat[0][2*x+1];
			c.mat[1][x]=a.mat[1][2*x+1];
		}

		// recursively perform inverse FFT on 
		// even and odd elements of "a"
		p = invFFT(b);
		q = invFFT(c);

		// compute first half of resultant
		// vector y
		// y = p + (1/w)*q (w is root of unity of i)
		for(int i=0;i<=((a.columns)/2)-1;i++)
		{
			// compute root of unity of ith element
			roo.rootun(i,a.columns);
			// get ith element of p and q
			p1.Re=p.mat[0][i];
			p1.Im=p.mat[1][i];
			q1.Re=q.mat[0][i];
			q1.Im=q.mat[1][i];

			y1.divide(q1,roo);
			y2.add(y1,p1);
			y.mat[0][i]=y2.Re;
			y.mat[1][i]=y2.Im;			
		}
		// compute second half of resultant polynomial y
		for(int i=(a.columns/2);i<a.columns;i++)
		{
			// generate root of unity of ith element
			roo.rootun(i,a.columns);
			// get i-(n/2) elements of p and q
			p1.Re=p.mat[0][i-((a.columns)/2)];
			p1.Im=p.mat[1][i-((a.columns)/2)];
			q1.Re=q.mat[0][i-((a.columns)/2)];
			q1.Im=q.mat[1][i-((a.columns)/2)];
			
			y1.divide(q1,roo);
			y2.add(y1,p1);
			
			y.mat[0][i]=y2.Re;
			y.mat[1][i]=y2.Im;
		}

		return y;
	}

}

//==============================================================================
// Main Class for applet Lab
//
//==============================================================================
public class Lab extends Applet
{
	Matrix a,at,b,bt,c,d; // Matrices used in the LU and FFT
	String in,in2;		  // Input strings
	Button Calculate;	  // LU decomposition button
	Button FFT;			  // FFT polynomial mult. button
	TextField input;	  // first input field
	TextField input2;	  // second input field
	TextField size;		  // size of LU matrix input field
	int flag;			  // global flag

	// Lab Class Constructor
	//--------------------------------------------------------------------------
	public Lab()
	{
	

	}


	// The init() method is called by the AWT when an applet is first loaded or
	// reloaded.  Override this method to perform whatever initialization your
	// applet needs, such as initializing data structures, loading images or
	// fonts, creating frame windows, setting the layout manager, or adding UI
	// components.
    //--------------------------------------------------------------------------
	public void init()
	{
        // generate applet layout
		setLayout(new FlowLayout());
		input = new TextField(30);
		input2 = new TextField(30);
		size = new TextField(3);
		Label sizelabel = new Label("Matrix size: ");
		Label inputlabel = new Label("Input Matrix: ");
		Label input2label = new Label("2nd Input Matrix: ");
		Calculate = new Button("LU");
		FFT = new Button("FFT");
		inputlabel.LEFT;
		sizelabel.LEFT;

		// size request
		add(sizelabel);
		add(size);
		// add first input field
		add(inputlabel);
		add(input);
		// calculation buttons
		add(Calculate);
		add(FFT);
		input2label.RIGHT;
		// add second input field
		add(input2label);
		add(input2);

		flag=2; // set flag so only layout is drawn
    	resize(320, 240);
	}

	// destroy() is called when when you applet is terminating and being unloaded.
	//-------------------------------------------------------------------------
	public void destroy()
	{
		
	}

	// Paint Method
	//--------------------------------------------------------------------------
	public void paint(Graphics g)
	{
		
		// If the user choose polynomial multiplication
		// then output result
		if(flag==3)
		{
			Float intern;
			int x =	0;
		
			// output information strings
			g.drawString("Polynomial result of \"fast\" polynomial multiplication by FFT",10,70);
			g.drawString("Results are coefficients in increasing order",10,100);
			g.drawString("from left to right, and from top to bottom",10, 115);
			// print out polynomial coefficients in increasing order
			for(int y=0;y<c.columns;y++)
			{ 
				intern = new Float(c.mat[0][y]);
				if(((y*80)-(x*400)) == 400)
					x = x+1;
				g.drawString(intern.toString(), 10 + (y*80) - (x*400),150 + x*15); 
			}
		}
		
		// If user choose LU Decomposition
		// then output the L and U matrices
		if(flag < 2)
		{
			// If no errors occured in LU decomp
			// then output the matrices
			if(flag == 0 && c.flag==0)
			{
				Float intern;
				g.drawString("Simple LU-Decomposition",10,70);
				g.drawString("L",10 + (c.rows - 1)*45,90);
			
				// output L by cycling through rows
				for(int x=0; x < c.rows; x++)
				{
					for(int y=0; y < c.rows; y++)
					{	intern = new Float(c.mat[x][y]);
						g.drawString(intern.toString(), 10 + y*90, 110 + 20*x);
					}	
				}	

				// output U by cycling through rows
				g.drawString("U",10 + (c.rows - 1)*45,110 + 20*(c.rows-1)+50);
				for(int x=0; x < c.rows; x++)
				{
					for(int y=c.rows; y < c.columns; y++)
					{ 
						intern = new Float(c.mat[x][y]);
						g.drawString(intern.toString(), 10 + (y-c.rows)*90,110 + 20*(c.rows-1)+75+20*x);
					}	
				}
			}
			// if user did not enter enough elements for 
			// matrix desired, then return error
			else if(flag == 1)
				g.drawString("Invalid number of elements",10,60);
			//if alpha = 0 encountered in simple LU decomposition
			// return message
			else if(c.flag == 1)
				g.drawString("alpha = 0 encountered",10,60);
		}
		
	}
	//		The start() method is called when the page containing the applet
	// first appears on the screen. The AppletWizard's initial implementation
	// of this method starts execution of the applet's thread.
	//--------------------------------------------------------------------------
	public void start()
	{
		
	}
	
	//		The stop() method is called when the page containing the applet is
	// no longer on the screen. The AppletWizard's initial implementation of
	// this method stops execution of the applet's thread.
	//--------------------------------------------------------------------------
	public void stop()
	{
	}


	public boolean action(Event event, Object arg)
	{
		
		// LU decomposition request by user
		if(event.target == Calculate)
		{
			StringTokenizer tok;
			int temp,to;
			MatrixOp k = new MatrixOp();
		
			Graphics g = this.getGraphics();
			flag = 0;
			// get the size of the matrix desired
			in = this.size.getText();
			if(in.length() != 1)
			{
				flag=1;
				repaint();
				return true;
			}
			temp = (Integer.valueOf(in)).intValue();
			// get elements of matrix
			in = this.input.getText();
			tok = new StringTokenizer(in);
			// make sure there are enough elements
			if((temp*temp) != tok.countTokens())
			{
				flag=1;
				repaint();
				return true;
			}
			// if there are enough elements
			// place them in matrix
			else
			{
				a = new Matrix(temp,temp);
				for(int x=0;x<temp;x++)
				{
					for(int y=0;y<temp;y++)
					{
					a.mat[x][y]= (Double.valueOf(tok.nextToken())).doubleValue();
					}
				}
				// LU decompose the matrix
				c = k.LU(a);			
				repaint();
				return true;
			}
		}


		// user chose FFT polynomial multiplication
		if(event.target == FFT)
		{
			double tempr,tempr2;
			int aa = 2, bb=2;
			StringTokenizer tok,tok2;
			int temp,temp2;
			MatrixOp k = new MatrixOp();
			Complex p = new Complex();
			Complex q = new Complex();
			Complex c1 = new Complex();
			

			Graphics g = this.getGraphics();
			flag = 0;

			// get coefficients of first polynomial
			in = this.input.getText();
			tok = new StringTokenizer(in);
			temp = tok.countTokens();
			tempr= temp;
		
			// find the next highest order of 2 which the
			// coefficients can fit in
			if((tempr) != aa && (tempr) > 1)
			{
				while(aa<(tempr))
					aa = aa * 2;
			}

			// get coefficients of second polynomial
			in2 = this.input2.getText();
			tok2 = new StringTokenizer(in2);
			temp2 = tok2.countTokens();
			tempr2= temp2;
		
			// find the next highest order or 2 in which
			// the second polynomial can fit in
			if((tempr2) != bb && (tempr2) > 1)
			{
				while(bb<(tempr2))
					bb = bb * 2;
			}
			
			// create two vectors that are twice the
			// size of the largest of the two polynomial vectors
			if(aa < bb)
			{
				a = new Matrix(2,(2 * bb));
				b = new Matrix(2,(2 * bb));
				d = new Matrix(2,(2 * bb));
			}
			else
			{
				a = new Matrix(2,(2 * aa));
				b = new Matrix(2,(2 * aa));
				d = new Matrix(2,(2 * aa));
			}
			
			// copy elements to the vectors for processing
			for(int x=0;x<(temp);x++)
			{
			a.mat[0][x]= (Double.valueOf(tok.nextToken())).doubleValue();	
			a.mat[1][x]= 0.0;
			}

			for(int x=0;x<(temp2);x++)
			{
			b.mat[0][x]= (Double.valueOf(tok2.nextToken())).doubleValue();	
			b.mat[1][x]= 0.0;
			}

			//Fast Fourier Transform both polynomials
			at = k.FFT(a);
			bt = k.FFT(b);

			// multiple respective elements of each transform
			// together to form the final fourier orders
			for(int i=0;i<at.columns;i++)
			{	
			p.Re=at.mat[0][i];
			p.Im=at.mat[1][i];
			q.Re=bt.mat[0][i];
			q.Im=bt.mat[1][i];

			c1.multiply(p,q);
			d.mat[0][i]=c1.Re;
			d.mat[1][i]=c1.Im;			
			}

			// inverse fast fourier transform the final
			// set of orders to get resultant polynomial
			c = k.invFFT(d);

			// normalize the polynomial coefficients
			for(int i=0;i<c.columns;i++)
			{
				c.mat[0][i]=(c.mat[0][i]) / c.columns;
				c.mat[1][i]=(c.mat[1][i]) / c.columns;
			}
			for(int x=0;x<c.columns;x++)
			{
				if(c.mat[0][x]<1.0e-08 && c.mat[0][x]>-1.0e-08) c.mat[0][x]=0.0;
				if(c.mat[1][x]<1.0e-08 && c.mat[1][x]>-1.0e-08) c.mat[1][x]=0.0;
			}
			// set flag to draw the result
			flag=3;
			repaint();
			return true;
		}
		return false;
	}


	// TODO: Place additional applet code here

}