// The inimmitable graph visualization explorer
public class Adventurer {
	public Graph g;
	public Adventurer( Graph base ) {
		g = base;
		make_all_basic_matrices();
	}

	public void make_all_basic_matrices() {
		int max = g.numbering.size();
		g.matrix = new Matrix[Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT];
		float[] threeArgs = new float[] { (float)Math.random()*max, (float)Math.random()*max, (float)Math.random()*max };
		for(int i = 0; i < Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT;i++) {
			Algebraicist.form(g, i);
		}
	}

	public Matrix make_random_filter(int rand_base_matrix1, int rand_base_matrix2, int rand_power, int rand_binary ) {
		//System.out.println("RF: " + rand_base_matrix1 + "," + rand_base_matrix2 + "," + rand_power + "," + rand_binary);
		Matrix f = new Matrix(g.matrix[rand_base_matrix1]);
		int max = g.numbering.size();
		float[] threeArgs = new float[] { (float)Math.random()*max, (float)Math.random()*max, (float)Math.random()*max };
		f = Arithmetician.fastMultiproduct(f, rand_power);
		f = Arithmetician.ApplyBinary(rand_binary, f, g.matrix[rand_base_matrix2], threeArgs);
		return f;
	}

	public Matrix make_random_Matrix_filter(int mod_or_not, int[] param_Base_Matrix, int[] param_Unary_Ops, int[] param_Binary_Combinators, int[] zeroOne, float[] threeArgs) {
		// We choose a filter length of max 5 matrices
		// We choose five matrices, and five unary filters to apply to each matrix
		// We choose 4 binary combinator ops to combine these 5 matrices
		// We collapse these ops right to left to construct a random filter
		// But it is only random in the sense that we choose our exploration of the filterspace randomly (or using a basic GA)
		// It is deterministic on the parameters passed in here
		// Which in GA are 'genetic code' / feature vector
		Matrix[] theFive = new Matrix[5];
		int i = 5;
		while(--i >= 0)
			theFive[i] = g.matrix[param_Base_Matrix[i]];
		if(mod_or_not == 0) {
			i = 5;
			while(--i >= 0)
				theFive[i] = Arithmetician.ApplyUnary(param_Unary_Ops[i], theFive[i], threeArgs, zeroOne );
			i = 4;
			Matrix result = theFive[4];
			while(--i >= 0)
				result = Arithmetician.ApplyBinary(param_Binary_Combinators[i], theFive[i], result, threeArgs);
			return result;
		}
		else {
			i = 5;
			while(--i >= 0)
				theFive[i] = Arithmetician.ApplyModUnary(param_Unary_Ops[i], theFive[i], threeArgs, zeroOne );
			i = 4;
			Matrix result = theFive[4];
			while(--i >= 0)
				result = Arithmetician.ApplyModBinary(param_Binary_Combinators[i], theFive[i], result, threeArgs);
			return result;
		}
	}
}

import java.util.HashMap;
import java.util.HashSet;
import java.util.Iterator;
import java.util.Set;

// He makes stuff
public class Algebraicist {
	// 5 primitive matrices

	public static Matrix form( Graph g, int matrixType ) {
		switch(matrixType) {
		case(Graph.MatrixTypeIndex.ADJACENCY):
			return g.matrix[Graph.MatrixTypeIndex.ADJACENCY] = adjacency(g);
		case(Graph.MatrixTypeIndex.DEGREE):
			return g.matrix[Graph.MatrixTypeIndex.DEGREE] = degree(g);
		case(Graph.MatrixTypeIndex.LAPLACIAN):
			return g.matrix[Graph.MatrixTypeIndex.LAPLACIAN] = laplacian(g.matrix[Graph.MatrixTypeIndex.DEGREE], g.matrix[Graph.MatrixTypeIndex.ADJACENCY], 0.5f);
		case(Graph.MatrixTypeIndex.COLLOCATION):
			return g.matrix[Graph.MatrixTypeIndex.COLLOCATION] = collocation(g.matrix[Graph.MatrixTypeIndex.ADJACENCY]);
		default:
			return new Matrix(Matrix.IDENTITY,g.numbering.size(),g.numbering.size());
		}
	}

	public static Matrix adjacency( Graph g ) {
		int vertices = g.edgesByVertexA.size();
		Set<Graph.Vertex> verticeSet = g.edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		Matrix am = new Matrix(Matrix.NONE, vertices, vertices);
		HashMap<Graph.Vertex,Integer> numbering = g.numbering();
		while(vertexIterator.hasNext()) {
			Graph.Vertex v = vertexIterator.next();
			int vnum = numbering.get(v);
			Set<Graph.Edge> vOutEdges = g.edgesByVertexA.get(v);
			Set<Graph.Edge> vInEdges = g.edgesByVertexB.get(v);
			Iterator<Graph.Edge> vOutIterator = vOutEdges.iterator();
			while(vOutIterator.hasNext()) {
				Graph.Edge e = vOutIterator.next();
				Graph.Vertex b = e.b;
				int bnum = numbering.get(b);
				am.entry[vnum][bnum] = 1.0f;
			}
			Iterator<Graph.Edge> vInIterator = vInEdges.iterator();
			while(vInIterator.hasNext()) {
				Graph.Edge e = vInIterator.next();
				Graph.Vertex a = e.a;
				int anum = numbering.get(a);
				am.entry[anum][vnum] = 1.0f;
			}
		}
		am.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return am;
	}

	public static Matrix collocation( Matrix a ) {
		int i = a.rows;
		int j = a.columns;
		Matrix h = new Matrix(Matrix.COLLOCATION, i, j);
		while(--i >= 0) {
			h.hentry[i][i] = neighbourintersection(i,i,a);
			while(--j > i) {
				h.hentry[i][j] = neighbourintersection(i,j,a);
				h.hentry[j ][i] = h.hentry[i][j];
			}  
			j = a.columns;                
		}
		return h;
	}

	public static HashSet neighbourintersection(int row1, int row2, Matrix a) {
		int column = a.columns;
		HashSet<Integer> ni = new HashSet<Integer>();
		while(--column >= 0) {
			float a1c = a.entry[row1][column];
			float a2c = a.entry[row2][column];
			if(a1c == 1.0f && a2c == 1.0f)
				ni.add(new Integer(column));  
		}
		return ni;
	}

	public static Matrix incidence( Graph g ) { // VxE
		int vertices = g.numbering.size();
	int edges = g.edgeNumbering.size();
	Matrix im = new Matrix(Matrix.NONE, vertices, edges );
	Set<Graph.Edge> edgeSet = g.edgeNumbering.keySet();
	Iterator<Graph.Edge> edgeIterator = edgeSet.iterator();
	while(edgeIterator.hasNext()) {
		Graph.Edge e = edgeIterator.next();
		int edgnum = g.edgeNumbering.get(e);
		Graph.Vertex a = e.a;
		Graph.Vertex b = e.b;
		int anum = g.numbering.get(a);
		int bnum = g.numbering.get(b);
		im.entry[anum][edgnum] = 1.0f;
		im.entry[bnum][edgnum] = 1.0f;
	}
	return im;
	}

	public static Matrix SignedIncidence( Graph g ) { // VxE
		int vertices = g.numbering.size();
		int edges = g.edgeNumbering.size();
		Matrix im = new Matrix(Matrix.NONE, vertices, edges );
		Set<Graph.Edge> edgeSet = g.edgeNumbering.keySet();
		Iterator<Graph.Edge> edgeIterator = edgeSet.iterator();
		while(edgeIterator.hasNext()) {
			Graph.Edge e = edgeIterator.next();
			int edgnum = g.edgeNumbering.get(e);
			Graph.Vertex a = e.a;
			Graph.Vertex b = e.b;
			int anum = g.numbering.get(a);
			int bnum = g.numbering.get(b);
			im.entry[anum][edgnum] = 1.0f;
			im.entry[bnum][edgnum] = -1.0f;
		}
		return im;
	}

	public static Matrix degree( Graph g ) {
		int vertices = g.edgesByVertexA.size();
		Matrix dg = new Matrix(Matrix.NONE, vertices, vertices);
		HashMap<Graph.Vertex,Integer> numbering = g.numbering();
		Set<Graph.Vertex> verticeSet = g.edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		g.volume = 0.0f;
		while(vertexIterator.hasNext()) {
			Graph.Vertex v = vertexIterator.next();
			int vnum = numbering.get(v);
			Set<Graph.Edge> vOutEdges = g.edgesByVertexA.get(v);
			Set<Graph.Edge> vInEdges = g.edgesByVertexB.get(v);
			g.volume += dg.entry[vnum][vnum] = vOutEdges.size() + vInEdges.size();
		}
		dg.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return dg;
	}

	public static Matrix edgeAdjacency( Graph g ) { // ExE
		Set<Graph.Vertex> verticeSet = g.edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		int edges = g.edgeNumbering.size();
		Matrix eda = new Matrix(Matrix.NONE, edges, edges);
		while(vertexIterator.hasNext()) {
			Graph.Vertex stupidFuckFace = vertexIterator.next();
			Set<Graph.Edge> fuckingMoron = g.edgesByVertexA.get(stupidFuckFace);
			Iterator<Graph.Edge> edgeIterator = fuckingMoron.iterator();
			while(edgeIterator.hasNext()) {
				Graph.Edge e = edgeIterator.next();
				int edgnum = g.edgeNumbering.get(e);
				Graph.Vertex a = e.a;
				Set<Graph.Edge> heIsStillAFuckFace = g.edgesByVertexA.get(a);
				Iterator<Graph.Edge> fuckIterator = heIsStillAFuckFace.iterator();
				while(fuckIterator.hasNext()) {
					Graph.Edge eneighbour = fuckIterator.next();
					int eneighbournum = g.edgeNumbering.get(eneighbour);
					if(edgnum != eneighbournum)
						eda.entry[edgnum][eneighbournum] = 1.0f;
				}
				if(a!=e.b) {
					Graph.Vertex b = e.b;
					heIsStillAFuckFace = g.edgesByVertexA.get(b);
					fuckIterator = heIsStillAFuckFace.iterator();
					while(fuckIterator.hasNext()) {
						Graph.Edge eneighbour = fuckIterator.next();
						int eneighbournum = g.edgeNumbering.get(eneighbour);
						if(edgnum != eneighbournum)
							eda.entry[edgnum][eneighbournum] = 1.0f;
					}	
				}
			}
		}
		eda.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return eda;
	}

	public static Matrix edgeDegree( Graph g ) { // ExE
		Set<Graph.Vertex> verticeSet = g.edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		int edges = g.edgeNumbering.size();
		Matrix edd = new Matrix(Matrix.NONE, edges, edges);
		while(vertexIterator.hasNext()) {
			Graph.Vertex stupidFuckFace = vertexIterator.next();
			Set<Graph.Edge> fuckingMoron = g.edgesByVertexA.get(stupidFuckFace);
			Iterator<Graph.Edge> edgeIterator = fuckingMoron.iterator();
			while(edgeIterator.hasNext()) {
				Graph.Edge e = edgeIterator.next();
				int edgnum = g.edgeNumbering.get(e);
				Graph.Vertex a = e.a;
				Set<Graph.Edge> heIsStillAFuckFace = g.edgesByVertexA.get(a);
				Iterator<Graph.Edge> fuckIterator = heIsStillAFuckFace.iterator();
				while(fuckIterator.hasNext()) {
					Graph.Edge eneighbour = fuckIterator.next();
					int eneighbournum = g.edgeNumbering.get(eneighbour);
					if(edgnum != eneighbournum)
						edd.entry[edgnum][edgnum] += 1.0f;
				}
				if(a!=e.b) {
					Graph.Vertex b = e.b;
					heIsStillAFuckFace = g.edgesByVertexA.get(b);
					fuckIterator = heIsStillAFuckFace.iterator();
					while(fuckIterator.hasNext()) {
						Graph.Edge eneighbour = fuckIterator.next();
						int eneighbournum = g.edgeNumbering.get(eneighbour);
						if(edgnum != eneighbournum)
							edd.entry[edgnum][edgnum] += 1.0f;
					}	
				}
			}
		}
		edd.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return edd;
	}


	// 5 secondary matrix operators

	public static Matrix complement( Matrix a ) {
		int i = a.rows;
		int j = a.columns;
		Matrix c = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if( a.entry[i][j] != 0.0f) 
					c.entry[i][j] = 0.0f;
				else
					c.entry[i][j] = 1.0f;
			}
			j = a.columns;
		}
		return c;
	}

	public static Matrix exactlyTwoStepAdjacency( Matrix a ) {
		int i = a.rows;
		int j = a.columns;
		int shortest = a.shortest;
		Matrix a2 = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				int k = shortest;
				while ( --k >= 0 ) {
					if( a.entry[i][k] != 0.0f && a.entry[k][j] != 0.0f )
						a2.entry[i][j] = a.entry[i][k] + a.entry[k][j];
				}
			}
			j = a.columns;
		}
		return a2;
	}

	public static Matrix exactlyOneStepDegree( Matrix d, Matrix a ) {
		int i = a.rows;
		int j = a.columns;
		Matrix d1 = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if( a.entry[i][j] != 0.0f) 
					d1.entry[i][i] += d.entry[j][j];
			}
			j = a.columns;
		}
		return d1;
	}

	public static Matrix exactlyTwoStepDegree( Matrix d, Matrix a2 ) {
		int i = a2.rows;
		int j = a2.columns;
		Matrix d2 = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if( a2.entry[i][j] != 0.0f) 
					d2.entry[i][i] += d.entry[j][j];
			}
			j = a2.columns;
		}
		return d2;
	}

	public static Matrix seidelAdjacency( Matrix a ) {
		int i = a.rows;
		int j = a.columns;
		Matrix sa = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if(i==j)
					sa.entry[i][i] = 0.0f;
				else if(a.entry[i][j] != 0.0f)
					sa.entry[i][j] = -1.0f;
				else
					sa.entry[i][j] = 1.0f;
			}
			j = a.columns;
		}
		return sa;
	}

	// 22 Matrix combinators 

	public static Matrix augmentedAdjacencyMatrix( Matrix incidence ) {
		return Arithmetician.productResult(incidence, incidence.transpose());
	}

	public static Matrix FanChungBoundaryOperator( Matrix SignedIncidence, Matrix degree ) {
		int i = SignedIncidence.rows;
		int j = SignedIncidence.columns;
		Matrix sd = Arithmetician.entrywiseSquareRootResult(degree);
		Matrix b = new Matrix(Matrix.NONE, j, i);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				b.entry[j][i] = SignedIncidence.entry[i][j]/sd.entry[i][i];
			}
			j = SignedIncidence.columns;
		}
		return b;
	}

	public static Matrix shareOfDegree( Graph g, Matrix degree ) {
		int i = degree.rows;
		int j = degree.columns;
		Matrix sd = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				sd.entry[i][j] = degree.entry[i][j]/g.volume;
			}
			j = degree.columns;
		}
		return sd;
	}

	public static Matrix minimumDistanceMatrix(Matrix connectivityMetric) {
		return Arithmetician.FloydWarshall(connectivityMetric, 1.0f, 50.0f, connectivityMetric.shortest*4.0f*50.0f, 0, Arithmetician.MIN );
	}


	public static Matrix minimumDistanceMatrixParameterized(Matrix connectivityMetric, float conMetParam) {
		// connectivityMetric.show();       
		return Arithmetician.FloydWarshall(connectivityMetric, conMetParam, 50.0f, connectivityMetric.shortest*4.0f*50.0f, 0, Arithmetician.MIN );
	}

	public static Matrix generalizedAdjacency( float x, float y, float z, Matrix a ) {
		//(xI +yA + z(J-I-A))
		return Arithmetician.pairwiseSumResult(
				Arithmetician.scalarMultiplyResult( 
						new Matrix(Matrix.IDENTITY, a.rows, a.columns), x), 
						Arithmetician.pairwiseSumResult(
								Arithmetician.scalarMultiplyResult(a, y),
								Arithmetician.scalarMultiplyResult(
										Arithmetician.pairwiseSubtractResult(
												new Matrix(Matrix.COUNTERIDENTITY, a.rows, a.columns),
												Arithmetician.pairwiseSumResult(new Matrix(Matrix.IDENTITY, a.rows, a.columns), a)),
												z)
								));

	}

	public static Matrix weightedDegree( Matrix degree, float weight ) {
		return Arithmetician.scalarMultiplyResult(degree, weight);
	}

	public static Matrix weightedAdjacency( Matrix adjacency, float weight ) {
		return Arithmetician.scalarMultiplyResult(adjacency, weight);
	}

	public static Matrix SpielmanLazywalkMatrix(Matrix degree, Matrix adjacency) {
		return Arithmetician.scalarMultiplyResult(
				Arithmetician.pairwiseSumResult(
						new Matrix(Matrix.IDENTITY, adjacency.rows, adjacency.columns),
						Arithmetician.productResult(adjacency, 
								Arithmetician.Inverse(degree))), 
								0.5f);
	}

	public static Matrix lazywalkMatrixSoftInverse(Matrix degree, Matrix adjacency) {
		return Arithmetician.scalarMultiplyResult(
				Arithmetician.pairwiseSumResult(
						new Matrix(Matrix.IDENTITY, adjacency.rows, adjacency.columns),
						Arithmetician.productResult(degree, 
								Arithmetician.Inverse(adjacency))), 
								2.0f);
	}

	public static Matrix nStepWalks( Matrix adjacency, int nsteps ) {
		return Arithmetician.fastMultiproduct(adjacency, nsteps);
	}

	public static Matrix transitionMatrixRandomWalk( Matrix a, Matrix degree ) {
		int i = a.rows;
		int j = a.columns;
		Matrix t = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if( a.entry[i][j] != 0.0f )
					t.entry[i][j] = 1.0f/degree.entry[i][i];
			}
			j = a.columns;
		}
		return t;
	}

	public static Matrix laplacian(Matrix degree, Matrix adjacency, float weight) {
		if(weight != 1.0f)
			//return weightByDegree(degree, Arithmetician.pairwiseSubtractResult(degree,Arithmetician.scalarMultiplyResult(
			//		 adjacency,
			//		weight)));
			return Arithmetician.pairwiseSubtractResult(degree,Arithmetician.scalarMultiplyResult(
					adjacency,
					weight));
		else
			return 	Arithmetician.pairwiseSubtractResult(degree, adjacency);
		//return 	weightByDegree(degree, Arithmetician.pairwiseSubtractResult(degree, adjacency));
	}

	public static Matrix laplaceOperator( Matrix d ) {
		int i = d.rows;
		int j = d.columns;
		Matrix lo = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if(i==j)
					lo.entry[i][i] = 1.0f;
				else
					lo.entry[i][j] = -1.0f/d.entry[i][j]; 
			}
			j = d.columns;
		}
		return lo;
	}

	public static Matrix weightByDegree(Matrix degree, Matrix base) {
		int i = base.rows;
		int j = base.columns;
		Matrix pm = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			float degi = degree.entry[i][i];
			while( --j >= 0 ) {
				if(i==j)
					pm.entry[i][i] = base.entry[i][i];
				else
					pm.entry[i][j] = base.entry[i][j]*degi*degree.entry[j][j]; 
			}
			j = base.columns;
		}  
		return pm;
	}

	public static Matrix SoftNormalizedLaplacian(Matrix degree, Matrix laplacian, float weight ) {
		Matrix entryWiseRoot = Arithmetician.entrywiseSquareRootResult(degree);
		Matrix entryWiseRootInverse = Arithmetician.entrywiseInversionResult(entryWiseRoot);
		if(weight != 1.0f)
			return Arithmetician.productResult(
					Arithmetician.productResult(entryWiseRootInverse, Arithmetician.scalarMultiplyResult(laplacian, weight)),
					entryWiseRoot);
		else
			return Arithmetician.productResult(
					Arithmetician.productResult(entryWiseRootInverse, laplacian),
					entryWiseRoot);
	}



	public static Matrix HardNormalizedLaplacian(Matrix degree, Matrix laplacian, float weight ) {
		Matrix[] rootAndInverse = Arithmetician.squareRootInverseSquareRootResult(degree);
		if(weight != 1.0f)
			return Arithmetician.productResult(
					Arithmetician.productResult(rootAndInverse[1], Arithmetician.scalarMultiplyResult(laplacian, weight)),
					rootAndInverse[0]);

		else
			return Arithmetician.productResult(
					Arithmetician.productResult(rootAndInverse[1], laplacian),
					rootAndInverse[0]);
	}

	public static Matrix eigenvectors( Matrix m ) { 
		Arithmetician.PowerIterationDeflate(m, m.shortest, 0.000000001, 0 );
		Matrix e = new Matrix(m.eigenvector);
		System.arraycopy(m.eigenvalue, 0, e.eigenvalue, 0, m.eigenvalue.length);
		return e;
	}

	public static void show_v (float[] v) {
		int i = v.length;
		while(--i >= 0)
			System.out.print(v[i] + " ");
		System.out.println();
	}


	public static Matrix commuteTime( Graph g ) {  
		Matrix nl = g.matrix[Graph.MatrixTypeIndex.LAPLACIAN];
		Matrix d = Arithmetician.entrywiseSquareRootResult(g.matrix[Graph.MatrixTypeIndex.DEGREE]);
		Matrix eigns = eigenvectors(nl);
		Matrix ct = new Matrix(Matrix.NONE, d.rows, d.columns);
		float vol = g.volume;
		int k = nl.shortest;
		while(--k >= 0 ) {
			int i = d.rows;
			while(--i >= 0) {
				int j = d.columns;
				while(--j >= 0) {
					float eki = eigns.entry[k][i]/d.entry[i][i];
					float ekj = eigns.entry[k][j]/d.entry[j][j];
					ct.entry[i][j] += vol*(eki-ekj)*(eki-ekj)/eigns.eigenvalue[k];
				}
			}
		}
		return ct;
	}

	public static Matrix deformedLaplacian(Matrix d, Matrix a, float deformation ) {
		//I - sA + s2(D - I)
		return Arithmetician.pairwiseSumResult(
				Arithmetician.pairwiseSubtractResult(
						new Matrix(Matrix.IDENTITY, a.rows, a.columns ),
						Arithmetician.scalarMultiplyResult(a, deformation)),
						Arithmetician.scalarMultiplyResult(
								Arithmetician.pairwiseSubtractResult(d, new Matrix(Matrix.IDENTITY, d.rows, d.columns)),
								deformation*deformation));

	}

	public static Matrix projectionMatrix( int rows, int columns ) {
		float ns = (float)Math.sqrt(rows*columns);
		float diag = (ns-1.0f)/ns;
		float nondiag = -1.0f/ns;
		int i = rows;
		int j = columns;
		Matrix pm = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if(i==j)
					pm.entry[i][i] = diag;
				else
					pm.entry[i][j] = nondiag; 
			}
			j = columns;
		}
		return pm;
	}

	public static Matrix projectedMatrix(Matrix embeddingMetric, float projectionCoefficient ) { // 0.5f
		// p = 0.5PEP
		return Arithmetician.productResult(
				Arithmetician.scalarMultiplyResult(projectionMatrix(embeddingMetric.rows, embeddingMetric.columns), projectionCoefficient),
				Arithmetician.productResult( embeddingMetric, projectionMatrix(embeddingMetric.rows, embeddingMetric.columns)));
	}

	public static Matrix RPaugmentedDegreeMatrix(Matrix d, Matrix mindistance) {
		Matrix mindistanceRoot = Arithmetician.entrywiseSquareRootResult(mindistance);
		int i = d.rows;
		int j = d.columns;
		Matrix augd = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if(i==j)
					augd.entry[i][i] = d.entry[i][i];
				else if ((1 << (int)(mindistanceRoot.entry[i][j])) != 0.0f )
					augd.entry[i][j] = d.entry[j][j] / (1 << (int)(mindistanceRoot.entry[i][j])) ;
			}
			j = d.columns;
		}
		return augd;
	}

	public static Matrix maximumDistanceMatrix(Matrix connectivityMetric) {
		return Arithmetician.FloydWarshall(connectivityMetric, 1.0f, 50.0f, connectivityMetric.shortest*4.0f*50.0f, 1, Arithmetician.MAX );
	}












}

import java.util.Arrays;

public final class Arithmetician {
	public static final int MIN = 0;
	public static final int MAX = 1;
	public int parametrization = (2*3*10*18*9*32);

	public class MatrixOperation {
		public static final int UNARY = 0;
		public static final int SCALAR = 1;
		public static final int BINARY = 2;
		public static final int MOD = 3;

		public class MatrixBinary {
			public static final int PAIRWISESUM = 0;
			public static final int PAIRWISESUBTRACT = 1;
			public static final int HADARMARDPRODUCT = 2;
			public static final int PRODUCTPROPAGATED = 3;
			public static final int PRODUCT = 4;
			public static final int HADAMARDDIVISION = 5;
			public static final int AND = 6;
			public static final int OR = 7;
			public static final int XOR = 8;
			public static final int KOR = 9;
			public static final int FROBENIUSINNERPRODUCT = 10;
			public static final int BINARY_FIELD_COUNT = 11;

		}
		public class MatrixUnary {
			public static final int ENTRYWISEINVERSION = 1;
			public static final int EXPONENTIALTRACE = 2;
			public static final int SQUAREROOTINVERSESQUAREROOT = 3;
			public static final int SOFTSQUAREROOT = 4;
			public static final int ENTRYWISENATURALLOG = 5;
			public static final int RANDOMFILL = 6;
			public static final int L0NORMALIZE = 7;
			public static final int L1NORMALIZE = 8;
			public static final int L2NORMALIZE = 9;
			public static final int LINFNORMALIZE = 10;
			public static final int FLOYDWARSHALL = 11;
			public static final int POWERITERATIONDEFLATE = 12;
			public static final int LUDOLITTLEDECOMPOSITION = 13;
			public static final int DETERMINANT = 14;
			public static final int INVERSE = 15;
			public static final int NOT = 16;
			public static final int COMPLEMENT = 17;
			public static final int UNARY_FIELD_COUNT = 18;
		}
		public class MatrixScalar {
			public static final int SCALARADD = 0;
			public static final int SCALARSUBTRACT = 1;
			public static final int SCALARMULTIPLY = 2;
			public static final int SCALARDIVIDE = 3;
			public static final int HADAMARDPOWER = 4;		
			public static final int FASTPROPAGATEDMULTIPRODUCT = 5;
			public static final int FASTMULTIPRODUCT = 6;
			public static final int ENTRYWISEMAPTORANGE = 7;
			public static final int THRESHOLD = 8;

		}

		public class Mod {
			public static final int UNARY = 0;
			public static final int SCALAR = 1;
			public static final int BINARY = 2;
			public class MatrixBinary {
				public static final int PAIRWISESUM = 0;
				public static final int PAIRWISESUBTRACT = 1;
				public static final int HADARMARDPRODUCT = 2;
				public static final int PRODUCTPROPAGATED = 3;
				public static final int PRODUCT = 4;
				public static final int HADAMARDDIVISION = 5;
				public static final int AND = 6;
				public static final int OR = 7;
				public static final int XOR = 8;
				public static final int KOR = 9;

			}
			public class MatrixUnary {
				public static final int FROBENIUSINNERPRODUCT = 0;
				public static final int ENTRYWISEINVERSION = 1;
				public static final int EXPONENTIALTRACE = 2;
				public static final int SQUAREROOTINVERSESQUAREROOT = 3;
				public static final int SOFTSQUAREROOT = 4;
				public static final int ENTRYWISENATURALLOG = 5;
				public static final int RANDOMFILL = 6;
				public static final int L0NORMALIZE = 7;
				public static final int L1NORMALIZE = 8;
				public static final int L2NORMALIZE = 9;
				public static final int LINFNORMALIZE = 10;
				public static final int FLOYDWARSHALL = 11;
				public static final int POWERITERATIONDEFLATE = 12;
				public static final int LUDOLITTLEDECOMPOSITION = 13;
				public static final int DETERMINANT = 14;
				public static final int INVERSE = 15;
				public static final int NOT = 16;
				public static final int COMPLEMENT = 17;
			}
			public class MatrixScalar {
				public static final int SCALARADD = 0;
				public static final int SCALARSUBTRACT = 1;
				public static final int SCALARMULTIPLY = 2;
				public static final int SCALARDIVIDE = 3;
				public static final int HADAMARDPOWER = 4;		
				public static final int FASTPROPAGATEDMULTIPRODUCT = 5;
				public static final int FASTMULTIPRODUCT = 6;
				public static final int ENTRYWISEMAPTORANGE = 7;
				public static final int THRESHOLD = 8;

			}
		}

	}

	public static void vectorScalarAdd( float[] v, float add ) {
		int i = v.length;
		while(--i >= 0)
			v[i]+=add;
	}

	public static Matrix ApplyScalar( int opNum, Matrix m, float[] threeArgs) {
		switch(opNum) {
		case(MatrixOperation.MatrixScalar.SCALARADD):
			return scalarAddResult(m, threeArgs[0]);
		case(MatrixOperation.MatrixScalar.SCALARSUBTRACT):
			return scalarSubtractResult(m, threeArgs[0]);
		case(MatrixOperation.MatrixScalar.SCALARMULTIPLY):
			return scalarMultiplyResult(m, threeArgs[0]);
		case(MatrixOperation.MatrixScalar.SCALARDIVIDE):
			return scalarDivideResult(m, threeArgs[0]);
		case(MatrixOperation.MatrixScalar.HADAMARDPOWER):		
			return hadamardPowerResult(m, (int)threeArgs[0]);
		case(MatrixOperation.MatrixScalar.FASTPROPAGATEDMULTIPRODUCT):
			fastPropagatedMultiproduct(m, (int)threeArgs[0]);
		return m;
		case(MatrixOperation.MatrixScalar.FASTMULTIPRODUCT):
			return fastMultiproduct(m, (int)threeArgs[0]);
		case(MatrixOperation.MatrixScalar.ENTRYWISEMAPTORANGE):
			return scaleToRangeResult(m, threeArgs[0], threeArgs[1]);
		case(MatrixOperation.MatrixScalar.THRESHOLD):
			return thresholdResult(m, threeArgs[0], (int)threeArgs[1]);
		default:
			return new Matrix(m);
		}
	}

	public static Matrix ApplyModScalar( int opNum, Matrix m, float[] threeArgs) {
		switch(opNum) {
		case(MatrixOperation.MatrixScalar.SCALARADD):
			return ModscalarAddResult(m, threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.SCALARSUBTRACT):
			return ModscalarSubtractResult(m, threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.SCALARMULTIPLY):
			return ModscalarMultiplyResult(m, threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.SCALARDIVIDE):
			return ModscalarDivideResult(m, threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.HADAMARDPOWER):		
			return ModhadamardPowerResult(m, (int)threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.FASTPROPAGATEDMULTIPRODUCT):
			ModfastPropagatedMultiproduct(m, (int)threeArgs[0], threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixScalar.FASTMULTIPRODUCT):
			return ModfastMultiproduct(m, (int)threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.ENTRYWISEMAPTORANGE):
			return ModscaleToRangeResult(m, threeArgs[0], threeArgs[1], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.THRESHOLD):
			return ModthresholdResult(m, threeArgs[0], (int)threeArgs[1], threeArgs[2]);
		default:
			return new Matrix(m);
		}
	}


	public  static Matrix ApplyUnary( int opNum, Matrix m, float[] threeArgs, int[] zeroOne) {
		switch(opNum) {
		case(MatrixOperation.MatrixUnary.ENTRYWISEINVERSION):
			return entrywiseInversionResult(m);
		case(MatrixOperation.MatrixUnary.EXPONENTIALTRACE):
			float result = exponentialTrace(m);
		return scalarMultiplyResult(m,result);
		case(MatrixOperation.MatrixUnary.SQUAREROOTINVERSESQUAREROOT):
			return squareRootInverseSquareRootResult(m)[0];
		case(MatrixOperation.MatrixUnary.SOFTSQUAREROOT):
			return entrywiseSquareRootResult(m);
		case(MatrixOperation.MatrixUnary.ENTRYWISENATURALLOG):
			return entrywiseLogarithmResult(m);
		case(MatrixOperation.MatrixUnary.RANDOMFILL):
			fill_random(m);
		return m;
		case(MatrixOperation.MatrixUnary.L0NORMALIZE):
			L0normalize(m);
		return m;
		case(MatrixOperation.MatrixUnary.L1NORMALIZE):
			L1normalize(m);
		return m;
		case(MatrixOperation.MatrixUnary.L2NORMALIZE):
			L2normalize(m);
		return m;
		case(MatrixOperation.MatrixUnary.LINFNORMALIZE):
			Linfnormalize(m);
		return m;
		case(MatrixOperation.MatrixUnary.FLOYDWARSHALL):
			return FloydWarshall(m, 0.00000001f, threeArgs[0], threeArgs[1], zeroOne[0], zeroOne[1]);
		case(MatrixOperation.MatrixUnary.POWERITERATIONDEFLATE):
			m.eigenvalue = new float[m.rows];
		m.eigenvector = new float[m.rows][m.columns];
		PowerIterationDeflate(m, m.shortest, (double)0.00000001f, 1); // 1 populates krylov subspace
		return new Matrix(m.eigenvector);
		case(MatrixOperation.MatrixUnary.LUDOLITTLEDECOMPOSITION):
			return LU_Dolittle_DecompositionResult(m);
		case(MatrixOperation.MatrixUnary.DETERMINANT):
			result = Determinant(m);
		return scalarMultiplyResult(m,result);
		case(MatrixOperation.MatrixUnary.INVERSE):
			return Inverse(m);
		case(MatrixOperation.MatrixUnary.NOT):
			return NOTResult(m);
		case(MatrixOperation.MatrixUnary.COMPLEMENT):
			return complement(m);
		default:
			return new Matrix(m);
		}
	}

	public static Matrix ApplyModUnary( int opNum, Matrix m, float[] threeArgs, int[] zeroOne) {
		switch(opNum) {
		case(MatrixOperation.MatrixUnary.ENTRYWISEINVERSION):
			return ModentrywiseInversionResult(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.EXPONENTIALTRACE):
			float result = ModexponentialTrace(m, threeArgs[2]);
		return ModscalarMultiplyResult(m, result, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.SQUAREROOTINVERSESQUAREROOT):
			return ModsquareRootInverseSquareRootResult(m, threeArgs[2])[0];
		case(MatrixOperation.MatrixUnary.SOFTSQUAREROOT):
			return ModentrywiseSquareRootResult(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.ENTRYWISENATURALLOG):
			return ModentrywiseLogarithmResult(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.RANDOMFILL):
			Modfill_random(m, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixUnary.L0NORMALIZE):
			ModL0normalize(m, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixUnary.L1NORMALIZE):
			ModL1normalize(m, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixUnary.L2NORMALIZE):
			ModL2normalize(m, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixUnary.LINFNORMALIZE):
			ModLinfnormalize(m, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixUnary.FLOYDWARSHALL):
			return ModFloydWarshall(m, 0.00000001f, threeArgs[0], threeArgs[1], zeroOne[0], zeroOne[1], threeArgs[2]);
		case(MatrixOperation.MatrixUnary.POWERITERATIONDEFLATE):
			ModPowerIterationDeflate(m, m.shortest, (double)0.00000001f, 1, threeArgs[2]); // 1 populates krylov subspace
		return new Matrix(m.eigenvector);
		case(MatrixOperation.MatrixUnary.LUDOLITTLEDECOMPOSITION):
			return ModLU_Dolittle_DecompositionResult(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.DETERMINANT):
			result = ModDeterminant(m, threeArgs[2]);
		return ModscalarMultiplyResult(m,result, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.INVERSE):
			return ModInverse(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.NOT):
			return ModNOTResult(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.COMPLEMENT):
			return complement(m);
		default:
			return new Matrix(m);
		}
	}


	public  static Matrix ApplyBinary( int opNum, Matrix m, Matrix m2, float[] threeArgs) {
		switch(opNum) {
		case(MatrixOperation.MatrixBinary.FROBENIUSINNERPRODUCT):
			float result = FrobeniusInnerProduct(m, m2);
		return scalarMultiplyResult(m,result);
		case(MatrixOperation.MatrixBinary.PAIRWISESUM):
			return pairwiseSumResult(m, m2);
		case(MatrixOperation.MatrixBinary.PAIRWISESUBTRACT):
			return pairwiseSubtractResult(m, m2);
		case(MatrixOperation.MatrixBinary.HADARMARDPRODUCT):
			return hadamardProductResult(m, m2);
		case(MatrixOperation.MatrixBinary.PRODUCTPROPAGATED):
			productPropagated(m, m2);
		return m;
		case(MatrixOperation.MatrixBinary.PRODUCT):
			return productResult(m, m2);
		case(MatrixOperation.MatrixBinary.HADAMARDDIVISION):
			return hadamardDivisionResult(m, m2);
		case(MatrixOperation.MatrixBinary.AND):
			return ANDResult(m, m2);
		case(MatrixOperation.MatrixBinary.OR):
			return ORResult(m, m2);
		case(MatrixOperation.MatrixBinary.XOR):
			return XORResult(m, m2);
		case(MatrixOperation.MatrixBinary.KOR):
			return KORResult(m, m2);
		default:
			return new Matrix(m);
		}
	}
	//Addition


	public  static Matrix ApplyModBinary( int opNum, Matrix m, Matrix m2, float[] threeArgs) {
		switch(opNum) {
		case(MatrixOperation.MatrixBinary.FROBENIUSINNERPRODUCT):
			float result = ModFrobeniusInnerProduct(m, m2, threeArgs[2]);
		return ModscalarMultiplyResult(m,result, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.PAIRWISESUM):
			return ModpairwiseSumResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.PAIRWISESUBTRACT):
			return ModpairwiseSubtractResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.HADARMARDPRODUCT):
			return ModhadamardProductResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.PRODUCTPROPAGATED):
			ModproductPropagated(m, m2, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixBinary.PRODUCT):
			return ModproductResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.HADAMARDDIVISION):
			return ModhadamardDivisionResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.AND):
			return ModANDResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.OR):
			return ModORResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.XOR):
			return ModXORResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.KOR):
			return ModKORResult(m, m2, threeArgs[2]);
		default:
			return new Matrix(m);
		}
	}

	public static void scalarAdd(Matrix accumulator, float addend) {
		if(addend == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestColumn = accumulator.columns;
		int i = accumulator.rows;
		int j = shortestColumn;
		while(--i >= 0) {
			while(--j >= 0)
				accumulator.entry[i][j] += addend;
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static void scalarAdd(Matrix accumulator, float addend, int row) {
		if(addend == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestColumn = accumulator.columns;
		int i = row;
		int j = shortestColumn;
		while(--j >= 0)
			accumulator.entry[i][j] += addend;
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix scalarAddResult(Matrix adder, float addend) {
		Matrix result;
		if(addend == Matrix.ZEROF) {
			result = new Matrix(adder);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = adder.rows;
		int shortestColumn = adder.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = adder.entry[i][j] + addend;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void pairwiseSum(Matrix accumulator, Matrix summand) {
		int shortestRow = accumulator.rows < summand.rows ? accumulator.rows : summand.rows;
		int shortestColumn = accumulator.columns < summand.columns ? accumulator.columns : summand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] += summand.entry[i][j];
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix pairwiseSumResult(Matrix adder, Matrix summand) {
		int shortestRow = adder.rows < summand.rows ? adder.rows : summand.rows;
		int shortestColumn = adder.columns < summand.columns ? adder.columns : summand.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = adder.entry[i][j] + summand.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Subtraction

	public static void scalarSubtract(Matrix accumulator, float subtractand) {
		if(subtractand == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] -= subtractand;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static void scalarSubtract(Matrix accumulator, float subtractand, int row) {
		if(subtractand == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestColumn = accumulator.columns;
		int i = row;
		int j = shortestColumn;
		while(--j >= 0)
			accumulator.entry[i][j] -= subtractand;
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix scalarSubtractResult(Matrix subtractor, float subtractand) {
		Matrix result;
		if(subtractand == Matrix.ZEROF) {
			result = new Matrix(subtractor);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = subtractor.rows;
		int shortestColumn = subtractor.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = subtractor.entry[i][j] - subtractand;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void pairwiseSubtract(Matrix accumulator, Matrix subtractand) {
		int shortestRow = accumulator.rows < subtractand.rows ? accumulator.rows : subtractand.rows;
		int shortestColumn = accumulator.columns < subtractand.columns ? accumulator.columns : subtractand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] -= subtractand.entry[i][j];
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix pairwiseSubtractResult(Matrix subtractor, Matrix subtractand) {
		int shortestRow = subtractor.rows < subtractand.rows ? subtractor.rows : subtractand.rows;
		int shortestColumn = subtractor.columns < subtractand.columns ? subtractor.columns : subtractand.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = subtractor.entry[i][j] - subtractand.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Multiplication

	public static void scalarMultiply(Matrix accumulator, float multiplier) {
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] *= multiplier;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static void scalarMultiply(Matrix accumulator, float multiplier, int row) {
		int shortestColumn = accumulator.columns;
		int i = row;
		int j = shortestColumn;
		while(--j >= 0)
			accumulator.entry[i][j] *= multiplier;
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix scalarMultiplyResult(Matrix multiplicand, float multiplier) {
		int shortestRow = multiplicand.rows;
		int shortestColumn = multiplicand.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = multiplicand.entry[i][j] * multiplier;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static float FrobeniusInnerProduct( Matrix multiplicand, Matrix multiplier ) {
		int shortestRow = multiplicand.rows < multiplier.rows ? multiplicand.rows : multiplier.rows;
		int shortestColumn = multiplicand.columns < multiplier.columns ? multiplicand.columns : multiplier.columns;
		int i = shortestRow;
		float accumulator = 0.0f;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator += multiplicand.entry[i][j]*multiplier.entry[i][j];
			}  
		}
		return accumulator;
	}

	public static void hadamardProduct(Matrix accumulator, Matrix multiplicand) {
		int shortestRow = accumulator.rows < multiplicand.rows ? accumulator.rows : multiplicand.rows;
		int shortestColumn = accumulator.columns < multiplicand.columns ? accumulator.columns : multiplicand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] *= multiplicand.entry[i][j];
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix hadamardProductResult(Matrix multiplier, Matrix multiplicand) {
		int shortestRow = multiplier.rows < multiplicand.rows ? multiplier.rows : multiplicand.rows;
		int shortestColumn = multiplier.columns < multiplicand.columns ? multiplier.columns : multiplicand.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = multiplier.entry[i][j] * multiplicand.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void productPropagated(Matrix accumulator, Matrix multiplier) {
		int shortestRow = accumulator.rows < multiplier.rows ? accumulator.rows : multiplier.rows;
		int shortestColumn = accumulator.columns < multiplier.columns ? accumulator.columns : multiplier.columns;
		int shortest = shortestRow < shortestColumn ? shortestRow : shortestColumn;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int k = shortest;
				float entrySum = 0.0f;
				while(--k >= 0) {
					entrySum += accumulator.entry[i][k]*multiplier.entry[k][j];
				}
				accumulator.entry[i][j] = entrySum;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix productResult(Matrix multiplicand, Matrix multiplier) { 
		int shortest = multiplicand.columns < multiplier.rows ? multiplicand.columns : multiplier.rows;
		int i = multiplicand.rows;
		int j = multiplier.columns;
		Matrix result = new Matrix(Matrix.NONE, i, j);
		while(--i >= 0) {
			j = multiplier.columns;
			while(--j >= 0) {
				int k = shortest;
				float entrySum = 0.0f;
				while(--k >= 0) {
					entrySum += multiplicand.entry[i][k]*multiplier.entry[k][j];
				}
				result.entry[i][j] = entrySum;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Division

	public static void entrywiseInversion( Matrix divisor ) {
		int shortestRow = divisor.rows;
		int shortestColumn = divisor.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float divisorij = divisor.entry[i][j];
				if(divisorij != Matrix.ZEROF)
					divisor.entry[i][j] = Matrix.IDENTITYF / divisorij;
				else
					divisor.entry[i][j] = Matrix.IDENTITYF;
			}  
		}
	}

	public static Matrix entrywiseInversionResult ( Matrix divisor ) {
		int shortestRow = divisor.rows;
		int shortestColumn = divisor.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float divisorij = divisor.entry[i][j];
				if(divisorij != Matrix.ZEROF)
					result.entry[i][j] = Matrix.IDENTITYF / divisorij;
				else
					result.entry[i][j] = Matrix.IDENTITYF;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void scalarDivide(Matrix accumulator, float divisor) {
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		if(divisor == Matrix.ZEROF) {
			Arrays.fill(accumulator.entry, Matrix.IDENTITYF);
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) 
				accumulator.entry[i][j] /= divisor;
		}
	}

	public static Matrix scalarDivideResult(Matrix dividend, float divisor) {
		int shortestRow = dividend.rows;
		int shortestColumn = dividend.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor != Matrix.ZEROF)
					result.entry[i][j] = dividend.entry[i][j] / divisor;
				else
					result.entry[i][j] = Matrix.IDENTITYF;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void hadamardDivision(Matrix accumulator, Matrix divisor) {
		int shortestRow = accumulator.rows < divisor.rows ? accumulator.rows : divisor.rows;
		int shortestColumn = accumulator.columns < divisor.columns ? accumulator.columns : divisor.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor.entry[i][j] != Matrix.ZEROF)
					accumulator.entry[i][j] /= divisor.entry[i][j];
				else
					accumulator.entry[i][j] = Matrix.ZEROF;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix hadamardDivisionResult(Matrix dividend, Matrix divisor) {
		int shortestRow = dividend.rows < divisor.rows ? dividend.rows : divisor.rows;
		int shortestColumn = dividend.columns < divisor.columns ? dividend.columns : divisor.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor.entry[i][j] != Matrix.ZEROF)
					result.entry[i][j] = dividend.entry[i][j] / divisor.entry[i][j];
				else
					result.entry[i][j] = Matrix.IDENTITYF;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Powers

	public static void hadamardPower(Matrix accumulator, int power) {
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		if(power == 0) {
			accumulator.make_ones();
			return;  
		}
		if(power == 1)
			return;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int p = power;
				while(--p > 0)
					accumulator.entry[i][j] *= accumulator.entry[i][j];
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix hadamardPowerResult(Matrix multiplier, int power) {

		int shortestRow = multiplier.rows;
		int shortestColumn = multiplier.columns;
		if(power == 1)
			return new Matrix(multiplier);
		Matrix result = new Matrix(Matrix.ONES, shortestRow, shortestColumn);
		if(power == 0)
			return result;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int p = power;
				while(--p >= 0)
					result.entry[i][j] *= multiplier.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static int[] radix(int num, int base) {
		int b = 1; int b1 = base;
		int rem = 0; int rem1 = num % b1;
		int bit = 0;
		int[] bits = new int[10];
		while(b <= num) {
			bits[bit] = (rem1 - rem)/b;
			b=b1;b1*=base;
			rem=rem1; rem1 = num % b1;
			bit++;
			if(bit >= bits.length) {
				int[] temp = new int[bit<<1];
				System.arraycopy(temp, 0, bits, 0, bit);
				bits = temp;
			}
		}
		return bits;
	}

	public static float exponentialTrace( Matrix m ) {
		float result = 1.0f;
		float trace = m.trace();
		int itrace = (int)trace;
		float diff = trace - itrace;
		int[] bits = radix(itrace, 2);
		int iteration = 0;
		float lastPow = 2.71828183f;
		int lastPowIndex = 0;
		while(++iteration < bits.length) {
			int bit = bits[iteration];
			if(bit == 1) {
				int squaresRequired = iteration - lastPowIndex;
				while(--squaresRequired >= 0)
					lastPow *= lastPow;
				result *= lastPow;
			}
		}
		float extra = -(float)Math.log(diff);
		if(extra != 0.0f)
			result *= extra;
		return result;
	}

	public static void polysquarePropagated(Matrix base, int times) {
		while(--times >= 0)
			productPropagated(base, base);
	}

	public static void fastPropagatedMultiproduct(Matrix base, int exponent) { 
		int[] operations = radix(exponent, 2);
		int lastProductIndex = 0, jmax = operations.length;
		Matrix partialPowers;
		if(operations[0] == 0)
			partialPowers = new Matrix(Matrix.IDENTITY, base.rows, base.columns);
		else
			partialPowers = new Matrix(base);
		for(int j = 1; j < jmax;j++) {
			int bit = operations[j];
			if(bit == 1) {
				polysquarePropagated(partialPowers, j-lastProductIndex);
				lastProductIndex = j;
				productPropagated(base, partialPowers);
			}
		}
	}

	public static Matrix polysquare(Matrix base, int times) {
		while(--times >= 0)
			base = productResult(base, base);
		return base;
	}

	public static Matrix fastMultiproduct(Matrix base, int exponent) { 
		int[] operations = radix(exponent, 2);
		int lastProductIndex = 0, jmax = operations.length;
		Matrix partialPowers;
		if(operations[0] == 0)
			partialPowers = new Matrix(Matrix.IDENTITY, base.rows, base.columns);
		else
			partialPowers = new Matrix(base);
		for(int j = 1; j < jmax;j++) {
			int bit = operations[j];
			if(bit == 1) {
				polysquare(partialPowers, j-lastProductIndex);
				lastProductIndex = j;
				base = productResult(base, partialPowers);
			}
		}
		return base;
	}

	public static Matrix squareRootReturnInverseSquareRoot(Matrix y) {
		Matrix z = new Matrix(Matrix.IDENTITY, y.rows, y.columns);
		int j = (int)Math.sqrt(y.rows) + 1;
		while(--j >= 0) {
			scalarMultiply(pairwiseSumResult(y,Inverse(z)), 0.5f);
			scalarMultiply(pairwiseSumResult(z,Inverse(y)), 0.5f);
		}
		y.updates &= (Matrix.all_changed^Matrix.entry_changed);
		z.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return z;
	}

	public static Matrix[] squareRootInverseSquareRootResult(Matrix x) {
		Matrix y = new Matrix(x);
		Matrix z = new Matrix(Matrix.IDENTITY, y.rows, y.columns);
		int j = (int)Math.sqrt(y.rows) + 1;
		while(--j >= 0) {
			scalarMultiply(pairwiseSumResult(y,Inverse(z)), 0.5f);
			scalarMultiply(pairwiseSumResult(z,Inverse(y)), 0.5f);
		}
		y.updates &= (Matrix.all_changed^Matrix.entry_changed);
		z.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return new Matrix[] { y, z };
	}

	public static void entrywiseSquareRoot(Matrix x) {
		int shortestRow = x.rows;
		int shortestColumn = x.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float xij = x.entry[i][j];
				float sign;
				if(xij != 0.0f)
					sign = xij/Math.abs(xij);
				else 
					sign = 1.0f;
				x.entry[i][j] = (float)( sign*Math.sqrt(Math.abs(xij)) );
			}  
		}
		x.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix entrywiseSquareRootResult(Matrix x) {
		int shortestRow = x.rows;
		int shortestColumn = x.columns;
		int i = shortestRow;
		Matrix result = new Matrix(Matrix.NONE, x.rows, x.columns);
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float xij = x.entry[i][j];
				float sign;
				if(xij != 0.0f)
					sign = xij/Math.abs(xij);
				else 
					sign = 1.0f;
				result.entry[i][j] = (float)( sign*Math.sqrt(Math.abs(xij)) );
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}	

	// Logarithm

	public static void entrywiseLogarithm( Matrix m ) {
		int shortestRow = m.rows;
		int shortestColumn = m.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float eij = m.entry[i][j];
				float sign;
				if(eij != 0.0f)
					sign = eij/Math.abs(eij);
				else 
					sign = 1.0f;
				m.entry[i][j] = (float)(sign*Math.log(eij));
			}  
		}
		m.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix entrywiseLogarithmResult( Matrix m ) {
		int shortestRow = m.rows;
		int shortestColumn = m.columns;
		int i = shortestRow;
		Matrix result = new Matrix(m);
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float xij = m.entry[i][j];
				float sign;
				if(xij != 0.0f)
					sign = xij/Math.abs(xij);
				else 
					sign = 1.0f;
				result.entry[i][j] = (float)(sign*Math.log(xij));
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Mapping

	public static void scaleToRange(Matrix m, float low, float high) {
		float[] minmax = m.minmax();
		float newrange = high-low;
		float oldrange = minmax[1]-minmax[0];
		float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
		scalarSubtract(m, minmax[0]);
		scalarMultiply(m, range );
		scalarAdd(m, low );
	}

	public static void scaleToRangeRows(Matrix m, float low, float high) {
		int i = m.rows;
		while(--i >= 0) {
			float[] minmax = m.minmaxRow(i);
			float newrange = high-low;
			float oldrange = minmax[1]-minmax[0];
			float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
			scalarSubtract(m, minmax[0], i);
			scalarMultiply(m, range, i );
			scalarAdd(m, low, i );
		}
	}

	public static void scaleToRangeRow(Matrix m, int row, float low, float high) {
		int i = row;
		float[] minmax = m.minmaxRow(i);
		float newrange = high-low;
		float oldrange = minmax[1]-minmax[0];
		float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
		scalarSubtract(m, minmax[0], i);
		scalarMultiply(m, range, i );
		scalarAdd(m, low, i );

	}

	public static Matrix scaleToRangeResult(Matrix m, float low, float high) {
		float[] minmax = m.minmax();
		float newrange = high-low;
		float oldrange = minmax[1]-minmax[0];
		float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
		Matrix result = new Matrix(Matrix.NONE, m.rows, m.columns);
		result = scalarSubtractResult(m, minmax[0]);
		scalarMultiply(result, range );
		scalarAdd(result, low );
		return result;
	}


	// Random

	public static void fill_random(Matrix m) {
		int i = m.rows;
		while(--i >= 0 ) {
			int j = m.columns;
			while(--j >= 0 )
				m.entry[i][j] = (float)Math.random();
		}
	}

	public static void fill_random( float[] v) {
		int j = v.length;
		while(--j >= 0 ) v[j] = (float)Math.random();
	}


	// Norms & normalization vectors and matrices : 0, 1, 2, inf

	public static void Linfnormalize( Matrix m ) {
		int i = m.rows;
		int j = m.columns;
		if( i == 0 || j == 0)
			return;
		float max = Math.abs(m.entry[0][0]);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				float eij = Math.abs(m.entry[i][j]);
				if(eij > max)
					max = eij;
			}
			j = m.columns;
		}
		i = m.rows;
		if(max != 0.0f) {
			while( --i >= 0 ) {
				j = m.columns;
				while( --j >= 0 ) {
					m.entry[i][j] /= max;
				}
			}
		}
	}

	public static void Linfnormalize( float[] v ) {
		int j = v.length;
		if(j == 0)
			return;
		float max = Math.abs(v[0]);
		while(--j > 0) {
			float vj = Math.abs(v[j]);
			if(vj > max)
				max = vj;
		}
		j = v.length;
		if(max == 0.0f)
			return;
		while(--j >= 0) v[j] /= max;
	}

	public static float Linfnorm( Matrix m ) {
		int i = m.rows;
		int j = m.columns;
		if( i == 0 || j == 0)
			return 0.0f;
		float max = Math.abs(m.entry[0][0]);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				float eij = Math.abs(m.entry[i][j]);
				if(eij > max)
					max = eij;
			}
			j = m.columns;
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[3] = max;
		return max;
	}

	public static float Linfnorm( float[] v ) {
		int j = v.length;
		if(j == 0)
			return 0.0f;
		float max = Math.abs(v[0]);
		while(--j > 0) {
			float vj = Math.abs(v[j]);
			if(vj > max)
				max = vj;
		}
		return max;
	}

	public static void L2normalize( Matrix m ) {
		int i = m.rows;
		float sumSquares = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sumSquares += m.entry[i][j]*m.entry[i][j];
			}
		}
		sumSquares = (float)Math.sqrt(sumSquares);
		i = m.rows;
		if(sumSquares == 0.0f)
			return;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sumSquares;
			}
		}
	}

	public static void L2normalize( float [] v ) {
		int j = v.length;
		float sumSquares = 0.0f;
		while(--j >= 0) sumSquares += v[j]*v[j];
		j = v.length;
		sumSquares = (float)Math.sqrt(sumSquares);
		if(sumSquares == 0.0f)
			return;
		while(--j >= 0) v[j] /= sumSquares;
	}

	public static float L2norm( Matrix m ) {
		int i = m.rows;
		float sumSquares = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sumSquares += m.entry[i][j];
			}
		}
		sumSquares = (float)Math.sqrt(sumSquares);
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[2] = sumSquares;
		return sumSquares;
	}

	public static float L2norm( float [] v ) {
		int j = v.length;
		float sumSquares = 0.0f;
		while(--j >= 0) sumSquares += v[j]*v[j];
		return (float)Math.sqrt(sumSquares);
	}

	public static void L1normalize( Matrix m ) {
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sum += Math.abs(m.entry[i][j]);
			}
		}
		i = m.rows;
		if(sum == 0.0f)
			return;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sum;
			}
		}
	}

	public static void L1normalize( float[] v ) {
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) sum += Math.abs(v[j]);
		j = v.length;
		if(sum == 0.0f)
			return;
		while(--j >= 0) v[j] /= sum;
	}

	public static float L1norm( Matrix m ) {
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sum += Math.abs(m.entry[i][j]);
			}
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[1] = sum;
		return sum;
	}

	public static float L1norm( float[] v ) {
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) sum += Math.abs(v[j]);
		return sum;
	}

	public static void L0normalize( Matrix m ) {
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				if(m.entry[i][j] != Matrix.ZEROF)
					sum += 1.0f;
			}
		}
		i = m.rows;
		if(sum == 0.0f)
			return;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sum;
			}
		}
	}

	public static void L0normalize( float[] v ) {
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) 
			if(v[j] != 0.0f)
				sum += 1.0f;
		j = v.length;
		if(sum == 0.0f)
			return;
		while(--j >= 0) v[j] /= sum;
	}

	public static float L0norm( Matrix m ) {
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				if(m.entry[i][j] != Matrix.ZEROF)
					sum += 1.0f;
			}
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[0] = sum;
		return sum;
	}

	public static float L0norm( float[] v ) {
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) 
			if(v[j] != 0.0f)
				sum += 1.0f;
		return sum;
	}

	// Algorithms

	public static Matrix FloydWarshall( Matrix connectivityMetric, float epsilon, float connected_bound, float disconnected_bound, int strictness, int minmax ) { 
		epsilon = Math.abs(epsilon);
		float connected = connected_bound;
		float disconnected = disconnected_bound;
		int i = connectivityMetric.rows;
		int maxspan = connectivityMetric.shortest;
		Matrix fw = new Matrix(Matrix.NONE, i, connectivityMetric.columns);
		while(--i >= 0) {
			int j = connectivityMetric.columns;
			int max = (minmax == MIN) ? 0 : i;
			while(--j >= max ) {
				float eij = connectivityMetric.entry[i][j];
				if( strictness == 0 && j < maxspan && i < maxspan ) {
					eij += connectivityMetric.entry[j][i];
				}
				if( eij != Matrix.ZEROF && eij >= -epsilon && eij <= epsilon )
					fw.entry[i][j] = connected;
				else
					fw.entry[i][j] = disconnected;
			}
		}
		int k = maxspan;
		while( --k >= 0 ) {      
			i = connectivityMetric.rows;
			while( --i >= 0 ) {
				float current_distance_part = fw.entry[i][k];
				int j = connectivityMetric.columns;
				while( -- j >= 0 ) {
					if(i==j)
						fw.entry[i][j] = 0.0f;
					else {
						float current_distance = fw.entry[i][j];
						float current_distance_whole;
						switch(minmax) {
						case(MIN):
							current_distance_whole = current_distance_part + fw.entry[k][j];
						fw.entry[i][j] = current_distance < current_distance_whole ? current_distance : current_distance_whole;
						break;
						case(MAX):
							current_distance_whole = current_distance_part < fw.entry[k][j] ? current_distance_part : fw.entry[k][j];
						fw.entry[i][j] = current_distance > current_distance_whole ? current_distance : current_distance_whole;
						break;
						}
					}
				}  
			}
		}
		i = connectivityMetric.rows;
		while( --i >= 0 ) {
			int j = connectivityMetric.columns;
			while( --j >= 0 ) 
				fw.entry[i][j] *= fw.entry[i][j];
		}
		return fw;
	}

	public static void PowerIterationDeflate(Matrix m, int eigens_to_compute, double epsilons, int populateKrylovSubspace ) { 
		int iteration = 0, i = -1;
		int maxIteration = m.shortest;
		m.krylovs = new float[populateKrylovSubspace][m.columns];
		while(++i < eigens_to_compute ) {
			Matrix next_eigen_approximation, current_eigen_approximation = new Matrix(Matrix.NONE, 1, m.columns );
			double prev, current = epsilons; float eigenvalue_approximation;
			fill_random( current_eigen_approximation );
			L2normalize( current_eigen_approximation );
			int iIteration = 0;
			do {
				prev = current;
				next_eigen_approximation = new Matrix(current_eigen_approximation);
				current_eigen_approximation = productResult(m, next_eigen_approximation.transpose()).transpose();
				if(iteration<populateKrylovSubspace)
					System.arraycopy(current_eigen_approximation.entry[0],0,m.krylovs[iteration++],0,m.columns);
				eigenvalue_approximation = FrobeniusInnerProduct(next_eigen_approximation, current_eigen_approximation);
				L2normalize(current_eigen_approximation);
				current = FrobeniusInnerProduct(next_eigen_approximation,current_eigen_approximation);
			} while( (current < (1-epsilons) && Math.abs( current/prev ) > (1 + epsilons)) && iIteration++ < maxIteration );
			m.eigenvalue[i] = eigenvalue_approximation;
	                Deflate(m, current_eigen_approximation, eigenvalue_approximation );
			int shortest = current_eigen_approximation.entry[0].length < m.eigenvector[i].length ? current_eigen_approximation.entry[0].length :  m.eigenvector[i].length;
			scalarMultiply(current_eigen_approximation, (float)Math.sqrt(Math.abs(eigenvalue_approximation)));
			System.arraycopy( current_eigen_approximation.entry[0], 0, m.eigenvector[i], 0, shortest );
		}
	}

	public static void show_v (float[] v) {
		int i = v.length;
		while(--i >= 0)
			System.out.print(v[i] + " ");
		System.out.println();
	}

	public static void Deflate( Matrix base, Matrix eigen_vector, float eigen_value ) {
		Matrix r1 = scalarMultiplyResult(eigen_vector, eigen_value);
		Matrix result = productResult(  eigen_vector.transpose(), r1 );
		pairwiseSubtract( base, result );
	}

	public static Matrix DeflateResult( Matrix base, Matrix eigen_vector, float eigen_value ) {
		Matrix result = productResult(  eigen_vector, scalarMultiplyResult(eigen_vector, eigen_value).transpose() );
		result = pairwiseSubtractResult( base, result );
		return result;
	}

	public static void LU_Dolittle_Decomposition( Matrix m ) { 
		int rows = m.rows; int columns = m.columns;
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				for(int k = 0; k < i-1; k++ )
					m.entry[i][j] -= m.entry[i][k]*m.entry[k][j];
			}
			for(int j = i+1; j < rows; j++ ) {
				for(int k = 0; k < i-1; k++ ) 
					m.entry[j][i] -= m.entry[j][k]*m.entry[k][i];
				if(m.entry[i][i] != 0.0f)
					m.entry[j][i] /= m.entry[i][i];
			}
		}
	}

	public static Matrix LU_Dolittle_DecompositionResult( Matrix m ) { 
		int rows = m.rows; int columns = m.columns;
		Matrix lowerupper = new Matrix(Matrix.IDENTITY, rows, columns );
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				lowerupper.entry[i][j] = m.entry[i][j];
				for(int k = 0; k < i-1; k++ )
					lowerupper.entry[i][j] -= lowerupper.entry[i][k]*lowerupper.entry[k][j];
			}
			for(int j = i+1; j < rows; j++ ) {
				lowerupper.entry[j][i] = m.entry[j][i];
				for(int k = 0; k < i-1; k++ ) 
					lowerupper.entry[j][i] -= lowerupper.entry[j][k]*lowerupper.entry[j][i];
				if(lowerupper.entry[i][i] != 0.0f)
					lowerupper.entry[j][i] /= lowerupper.entry[i][i];
			}
		}
		return lowerupper;
	}

	public static Matrix[] LU_Dolittle_DecompositionResultPair( Matrix m ) { 
		int rows = m.rows; int columns = m.columns;
		Matrix lower = new Matrix(Matrix.IDENTITY, rows, columns );
		Matrix upper = new Matrix(Matrix.NONE, rows, columns );
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				upper.entry[i][j] = m.entry[i][j];
				for(int k = 0; k < i-1; k++ )
					upper.entry[i][j] -= lower.entry[i][k]*upper.entry[k][j];
			}
			for(int j = i+1; j < rows; j++ ) {
				lower.entry[j][i] = m.entry[j][i];
				for(int k = 0; k < i-1; k++ ) 
					lower.entry[j][i] -= lower.entry[j][k]*upper.entry[j][i];
				if(lower.entry[i][i] != 0.0f)
					lower.entry[j][i] /= upper.entry[i][i];
			}
		}
		return new Matrix[] { lower, upper };
	}

	public static float diagonalDet( Matrix triangular ) {
		float det = 1.0f;
		int i = triangular.shortest;
		while(--i >= 0 )
			det *= triangular.entry[i][i];
		triangular.det = det;
		return det;
	}

	public static float Determinant( Matrix m ) { 
		Matrix[] lowerupper = LU_Dolittle_DecompositionResultPair(m);
		float detLower = diagonalDet(lowerupper[0]);
		float detUpper = diagonalDet(lowerupper[1]);
		m.det = detLower*detUpper;
		return m.det;
	}

	public static Matrix Inverse( Matrix a ) {
		Matrix m = new Matrix(a);
		Matrix w = new Matrix( Matrix.NONE, m.rows, m.columns );
		Matrix l = new Matrix( Matrix.NONE, m.rows, m.columns );
		Matrix u = new Matrix( Matrix.NONE, m.rows, m.columns );
		for(int k = 0; k < m.shortest; k++ ) {
			l.entry[k][k] = Matrix.IDENTITY;
			for(int i = k+1; i < m.shortest; i++ ) {
				float coeff = Matrix.IDENTITYF; 
				if(m.entry[k][k] != 0.0f)
					coeff = m.entry[i][k]/m.entry[k][k];
				l.entry[i][k] = coeff;
				for(int j = k+1; j < m.shortest; j++ )
					m.entry[i][j] -= coeff*m.entry[k][j];
			}
		}
		for(int j = 0; j < m.shortest; j++ ) {
			for(int i = 0; i < j; i++ ) {
				u.entry[i][j] = m.entry[i][j];
			}
		}
		float b[] = new float[m.shortest];
		float d[] = new float[m.shortest];
		float x[] = new float[m.shortest];
		for(int k = 0; k < m.columns; k++ ) {

			for(int i = 1; i < m.shortest; i++ ) {
				b[k] = 1.0f;
				d[1] = b[1];
				for(int j = 1; j < i-1; j++ ) {
					d[i] = b[i];
					d[i] -= l.entry[i][j]*d[j];
				}
			}
			if(u.entry[m.shortest-1][m.shortest-1] != 0.0f)
				x[m.shortest-1] = d[m.shortest-1]/u.entry[m.shortest-1][m.shortest-1];
			int i = m.shortest;
			while(--i >= 0 ) {
				x[i] = d[i];
				int j = m.shortest;
				while(--j > i ) {
					x[i] = x[i] - u.entry[i][j]*x[j];
				}
				if(u.entry[i][i] != 0.0f) 
					x[i] /= u.entry[i][i]; 
			}
			for(i = 0; i < m.shortest; i++ ) {
				w.entry[i][k] = x[i];
			}
			b[k] = 0.0f;
		}
		return w;
	}

	// Thresholding

	public static void threshold ( Matrix a, float threshold, int directed ) {
		int shortestRow = a.rows;
		int shortestColumn = a.columns ;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(directed==0)
					aij = Math.abs(a.entry[i][j]);
				if(aij < threshold)
					a.entry[i][j] = 0.0f;
				else
					a.entry[i][j] = 1.0f;
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix thresholdResult ( Matrix a, float threshold, int directed ) {
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(directed==0)
					aij = Math.abs(a.entry[i][j]);
				if(aij < threshold)
					result.entry[i][j] = 0.0f;
				else
					result.entry[i][j] = 1.0f;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Pair/Entry-Wise Fuzzy Binary arithmetic
	// a AND 0 = 0
	// a or 0 = a
	// a XOR 0 = a
	// a XOR a = 0
	// a AND b = a*b
	// a OR b = a+b
	// a XOR b = (a/b + b/a)
	// a KOR b = -ab + b/a 
	// NOT a = 0
	// a NAND B = a/b
	// a NOR b = a + 1/b

	public static void AND( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				a.entry[i][j] = a.entry[i][j]*b.entry[i][j];
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ANDResult ( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = a.entry[i][j]*b.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void OR( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				a.entry[i][j] = a.entry[i][j]+b.entry[i][j];
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ORResult ( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = a.entry[i][j]+b.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void XOR( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij == 0) {
					if(bij == 0 )
						a.entry[i][j] = 0.0f;
					else
						a.entry[i][j] = bij;
				}
				else {
					if(bij != 0 )
						a.entry[i][j] = aij/bij + bij/aij;
				}
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix XORResult ( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij == 0) {
					if(bij == 0 )
						result.entry[i][j] = 0.0f;
					else
						result.entry[i][j] = bij;
				}
				else {
					if(bij != 0 )
						result.entry[i][j] = aij/bij + bij/aij;
				}
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void KOR( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij != 0)
					a.entry[i][j] = -aij*bij + bij/aij;
			}
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix KORResult ( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij != 0)
					result.entry[i][j] = -aij*bij + bij/aij;
			}
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void NOT ( Matrix a ) {
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(aij != 0)
					a.entry[i][j] = 1/-aij;
			}
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix NOTResult ( Matrix a ) {
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(aij != 0)
					result.entry[i][j] = 1.0f/-aij;
			}
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// General

	public static Matrix complement ( Matrix m ) {
		float[] mm = m.minmax();
		float min = mm[0];
		float max = mm[1];
		float avg = (min+max)/2.0f;
		Matrix result = new Matrix(Matrix.NONE, m.rows, m.columns);
		int shortestRow = m.rows;
		int shortestColumn = m.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float mij = m.entry[i][j];
				if(mij < avg)
					result.entry[i][j] = max;
				else
					result.entry[i][j] = min;
			}
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Modular arithmetic 

	//Addition

	public static void ModscalarAdd(Matrix accumulator, float addend, float modulus ) {
		if(addend == Matrix.ZEROF || modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] += addend;
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModscalarAddResult(Matrix adder, float addend, float modulus ) {
		Matrix result;
		if(addend == Matrix.ZEROF || modulus == Matrix.ZEROF) {
			result = new Matrix(adder);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = adder.rows;
		int shortestColumn = adder.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = adder.entry[i][j] + addend;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModpairwiseSum(Matrix accumulator, Matrix summand, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows < summand.rows ? accumulator.rows : summand.rows;
		int shortestColumn = accumulator.columns < summand.columns ? accumulator.columns : summand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] += summand.entry[i][j];
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModpairwiseSumResult(Matrix adder, Matrix summand, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			adder.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return adder;
		}int shortestRow = adder.rows < summand.rows ? adder.rows : summand.rows;
		int shortestColumn = adder.columns < summand.columns ? adder.columns : summand.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = adder.entry[i][j] + summand.entry[i][j];
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Subtraction

	public static void ModscalarSubtract(Matrix accumulator, float subtractand, float modulus ) {
		if(modulus == Matrix.ZEROF || subtractand == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] -= subtractand;
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModscalarSubtractResult(Matrix subtractor, float subtractand, float modulus ) {
		Matrix result;
		if(modulus == Matrix.ZEROF || subtractand == Matrix.ZEROF) {
			result = new Matrix(subtractor);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = subtractor.rows;
		int shortestColumn = subtractor.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = subtractor.entry[i][j] - subtractand;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModpairwiseSubtract(Matrix accumulator, Matrix subtractand, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows < subtractand.rows ? accumulator.rows : subtractand.rows;
		int shortestColumn = accumulator.columns < subtractand.columns ? accumulator.columns : subtractand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] -= subtractand.entry[i][j];
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModpairwiseSubtractResult(Matrix subtractor, Matrix subtractand, float modulus ) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(subtractor);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = subtractor.rows < subtractand.rows ? subtractor.rows : subtractand.rows;
		int shortestColumn = subtractor.columns < subtractand.columns ? subtractor.columns : subtractand.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = subtractor.entry[i][j] - subtractand.entry[i][j];
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Multiplication

	public static void ModscalarMultiply(Matrix accumulator, float multiplier, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] *= multiplier;
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModscalarMultiplyResult(Matrix multiplicand, float multiplier, float modulus ) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(multiplicand);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = multiplicand.rows;
		int shortestColumn = multiplicand.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = multiplicand.entry[i][j] * multiplier;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static float ModFrobeniusInnerProduct( Matrix multiplicand, Matrix multiplier, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return modulus;
		}
		int shortestRow = multiplicand.rows < multiplier.rows ? multiplicand.rows : multiplier.rows;
		int shortestColumn = multiplicand.columns < multiplier.columns ? multiplicand.columns : multiplier.columns;
		int i = shortestRow;
		float accumulator = 0.0f;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator += multiplicand.entry[i][j]*multiplier.entry[i][j];
			}  
		}
		return accumulator % modulus;
	}

	public static void ModhadamardProduct(Matrix accumulator, Matrix multiplicand, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows < multiplicand.rows ? accumulator.rows : multiplicand.rows;
		int shortestColumn = accumulator.columns < multiplicand.columns ? accumulator.columns : multiplicand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] *= multiplicand.entry[i][j];
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModhadamardProductResult(Matrix multiplier, Matrix multiplicand, float modulus ) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(multiplicand);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = multiplier.rows < multiplicand.rows ? multiplier.rows : multiplicand.rows;
		int shortestColumn = multiplier.columns < multiplicand.columns ? multiplier.columns : multiplicand.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = multiplier.entry[i][j] * multiplicand.entry[i][j];
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModproductPropagated(Matrix accumulator, Matrix multiplier, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows < multiplier.rows ? accumulator.rows : multiplier.rows;
		int shortestColumn = accumulator.columns < multiplier.columns ? accumulator.columns : multiplier.columns;
		int shortest = shortestRow < shortestColumn ? shortestRow : shortestColumn;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int k = shortest;
				float entrySum = 0.0f;
				while(--k >= 0) {
					entrySum += accumulator.entry[i][k]*multiplier.entry[k][j];
				}
				accumulator.entry[i][j] = entrySum;
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModproductResult(Matrix multiplicand, Matrix multiplier, float modulus ) { 
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(multiplicand);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortest = multiplicand.columns < multiplier.rows ? multiplicand.columns : multiplier.rows;
		int i = multiplicand.rows;
		int j = multiplier.columns;
		result = new Matrix(Matrix.NONE, i, j);
		while(--i >= 0) {
			j = multiplier.columns;
			while(--j >= 0) {
				int k = shortest;
				float entrySum = 0.0f;
				while(--k >= 0) {
					entrySum += multiplicand.entry[i][k]*multiplier.entry[k][j];
				}
				result.entry[i][j] = entrySum;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Division

	public static void ModentrywiseInversion( Matrix divisor, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			divisor.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = divisor.rows;
		int shortestColumn = divisor.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float divisorij = divisor.entry[i][j];
				if(divisorij != Matrix.ZEROF)
					divisor.entry[i][j] = Matrix.IDENTITYF / divisorij;
				else
					divisor.entry[i][j] = Matrix.IDENTITYF;
				divisor.entry[i][j] %= modulus;
			}  
		}
	}

	public static Matrix ModentrywiseInversionResult ( Matrix divisor, float modulus ) {
		int shortestRow = divisor.rows;
		int shortestColumn = divisor.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float divisorij = divisor.entry[i][j];
				if(divisorij != Matrix.ZEROF)
					result.entry[i][j] = Matrix.IDENTITYF / divisorij;
				else
					result.entry[i][j] = Matrix.IDENTITYF;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModscalarDivide(Matrix accumulator, float divisor, float modulus) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		if(divisor == Matrix.ZEROF) {
			Arrays.fill(accumulator.entry, Matrix.IDENTITYF);
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) 
				accumulator.entry[i][j] /= divisor;
		}
	}

	public static Matrix ModscalarDivideResult(Matrix dividend, float divisor, float modulus) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(dividend);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = dividend.rows;
		int shortestColumn = dividend.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor != Matrix.ZEROF)
					result.entry[i][j] = dividend.entry[i][j] / divisor;
				else
					result.entry[i][j] = Matrix.IDENTITYF;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModhadamardDivision(Matrix accumulator, Matrix divisor, float modulus) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows < divisor.rows ? accumulator.rows : divisor.rows;
		int shortestColumn = accumulator.columns < divisor.columns ? accumulator.columns : divisor.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor.entry[i][j] != Matrix.ZEROF)
					accumulator.entry[i][j] /= divisor.entry[i][j];
				else
					accumulator.entry[i][j] = Matrix.ZEROF;
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModhadamardDivisionResult(Matrix dividend, Matrix divisor, float modulus) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(dividend);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = dividend.rows < divisor.rows ? dividend.rows : divisor.rows;
		int shortestColumn = dividend.columns < divisor.columns ? dividend.columns : divisor.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor.entry[i][j] != Matrix.ZEROF)
					result.entry[i][j] = dividend.entry[i][j] / divisor.entry[i][j];
				else
					result.entry[i][j] = Matrix.ZEROF;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Powers

	public static void ModhadamardPower(Matrix accumulator, int power, float modulus) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		if(power == 0) {
			accumulator.make_ones();
			return;  
		}
		if(power == 1)
			return;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int p = power;
				while(--p > 0) {
					accumulator.entry[i][j] *= accumulator.entry[i][j];
					accumulator.entry[i][j] %= modulus;
				}
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModhadamardPowerResult(Matrix multiplier, int power, float modulus) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(multiplier);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = multiplier.rows;
		int shortestColumn = multiplier.columns;
		if(power == 1)
			return new Matrix(multiplier);
		result = new Matrix(Matrix.ONES, shortestRow, shortestColumn);
		if(power == 0)
			return result;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int p = power;
				while(--p >= 0) {
					result.entry[i][j] *= multiplier.entry[i][j];
					result.entry[i][j] %= modulus;
				}
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static int[] Modradix(int num, int base, float modulus) {
		if(modulus == Matrix.ZEROF) {
			return new int[10];
		}
		int b = 1; int b1 = base;
		int rem = 0; int rem1 = num % b1;
		int bit = 0;
		int[] bits = new int[10];
		while(b <= num) {
			bits[bit] = ((rem1 - rem)/b) % (int)modulus;
			b=b1;b1*=base;
			rem=rem1; rem1 = num % b1;
			bit++;
			if(bit >= bits.length) {
				int[] temp = new int[bit<<1];
				System.arraycopy(bits, 0, temp, 0, bit);
				bits = temp;
			}
		}
		return bits;
	}

	public static float ModexponentialTrace( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF)
			return modulus;
		float result = 1.0f;
		float trace = m.trace();
		int itrace = (int)trace;
		float diff = trace - itrace;
		int[] bits = radix(itrace, 2);
		int iteration = 0;
		float lastPow = 2.71828183f;
		int lastPowIndex = 0;
		while(iteration++ < bits.length) {
			int bit = bits[iteration];
			if(bit == 1) {
				int squaresRequired = iteration - lastPowIndex;
				while(--squaresRequired >= 0)
					lastPow = (lastPow * lastPow) % modulus;
				result = (result * lastPow) % modulus;
			}
		}
		float extra = -(float)Math.log(diff);
		if(extra != 0.0f)
			result *= extra;
		return result % modulus;
	}

	public static void mod(Matrix base, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			base.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = base.rows;
		int shortestColumn = base.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				base.entry[i][j] %= modulus;
			}  
		}
		base.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix modResult( Matrix base, float modulus ) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(base);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = base.rows;
		int shortestColumn = base.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0)
				result.entry[i][j] %= modulus;
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModpolysquarePropagated(Matrix base, int times, float modulus) {
		if(modulus == Matrix.ZEROF) {
			base.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		while(--times >= 0) {
			productPropagated(base, base);
			mod(base, modulus);
		}
	}

	public static void ModfastPropagatedMultiproduct(Matrix base, int exponent, float modulus) { 
		if(modulus == Matrix.ZEROF) {
			base.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int[] operations = radix(exponent, 2);
		int lastProductIndex = 0, jmax = operations.length;
		Matrix partialPowers;
		if(operations[0] == 0)
			partialPowers = new Matrix(Matrix.IDENTITY, base.rows, base.columns);
		else
			partialPowers = new Matrix(base);
		for(int j = 1; j < jmax;j++) {
			int bit = operations[j];
			if(bit == 1) {
				ModpolysquarePropagated(partialPowers, j-lastProductIndex, modulus);
				lastProductIndex = j;
				ModproductPropagated(base, partialPowers, modulus);
				mod(base,modulus);
			}
		}
	}

	public static Matrix Modpolysquare(Matrix base, int times, float modulus) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(base);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		while(--times >= 0) {
			base = productResult(base, base);
			base = modResult(base, modulus);
		}
		base.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return base;
	}

	public static Matrix ModfastMultiproduct(Matrix base, int exponent, float modulus) { 
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(base);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int[] operations = radix(exponent, 2);
		int lastProductIndex = 0, jmax = operations.length;
		Matrix partialPowers;
		if(operations[0] == 0)
			partialPowers = new Matrix(Matrix.IDENTITY, base.rows, base.columns);
		else
			partialPowers = new Matrix(base);
		for(int j = 1; j < jmax;j++) {
			int bit = operations[j];
			if(bit == 1) {
				Modpolysquare(partialPowers, j-lastProductIndex, modulus);
				lastProductIndex = j;
				base = ModproductResult(base, partialPowers, modulus);
				base = modResult(base,modulus);
			}
		}
		return base;
	}

	public static Matrix ModsquareRootReturnInverseSquareRoot(Matrix y, float modulus) {
		if(modulus == Matrix.ZEROF) {
			y.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return ModInverse(y, modulus);
		}
		Matrix z = new Matrix(Matrix.IDENTITY, y.rows, y.columns);
		int j = (int)Math.sqrt(y.rows) + 1;
		while(--j >= 0) {
			scalarMultiply(pairwiseSumResult(y,Inverse(z)), 0.5f);
			scalarMultiply(pairwiseSumResult(z,Inverse(y)), 0.5f);
		}
		y.updates &= (Matrix.all_changed^Matrix.entry_changed);
		z.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return z;
	}

	public static Matrix[] ModsquareRootInverseSquareRootResult(Matrix x, float modulus) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(x);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return new Matrix[] {result, ModInverse(result, modulus)};
		}
		Matrix y = new Matrix(x);
		Matrix z = new Matrix(Matrix.IDENTITY, y.rows, y.columns);
		int j = (int)Math.sqrt(y.rows) + 1;
		while(--j >= 0) {
			scalarMultiply(pairwiseSumResult(y,Inverse(z)), 0.5f);
			scalarMultiply(pairwiseSumResult(z,Inverse(y)), 0.5f);
		}
		y.updates &= (Matrix.all_changed^Matrix.entry_changed);
		z.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return new Matrix[] { y, z };
	}

	public static void ModentrywiseSquareRoot(Matrix x, float modulus) {
		if(modulus == Matrix.ZEROF)
			return;
		int shortestRow = x.rows;
		int shortestColumn = x.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float xij = x.entry[i][j] % modulus;
				float sign = xij/Math.abs(xij);
				x.entry[i][j] = (float)( sign*Math.sqrt(Math.abs(xij)) );
			}  
		}
		x.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModentrywiseSquareRootResult(Matrix x, float modulus) {
		if(modulus == Matrix.ZEROF)
			return new Matrix(x);
		int shortestRow = x.rows;
		int shortestColumn = x.columns;
		int i = shortestRow;
		Matrix result = new Matrix(Matrix.NONE, x.rows, x.columns);
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float xij = x.entry[i][j] % modulus;
				float sign = xij/Math.abs(xij);
				result.entry[i][j] = (float)( sign*Math.sqrt(Math.abs(xij)) );
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Logarithm

	public static void ModentrywiseLogarithm( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int shortestRow = m.rows;
		int shortestColumn = m.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float eij = m.entry[i][j];
				float sign = eij/Math.abs(eij);
				m.entry[i][j] = (float)(sign*Math.log(eij));
			}  
		}
		m.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModentrywiseLogarithmResult( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return m;
		}
		int shortestRow = m.rows;
		int shortestColumn = m.columns;
		int i = shortestRow;
		Matrix result = new Matrix(m);
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float eij = m.entry[i][j];
				float sign = eij/Math.abs(eij);
				result.entry[i][j] = (float)(sign*Math.log(eij));
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}


	// Mapping 

	public static void ModscaleToRange(Matrix m, float low, float high, float modulus) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		float[] minmax = m.minmax();
		float newrange = high-low;
		float oldrange = minmax[1]-minmax[0];
		float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
		mod(m,modulus);
		scalarSubtract(m, minmax[0]);
		scalarMultiply(m, range );
		scalarAdd(m, low );

	}

	public static Matrix ModscaleToRangeResult(Matrix m, float low, float high, float modulus) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(m);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		float[] minmax = m.minmax();
		float newrange = high-low;
		float oldrange = minmax[1]-minmax[0];
		float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
		Matrix result = new Matrix(Matrix.NONE, m.rows, m.columns);
		result = ModscalarSubtractResult(m, minmax[0], modulus);
		ModscalarMultiply(result, range, modulus );
		ModscalarAdd(result, low, modulus );
		return result;
	}


	// Random

	public static void Modfill_random(Matrix m, float modulus) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = m.rows;
		while(--i >= 0 ) {
			int j = m.columns;
			while(--j >= 0 )
				m.entry[i][j] = (float)Math.random() % modulus;
		}
	}

	public static void Modfill_random( float[] v, float modulus) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int j = v.length;
		while(--j >= 0 ) v[j] = (float)Math.random() % modulus;
	}


	// Norms & normalization vectors and matrices : 0, 1, 2, inf

	public static void ModLinfnormalize( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = m.rows;
		int j = m.columns;
		if( i == 0 || j == 0)
			return;
		float max = Math.abs(m.entry[0][0]);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				float eij = Math.abs(m.entry[i][j]);
				if(eij > max)
					max = eij;
			}
			j = m.columns;
		}
		i = m.rows;
		while( --i >= 0 ) {
			j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= max;
			}
		}
	}

	public static void ModLinfnormalize( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int j = v.length;
		if(j == 0)
			return;
		float max = Math.abs(v[0]);
		while(--j > 0) {
			float vj = Math.abs(v[j]);
			if(vj > max)
				max = vj;
		}
		j = v.length;
		while(--j >= 0) v[j] /= max;
	}

	public static float ModLinfnorm( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		int i = m.rows;
		int j = m.columns;
		if( i == 0 || j == 0)
			return 0.0f;
		float max = Math.abs(m.entry[0][0]);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				float eij = Math.abs(m.entry[i][j]);
				if(eij > max)
					max = eij;
			}
			j = m.columns;
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[3] = max;
		return max;
	}

	public static float ModLinfnorm( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return modulus;
		}
		int j = v.length;
		if(j == 0)
			return 0.0f;
		float max = Math.abs(v[0]);
		while(--j > 0) {
			float vj = Math.abs(v[j]);
			if(vj > max)
				max = vj;
		}
		return max;
	}

	public static void ModL2normalize( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = m.rows;
		float sumSquares = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sumSquares += m.entry[i][j]*m.entry[i][j];
			}
		}
		sumSquares = (float)Math.sqrt(sumSquares);
		i = m.rows;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sumSquares;
			}
		}
	}

	public static void ModL2normalize( float [] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int j = v.length;
		float sumSquares = 0.0f;
		while(--j >= 0) sumSquares += v[j]*v[j];
		j = v.length;
		sumSquares = (float)Math.sqrt(sumSquares);
		while(--j >= 0) v[j] /= sumSquares;
	}

	public static float ModL2norm( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		int i = m.rows;
		float sumSquares = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sumSquares += (m.entry[i][j]*m.entry[i][j])%modulus;
			}
		}
		sumSquares = (float)Math.sqrt(sumSquares);
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[2] = sumSquares;
		return sumSquares;
	}

	public static float ModL2norm( float [] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return modulus;
		}
		int j = v.length;
		float sumSquares = 0.0f;
		while(--j >= 0) sumSquares += (v[j]*v[j])%modulus;
		return (float)Math.sqrt(sumSquares);
	}

	public static void ModL1normalize( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sum += Math.abs(m.entry[i][j]);
			}
		}
		i = m.rows;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sum;
			}
		}
	}

	public static void ModL1normalize( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) sum += Math.abs(v[j]);
		j = v.length;
		while(--j >= 0) v[j] /= sum;
	}

	public static float ModL1norm( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sum += Math.abs(m.entry[i][j]);
			}
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[1] = sum;
		return sum;
	}

	public static float ModL1norm( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return modulus;
		}
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) sum += Math.abs(v[j]);
		return sum;
	}

	public static void ModL0normalize( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				if(m.entry[i][j] != Matrix.ZEROF)
					sum += 1.0f;
			}
		}
		i = m.rows;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sum;
			}
		}
	}

	public static void ModL0normalize( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) 
			if(v[j] != 0.0f)
				sum += 1.0f;
		j = v.length;
		while(--j >= 0) v[j] /= sum;
	}

	public static float ModL0norm( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				if(m.entry[i][j] != Matrix.ZEROF)
					sum += 1.0f;
			}
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[0] = sum;
		return sum;
	}

	public static float ModL0norm( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return modulus;
		}
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) 
			if(v[j] != 0.0f)
				sum += 1.0f;
		return sum;
	}

	// Algorithms

	public static Matrix ModFloydWarshall( Matrix base, float epsilon, float connected_bound, float disconnected_bound, int strictness, int minmax, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(base);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		epsilon = Math.abs(epsilon);
		float connected = connected_bound;
		float disconnected = disconnected_bound;
		int i = base.rows;
		int maxspan = base.shortest;
		Matrix fw = new Matrix(Matrix.NONE, i, base.columns);
		while(--i >= 0) {
			int j = base.columns;
			while(--j >= 0 ) {
				float eij = base.entry[i][j];
				if( strictness == 0 && j < maxspan && i < maxspan ) {
					eij += base.entry[j][i];
				}
				if( eij > -epsilon && eij < epsilon )
					fw.entry[i][j] = connected;
				else
					fw.entry[i][j] = disconnected;
			}
		}
		int k = maxspan;
		while( --k >= 0 ) {      
			i = base.rows;
			while( --i >= 0 ) {
				float current_distance_part = fw.entry[i][k];
				int j = base.columns;
				while( -- j >= 0 ) {
					if(i==j)
						fw.entry[i][j] = 0.0f;
					else {
						float current_distance = fw.entry[i][j];
						float current_distance_whole = current_distance_part + fw.entry[k][j];
						switch(minmax) {
						case(MIN):
							fw.entry[i][j] = current_distance < current_distance_whole ? current_distance : current_distance_whole;
						break;
						case(MAX):
							fw.entry[i][j] = current_distance > current_distance_whole ? current_distance : current_distance_whole;
							break;
						}
						fw.entry[i][j] %= modulus;
					}
				}  
			}
		}
		i = base.rows;
		while( --i >= 0 ) {
			int j = base.columns;
			while( --j >= 0 ) 
				fw.entry[i][j] *= fw.entry[i][j];
		}
		return fw;
	}

	public static void ModPowerIterationDeflate(Matrix m, int eigens_to_compute, double epsilons, int populateKrylovSubspace, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int iteration = 0;
		m.krylovs = new float[populateKrylovSubspace][m.columns];
		while(--eigens_to_compute >= 0) {
			Matrix next_eigen_approximation, current_eigen_approximation = new Matrix(Matrix.NONE, 1, m.columns );
			double prev, current = epsilons; float eigenvalue_approximation;
			fill_random( current_eigen_approximation );
			L2normalize( current_eigen_approximation );
			do {
				prev = current;
				next_eigen_approximation = new Matrix(current_eigen_approximation);
				current_eigen_approximation = ModproductResult(m, next_eigen_approximation, modulus);
				if(iteration<populateKrylovSubspace)
					System.arraycopy(current_eigen_approximation.entry[0],0,m.krylovs[iteration++],0,m.columns);
				eigenvalue_approximation = ModFrobeniusInnerProduct(next_eigen_approximation, current_eigen_approximation, modulus);
				ModL2normalize(current_eigen_approximation, modulus);
				current = ModFrobeniusInnerProduct(next_eigen_approximation,current_eigen_approximation, modulus);
			} while( (current < (1-epsilons) && Math.abs( current/prev ) > (1 + epsilons)) );
			ModscalarMultiply(current_eigen_approximation, (float)Math.sqrt(Math.abs(eigenvalue_approximation)), modulus);
			m.eigenvalue[eigens_to_compute] = eigenvalue_approximation;
			System.arraycopy( current_eigen_approximation.entry[0], 0,  m.eigenvector[eigens_to_compute],0, m.columns );
			if(eigens_to_compute > 0)
				m = ModDeflateResult(m, current_eigen_approximation, eigenvalue_approximation, modulus );
		}
	}

	public static Matrix ModDeflateResult( Matrix base, Matrix eigen_vector, float eigen_value, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(base);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		Matrix result = ModproductResult(  eigen_vector, scalarMultiplyResult(eigen_vector, eigen_value).transpose(), modulus );
		result = ModpairwiseSubtractResult( base, result, modulus );
		return result;
	}

	public static void ModLU_Dolittle_Decomposition( Matrix m, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int rows = m.rows; int columns = m.columns;
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				for(int k = 0; k < i-1; k++ ) {
					m.entry[i][j] -= m.entry[i][k]*m.entry[k][j];
					m.entry[i][j] %= modulus;
				}
			}
			for(int j = i+1; j < rows; j++ ) {
				for(int k = 0; k < i-1; k++ ) 
					m.entry[j][i] -= m.entry[j][k]*m.entry[k][i];
				m.entry[j][i] /= m.entry[i][i];
				m.entry[j][i] %= modulus;
			}
		}
	}

	public static Matrix ModLU_Dolittle_DecompositionResult( Matrix m, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(m);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int rows = m.rows; int columns = m.columns;
		Matrix lowerupper = new Matrix(Matrix.IDENTITY, rows, columns );
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				lowerupper.entry[i][j] = m.entry[i][j];
				for(int k = 0; k < i-1; k++ ) {
					lowerupper.entry[i][j] -= lowerupper.entry[i][k]*lowerupper.entry[k][j];
					lowerupper.entry[i][j] %= modulus;
				}
			}
			for(int j = i+1; j < rows; j++ ) {
				lowerupper.entry[j][i] = m.entry[j][i];
				for(int k = 0; k < i-1; k++ ) 
					lowerupper.entry[j][i] -= lowerupper.entry[j][k]*lowerupper.entry[j][i];
				lowerupper.entry[j][i] /= lowerupper.entry[i][i];
				lowerupper.entry[j][i] %= modulus;
			}
		}
		return lowerupper;
	}

	public static Matrix[] ModLU_Dolittle_DecompositionResultPair( Matrix m, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(m);
			Matrix id = new Matrix(Matrix.IDENTITY, m.rows, m.columns );
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			id.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return new Matrix[] { result, id };
		}
		int rows = m.rows; int columns = m.columns;
		Matrix lower = new Matrix(Matrix.IDENTITY, rows, columns );
		Matrix upper = new Matrix(Matrix.NONE, rows, columns );
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				upper.entry[i][j] = m.entry[i][j];
				for(int k = 0; k < i-1; k++ ) {
					upper.entry[i][j] -= lower.entry[i][k]*upper.entry[k][j];
					upper.entry[i][j] %= modulus;
				}
			}
			for(int j = i+1; j < rows; j++ ) {
				lower.entry[j][i] = m.entry[j][i];
				for(int k = 0; k < i-1; k++ ) 
					lower.entry[j][i] -= lower.entry[j][k]*upper.entry[j][i];
				lower.entry[j][i] /= upper.entry[i][i];
				lower.entry[j][i] %= modulus;
			}
		}
		return new Matrix[] { lower, upper };
	}

	public static float ModdiagonalDet( Matrix triangular, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			triangular.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		float det = 1.0f;
		int i = triangular.shortest;
		while(--i >= 0 )
			det = (det * triangular.entry[i][i]) % modulus;
		triangular.det = det;
		return det;
	}

	public static float ModDeterminant( Matrix m, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		Matrix[] lowerupper = ModLU_Dolittle_DecompositionResultPair(m, modulus);
		float detLower = ModdiagonalDet(lowerupper[0], modulus);
		float detUpper = ModdiagonalDet(lowerupper[1], modulus);
		m.det = (detLower*detUpper) % modulus;
		return m.det;
	}

	public static Matrix ModInverse( Matrix a, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		Matrix m = new Matrix(a);
		Matrix w = new Matrix( Matrix.NONE, m.rows, m.columns );
		Matrix l = new Matrix( Matrix.NONE, m.rows, m.columns );
		Matrix u = new Matrix( Matrix.NONE, m.rows, m.columns );
		for(int k = 0; k < m.shortest; k++ ) {
			l.entry[k][k] = Matrix.IDENTITY;
			for(int i = k+1; i < m.shortest; i++ ) {
				float coeff = m.entry[i][k]/m.entry[k][k];
				l.entry[i][k] = coeff;
				for(int j = k+1; j < m.shortest; j++ ) {
					m.entry[i][j] -= coeff*m.entry[k][j];
					m.entry[i][j] %= modulus;
				}
			}
		}
		for(int j = 0; j < m.shortest; j++ ) {
			for(int i = 0; i < j; i++ ) {
				u.entry[i][j] = m.entry[i][j] % modulus;
			}
		}
		float b[] = new float[m.shortest];
		float d[] = new float[m.shortest];
		float x[] = new float[m.shortest];
		for(int k = 0; k < m.columns; k++ ) {
			d[1] = b[1];
			b[k] = 1.0f;
			for(int i = 1; i < m.shortest; i++ ) {
				for(int j = 1; j < i-1; j++ ) {
					d[i] = b[i];
					d[i] -= l.entry[i][j]*d[j];
					d[i] %= modulus;
				}
			}
			x[m.shortest-1] = d[m.shortest-1]/u.entry[m.shortest-1][m.shortest-1];
			for(int i = m.shortest; i >= 0; i-- ) {
				x[i] = d[i];
				for(int j = m.shortest; j > i; j-- ) {
					x[i] = x[i] - u.entry[i][j]*x[j];
					x[i] %= modulus;
				}
				x[i] /= u.entry[i][i]; 
			}
			for(int i = 0; i < m.shortest; i++ ) {
				w.entry[i][k] = x[i];
				w.entry[i][k] %= modulus;
			}
			b[k] = 0.0f;
		}
		return w;
	}

	// Thresholding

	public static void Modthreshold ( Matrix a, float threshold, int directed, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows;
		int shortestColumn = a.columns ;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(directed==0)
					aij = Math.abs(a.entry[i][j]);
				if((aij % modulus ) < threshold)
					a.entry[i][j] = 0.0f;
				else
					a.entry[i][j] = 1.0f;
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModthresholdResult ( Matrix a, float threshold, int directed, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(directed==0)
					aij = Math.abs(a.entry[i][j]);
				if((aij % modulus) < threshold)
					result.entry[i][j] = 0.0f;
				else
					result.entry[i][j] = 1.0f;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Pair/Entry-Wise Fuzzy Binary arithmetic
	// a AND 0 = 0
	// a or 0 = a
	// a XOR 0 = a
	// a XOR a = 0
	// a AND b = a*b
	// a OR b = a+b
	// a XOR b = (a/b + b/a)
	// a KOR b = -ab + b/a 
	// NOT a = 0
	// a NAND B = a/b
	// a NOR b = a + 1/b

	public static void ModAND( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				a.entry[i][j] = a.entry[i][j]*b.entry[i][j];
				a.entry[i][j] %= modulus;
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModANDResult ( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = a.entry[i][j]*b.entry[i][j];
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModOR( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				a.entry[i][j] = (a.entry[i][j]+b.entry[i][j]) % modulus;
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModORResult ( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = a.entry[i][j]+b.entry[i][j];
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModXOR( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij == 0) {
					if(bij == 0 )
						a.entry[i][j] = 0.0f;
					else
						a.entry[i][j] = bij;
				}
				else {
					if(bij != 0 )
						a.entry[i][j] = aij/bij + bij/aij;
				}
				a.entry[i][j] %= modulus;
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModXORResult ( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij == 0) {
					if(bij == 0 )
						result.entry[i][j] = 0.0f;
					else
						result.entry[i][j] = bij;
				}
				else {
					if(bij != 0 )
						result.entry[i][j] = aij/bij + bij/aij;
				}
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModKOR( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij != 0) {
					a.entry[i][j] = -aij*bij + bij/aij;
					a.entry[i][j] %= modulus;
				}
			}
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModKORResult ( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij != 0) {
					result.entry[i][j] = -aij*bij + bij/aij;
					result.entry[i][j] %= modulus;
				}
			}
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModNOT ( Matrix a, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(aij != 0) {
					a.entry[i][j] = 1/-aij;
					a.entry[i][j] %= modulus;
				}
			}
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModNOTResult ( Matrix a, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(aij != 0) {
					result.entry[i][j] = 1.0f/-aij;
					result.entry[i][j] %= modulus;
				}
			}
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

}

import processing.core.*; 

import java.util.ArrayList; 

public class AspektWorld extends PApplet {
	public static final int BIG_STROKE = 3;
	public static final int MIDDLE_STROKE = 1 ;
	public static final int SMALL_STROKE = 1;
	// private static final long serialVersionUID = 61803L;
	public final static int FORWARD = 1;
	public final static int BACKWARD = 2;
	public final static int RIGHTWARD = 4;
	public final static int LEFTWARD = 8;
	public final static int UPWARD = 16;
	public final static int DOWNWARD = 32;
	static public void main(String args[]) {
		PApplet.main(new String[] { "--bgcolor=#D4D0C8", "AspektWorld" });
	}
	public int direction = 0;
	public int last_direction = 0;
	public float lmx,lmy,dx,dy;
	public final static float MAX_SPEED = 60.0f;
	public final static float MIN_SPEED = 1.0f;
	public float speed = 60.0f;
	public float wd = 1600;
	public float hh = 900;
	public float dp = 1000;
	public PVector forward_heading;
	public PVector back_heading,left_heading,right_heading,up_heading,down_heading;
	public PVector straight_up, straight_left, straight_down, straight_right, straight_forward, straight_back;
	public int appletMode=1;
	public int diameter = 100;
	public float scalar;
	public int longest;
	public double ttheta;
	// public Rotation r;
	public int double_click;
	public float sbxs,sbys;
	public ArrayList<Line> lines;
	public int edgeColor = 1;
	public int mouseLookOn = 1;
	public int verticesOn = 0;
        public int edgesOn = 1;
	public float edgesOnC = 1;
	public float rotations = 0.0f;
	public int rotOn = 1;
	public int textOn = 1;
	public int fixedOn = 1;
        public int animateOn = 1;
        public int faster = 0;


	public Graph.Vertex hoverNode;
	public PMatrix3D currCameraMatrix;

	public PGraphics3D g3;
	public int in_a_contro2222l = 0;
	public int basics_panel_active = 0;
	public int matrix_panel_active = 0;
	public int ev_panel_active = 0;
	public int instructions_panel_active = 0;
	public int ev1_sz, ev2_sz, ev3_sz, ev4_sz, ev5_sz;
	public int ev_sz = 138, ev_sz_actv = 150;
	public int mpa_sz, mpl_sz, mpd_sz;
	public int m_sz = 100, m_sz_actv = 168;
	public int mpa_sd, mpl_sd, mpd_sd;
	public float m_actv_delta = (m_sz_actv-m_sz)/2.0f;
	public int counter = 0;
	public PFont face; 
	public float th_constant;
	public Visionary cs;
	public int droplets;
	public Rain[] drops;

	PVector pos;


	public Graph g1;
	public int layoutOn = 1;

	//public void attach_graph_to_local_coordinate_master() {
	//	int i = g1.v.length;
	//	while(--i>=0)
	//		g1.v[i].posa = pos;
	//	i = g1.e.length;
	//}



	@Override
	public void draw() {
                if(animateOn==1) {
                  if(faster==1)
                    edgesOnC += 0.3;
                  else
                    edgesOnC += 0.1;
                  if(edgesOnC >= 64.0f)
                    edgesOnC = 0.0f;
                  edgesOn = (int)edgesOnC;
                }
		//in_a_control = 0;
		background(0);
		pushMatrix();
		translate(pos.x,pos.y,pos.z);
		//draw_grid();
		g1.draw();
		//rotateY(rotations);
		//if(rotOn==1)
		//  rotations-=0.15f;
		//int i=droplets;        
		//while(--i>=0)
		//  drops[i].draw();
		//if(rotations <= -2*PI)
		//  rotations = 0.0f;
		popMatrix();
		//currCameraMatrix = new PMatrix3D(g3.camera);
		//camera();
		//if(matrix_panel_active != 0)
		//	draw_matrix_panel();
		//draw_controls();
		//g3.camera = currCameraMatrix;
		update_heading();
	}

	public void draw_grid() {

		int linum = lines.size();
		for(int i = 0; i < linum;i++) {
			Line lin = lines.get(i);
			lin.draw();
		}  
	}


	public void dvals() {
		pos = new PVector(3600,300,1000);
		lmx = 0.0f; lmy = 0.0f;
		//int boxnum = boxes.size();
		forward_heading = new PVector(1000,-300,500);
		straight_up = new PVector(0,1,0);
		straight_down = new PVector(0,-1,0);
		straight_left = new PVector(-1,0,0);
		straight_right = new PVector(1,0,0);
		straight_forward = new PVector(0,0,1);
		straight_back = new PVector(0,0,-1);
		make_normals();
		//cam.lookAt(forward_heading.x, forward_heading.y, forward_heading.z, 1.0f);

		camera(0,0,0, forward_heading.x, forward_heading.y, forward_heading.z, 0.0f, 1.0f, 0.0f);
	}


	public void input() {
		float mx = mouseX;
		float my = mouseY;
		//if(result < height*height/4 ) {
		if(mouseLookOn==1)
			cursor(CROSS);  
		else
			cursor(HAND);
		mx = map(mx, 0, width, -1, 1);
		my = map(my, 0, height, -1, 1); 
		//}
		//else {
		// cursor(HAND);
		// return; 
		//}
		float xang = mx * Math.abs(mx) / (23*PI);
		float yang = my * Math.abs(my) / (12*PI);
		if(mouseLookOn==1) {
			right_heading.mult(speed*sin(xang));
			forward_heading.mult(cos(xang));
			forward_heading.add(right_heading);
			straight_up.mult(speed*sin(yang));
			forward_heading.mult(cos(yang));
			forward_heading.add(straight_up);
		}
		//cam.lookAt(forward_heading.x, forward_heading.y, forward_heading.z, 1.0f);
		camera(0,0,0, forward_heading.x, forward_heading.y, forward_heading.z, 0.0f, 1.0f, 0.0f);
		forward_heading.normalize();
		update_normals();
		forward_heading.mult(speed);
	}


	@Override
	public void keyPressed(){
		if(key != CODED)
			switch(key) {
			case('w'):case('W'):direction |=FORWARD; last_direction = FORWARD; break;
			case('d'):case('D'):direction |=RIGHTWARD; last_direction = RIGHTWARD; break;
			case('s'):case('S'):direction |=BACKWARD; last_direction = BACKWARD; break;
			case('a'):case('A'):direction |=LEFTWARD; last_direction = LEFTWARD; break;
			case(' '):direction |=UPWARD; last_direction = UPWARD; break;
			case('c'):case('C'):direction |=DOWNWARD; last_direction = DOWNWARD; break;

			case('t'):case('T'):
				textOn ^= 1;
			break;
                        case('P'):
                            faster ^= 1;
                            break;
                       
                        case('p'):
                            animateOn ^=1;
                            edgesOnC = (float)edgesOn;
                            break;
                        case('['):
                                edgesOn -= 1;
                                if(edgesOn <= -1)
                                  edgesOn = 63;
                                break;
			
                        case(']'):
                                edgesOn += 1;
                                if(edgesOn >= 64)
                                  edgesOn = 0;
                                break;
			case('3'):
                                //edgesOn = 1;
				layoutOn ^= 1;
			if(layoutOn == 0) 
				cs.randomizeLayoutColorSize();
			else
				cs.squaresTrianglesCirclesLayout();
			break;
			case('`'):case('~'):
				mouseLookOn ^= 1;
			break;
			case('v'):case('V'):
				verticesOn ^= 1;
			break;
			case('r'):case('R'):
				rotOn ^= 1;
			}
		else {
			switch(keyCode) {
			case(SHIFT):
				if(speed > MIN_SPEED) {
					speed = MIN_SPEED; 
					break;
				}
			}
		}
	}

	@Override
	public void keyReleased(){  
		if(key != CODED)
			switch(key) {
			case('w'):case('W'):direction ^=FORWARD; last_direction = FORWARD; break;
			case('d'):case('D'):direction ^=RIGHTWARD; last_direction = RIGHTWARD; break;
			case('s'):case('S'):direction ^=BACKWARD; last_direction = BACKWARD; break;
			case('a'):case('A'):direction ^=LEFTWARD; last_direction = LEFTWARD; break;
			case(' '):direction ^=UPWARD; last_direction = UPWARD; break;
			case('c'):case('C'):direction ^=DOWNWARD; last_direction = DOWNWARD; break;

			}  
		else
			switch(keyCode) {
			case(SHIFT):
				if(speed < MAX_SPEED) {
					speed = MAX_SPEED; 
					break;
				}
			}
	}

	public void make_grid() {
		float spacing = 100;
		lines = new ArrayList((int)(2*(wd+dp)/spacing + 1));
		stroke(200,255,255,255);
		for(float i = -wd*2; i <= wd*2; i+=spacing) {
			Line nl = new Line(this, i,0,-2*wd,i,0,2*wd);
			nl.posa = pos;
			lines.add(nl);
		}
		for(float i = -wd*2; i <= wd*2; i+= spacing) {
			Line nl = new Line(this, -2*wd,0,i,2*wd,0,i);
			nl.posa = pos;
			lines.add(nl);
		}  
		stroke(0,0,0,255);
	}

	public void make_normals() {
		back_heading = new PVector(-forward_heading.x, -forward_heading.y, -forward_heading.z);
		up_heading = new PVector(0,speed,0);
		down_heading = new PVector(0,-speed,0);
		right_heading = forward_heading.cross(up_heading);
		left_heading = forward_heading.cross(down_heading);
	}



	@Override
	public void mouseClicked() {
		if(mouseEvent.getClickCount() == 1) { 
			fixedOn ^= 1;
			cs.randomGraph();
			cs.randomizeLayoutColorSize();
			g1 = cs.g1;
			if(layoutOn==1)
				cs.squaresTrianglesCirclesLayout();
			if(fixedOn == 1) {
				//attach_graph_to_local_coordinate_master();
				make_normals();
				make_grid();
				update_heading();
			}
		}
		else if(mouseEvent.getClickCount() == 3) {
			dvals();
			
		}



	}


	@Override
	public void setup() {
		textSize(15);
		th_constant = textAscent();
		//textFont(face);
		int xw = (screenWidth*15)/16;
		int yh = (screenHeight*7)/8;
		size(900,568,P3D);
		sbxs = width/50;
		sbys = -height/50;
		rectMode(CENTER);
		colorMode(HSB,360,255,255,255);
		longest = width < height ? 1 : 0;
		scalar = width < height ? width : height;
		//tthea = tan(scalar);
		g3 = (PGraphics3D)g;
		background(0);
		cs = new Visionary(this);
		cs.randomGraph();
		cs.randomizeLayoutColorSize();
		g1 = cs.g1;
		cs.squaresTrianglesCirclesLayout();
		//drops = cs.makeRain();
		//droplets = drops.length;
		dvals();
		make_normals();
		make_grid();
		update_heading();
		//smooth();
		// not for applet
	}

	public void setup_matrix_panel() {
		mpa_sz = m_sz;
		mpa_sd = 0;
		mpl_sz = m_sz;
		mpl_sd = 0;
		mpd_sz = m_sz_actv;
		mpd_sd = 1;
	}

	public void update_heading() {
		input();
		update_v();
	}
	public void update_normals() {
		back_heading.set(-forward_heading.x, -forward_heading.y, -forward_heading.z); 
		up_heading.set(0,speed,0);
		down_heading.set(0,-speed,0);
		straight_left.set(-1,0,0);
		straight_right.set(1,0,0);
		straight_forward.set(0,0,1);
		straight_back.set(0,0,-1);
		straight_up.set(0,1,0);
		straight_down.set(0,-1,0);
		right_heading = forward_heading.cross(straight_up);
		left_heading = forward_heading.cross(straight_down);
		back_heading.mult(speed);
		right_heading.mult(speed);
		left_heading.mult(speed);
	}



	public void update_selected() {
		int ns = g1.v.length;
		for(int i = 0; i < ns;i++) {
			Graph.Vertex n = g1.v[i];
			if(n == null)
				continue;
		}
	}

	public void update_v() {
		switch(direction) {
		case FORWARD:pos.sub(forward_heading); break;
		case BACKWARD:pos.sub(back_heading); break;
		case LEFTWARD:pos.sub(left_heading); break;
		case RIGHTWARD:pos.sub(right_heading); break;
		case UPWARD:pos.add(up_heading); break;
		case DOWNWARD:pos.add(down_heading); break;
		case FORWARD|LEFTWARD:pos.sub(forward_heading);pos.sub(left_heading); break;
		case FORWARD|RIGHTWARD:pos.sub(forward_heading);pos.sub(right_heading);  break;
		case BACKWARD|LEFTWARD:pos.sub(back_heading);pos.sub(left_heading);  break;
		case BACKWARD|RIGHTWARD:pos.sub(back_heading);pos.sub(right_heading);  break;
		case UPWARD|FORWARD:pos.sub(forward_heading);pos.add(up_heading);  break;
		case UPWARD|BACKWARD:pos.sub(back_heading);pos.add(up_heading);  break;
		case UPWARD|LEFTWARD:pos.add(up_heading);pos.sub(left_heading);  break;
		case UPWARD|RIGHTWARD:pos.add(up_heading);pos.sub(right_heading);  break;
		case DOWNWARD|FORWARD:pos.sub(forward_heading);pos.add(down_heading);  break;
		case DOWNWARD|BACKWARD:pos.sub(back_heading);pos.add(down_heading);  break;
		case DOWNWARD|LEFTWARD:pos.add(down_heading);pos.sub(left_heading);  break;
		case DOWNWARD|RIGHTWARD:pos.add(down_heading);pos.sub(right_heading);  break;
		case UPWARD|FORWARD|RIGHTWARD:pos.sub(forward_heading);pos.add(up_heading);pos.sub(right_heading);break;
		case UPWARD|FORWARD|LEFTWARD:pos.sub(forward_heading);pos.add(up_heading);pos.sub(left_heading);  break;
		case UPWARD|BACKWARD|RIGHTWARD:pos.sub(back_heading);pos.add(up_heading);pos.sub(right_heading);  break;
		case UPWARD|BACKWARD|LEFTWARD:pos.sub(back_heading);pos.add(up_heading);pos.sub(left_heading);  break;
		case DOWNWARD|FORWARD|LEFTWARD:pos.sub(forward_heading);pos.add(down_heading);pos.sub(left_heading);  break;
		case DOWNWARD|FORWARD|RIGHTWARD:pos.sub(forward_heading);pos.add(down_heading);pos.sub(right_heading);  break;
		case DOWNWARD|BACKWARD|LEFTWARD:pos.sub(back_heading);pos.add(down_heading);pos.sub(left_heading); break;
		case DOWNWARD|BACKWARD|RIGHTWARD:pos.sub(back_heading);pos.add(down_heading);pos.sub(right_heading);  break;
		case LEFTWARD|RIGHTWARD: 
			if(last_direction==LEFTWARD)
				pos.sub(left_heading);
			else
				pos.sub(right_heading);  
			break;
		case FORWARD|BACKWARD: 
			if(last_direction==FORWARD)
				pos.sub(forward_heading);
			else
				pos.sub(back_heading);
			break;
		case UPWARD|DOWNWARD:
			if(last_direction==UPWARD)
				pos.add(up_heading);
			else
				pos.add(down_heading);
			break;
		} 
		if(pos.y < 300)
			pos.y = 300;
	}
}

import java.util.HashMap;
import java.util.HashSet;
import java.util.ArrayList;
import java.util.Iterator;
import java.util.Set;

import processing.core.PVector;

class Rain extends Graph {
	public int resetCount;
	public PVector dpos;
	public AspektWorld aw;
	public Rain (AspektWorld awa) {
		super(awa);
		aw = awa;
		reset();
	}

	public void reset() {
		resetCount = 120+(int)aw.random(120);
		int vn = 3 + (int)aw.random(4);
		int i = 0;
		Graph.Vertex vl = this.new Vertex();
		vl.color = aw.color(200,200,200,200);
		vl.radius = 2;
		addVertex(vl);
		ArrayList<Graph.Vertex> knownVertices = new ArrayList<Graph.Vertex>();
		knownVertices.add(vl);
		dpos = new PVector(aw.random(4*aw.width)-2*aw.width,18000+aw.random(4000),aw.random(4*aw.width)-2*aw.width);
		while(i < vn) {
			Graph.Vertex vr = this.new Vertex();
			addVertex(vr);
			vr.color = aw.color(200,200,200,200);
			vr.radius = 2;
			vr.pos = new PVector(aw.random(300)-150,aw.random(300)-150,aw.random(300)-150);
			addEdge(new Edge(vl,vr));
			int randomKV = (int)aw.random(i);
			vl = (Graph.Vertex)knownVertices.get(randomKV);
			knownVertices.add(vr);
			i++;
		}
		crystallizeGraph();  
	}

	public void draw() {
		aw.pushMatrix();
		aw.translate(-dpos.x,-dpos.y,-dpos.z);
		super.draw();
		aw.popMatrix();
		if(dpos.y > 111) 
			dpos.y -= 10*v.length*v.length;
		else
			resetCount--;
		if(resetCount <= 0)
			reset();
	}



}

// The local image of the global concept
public class Graph { // implements Comparable<Graph> .... just joking
         
	public class MatrixTypeIndex {
		public static final int ADJACENCY = 0;
		public static final int DEGREE = 1;
		public static final int LAPLACIAN = 2;
		public static final int COLLOCATION = 3;
		public static final int MATRIX_FIELD_COUNT = 4;
	}
	public class EdgeBattle { 
		public static final int ESCHER = -8;
		public static final int ONEFROMTHREE = -5;
		public static final int TOTALFAIL = -2;
		public static final int DEPENDSWHICHWAYYOULOOKATIT = 1;
		public static final int PARALLEL = 4;
	}
	public class Vertex {
		public int id; public PVector pos, posa; public int color; public float radius; public int hover;
		public Vertex() {
			id = 0; 
			hover = 0;
			pos = new PVector(0,0,0);
		}
		public void draw() {
			aspektWorld.pushMatrix();
			aspektWorld.fill(color);
			aspektWorld.stroke(color);
			aspektWorld.strokeWeight(aspektWorld.SMALL_STROKE);
			aspektWorld.translate(pos.x,pos.y,pos.z);
			//aspektWorld.rotateY(aspektWorld.rotations);
			aspektWorld.box(radius);
			aspektWorld.popMatrix();
		}
	}
	public class Edge {
		public final int edgeColorB = aspektWorld.color(141,250,250,250); 
		public Vertex a; public Vertex b; public int[] strength; 
		public Edge( Vertex av, Vertex bv ) {
			a = av; b = bv; strength = new int[6];
		}
		public void draw() {
			if(a.hover==1) {
				aspektWorld.strokeWeight(aspektWorld.BIG_STROKE);
				aspektWorld.stroke(edgeColorB);
			}
			else {
				aspektWorld.strokeWeight(aspektWorld.SMALL_STROKE);
				aspektWorld.stroke(a.color^((int)aspektWorld.alpha(a.color)<<24)^(230<<24));
			}
			aspektWorld.line(a.pos.x,a.pos.y,a.pos.z,
					b.pos.x,b.pos.y,b.pos.z);
		}
	}
	public HashMap<Vertex, Set<Edge>> edgesByVertexA, edgesByVertexB;
	public HashMap<Vertex,Integer> numbering;
	public HashMap<Edge,Integer> edgeNumbering;
	public Matrix[] matrix;
	public Vertex[] v;
	public Edge[] e;
	public float volume;
	private final AspektWorld aspektWorld;
        public final int red,yellow,green,blue,gray;

	public Graph(AspektWorld aspektWorld) {
		this.aspektWorld = aspektWorld;
red = aspektWorld.color(0,200,200,200);
yellow = aspektWorld.color(60,200,200,200);
green = aspektWorld.color(119,255,250,200);
blue = aspektWorld.color(200,200,200,200);
gray = aspektWorld.color(0,0,200,100);
        
		edgesByVertexA = new HashMap<Vertex, Set<Edge>>();
		edgesByVertexB = new HashMap<Vertex, Set<Edge>>();
		matrix = new Matrix[MatrixTypeIndex.MATRIX_FIELD_COUNT];
	}

	public HashMap<Graph.Vertex,Integer> numbering() {
		if(numbering != null)
			return numbering;
		int vertices = edgesByVertexA.size();
		Set<Graph.Vertex> verticeSet = edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		numbering = new HashMap<Graph.Vertex,Integer> ();
		int i = vertices;
		while(--i >= 0)
			numbering.put(vertexIterator.next(), new Integer(i));
		return numbering;
	}

	public HashMap<Graph.Edge, Integer> edgeNumbering() {
		if(edgeNumbering != null)
			return edgeNumbering;
		Set<Graph.Vertex> verticeSet = edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		edgeNumbering = new HashMap<Edge,Integer> ();
		int edges = 0;
		while(vertexIterator.hasNext()) {
			Graph.Vertex stupidFuckFace = vertexIterator.next();
			Set<Graph.Edge> fuckingMoron = edgesByVertexA.get(stupidFuckFace);
			Iterator<Graph.Edge> edgeIterator = fuckingMoron.iterator();
			while(edgeIterator.hasNext()) {
				Edge e = edgeIterator.next();
				if(!edgeNumbering.containsKey(e))
					edgeNumbering.put(e,edges++);
			}
		}
		return edgeNumbering;
	}

	public void addVertex( Vertex  a ) {
		if(edgesByVertexA.containsKey(a))
			return;
		edgesByVertexA.put( a, new HashSet<Edge>() );
		edgesByVertexB.put( a, new HashSet<Edge>() );
	}

	public void cutVertex( Vertex a ) {
		edgesByVertexA.remove( a );
		edgesByVertexB.remove( a );
	}

	public void addEdge( Edge e ) {
		edgesByVertexA.get(e.a).add(e);
		edgesByVertexB.get(e.b).add(e);
	}

	public void cutEdge( Edge e ) {
		edgesByVertexA.get(e.a).remove(e);
		edgesByVertexB.get(e.b).remove(e);
	}

	public float degree ( Vertex v ) {
		int index = numbering.get(v);
		return matrix[Graph.MatrixTypeIndex.DEGREE].entry[index][index];  
	}

	public void crystallizeGraph() {
		numbering = null;
		edgeNumbering = null;
		numbering();
		edgeNumbering();
		Set<Vertex> vertices = numbering.keySet();
		v = new Vertex[vertices.size()];
		vertices.toArray(v);
		int j = v.length;
		while(--j>=0) {
			Graph.Vertex vj = v[j];
			numbering.put(vj, new Integer(j));
		}
		Set<Edge> edges = edgeNumbering.keySet();
		e = new Edge[edges.size()];
		edges.toArray(e);
		int i = Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT;
		j = 0;
		while(j < i) {
			Algebraicist.form(this, j++);  
		}
	}

	public void draw() {
		int i = e.length;
		if(aspektWorld.layoutOn == 1) {
			if(aspektWorld.edgesOn<32) {
				while(--i >= 0) {
				    if((aspektWorld.edgesOn&1)==1 && e[i].a.color == red) {
				      e[i].draw();
				    }
				    else if((aspektWorld.edgesOn&2)==2 && e[i].a.color == yellow) {
				      e[i].draw();
				    }
				    else if((aspektWorld.edgesOn&4)==4 && e[i].a.color == green) {
				      e[i].draw();
				    }
				    else if((aspektWorld.edgesOn&8)==8 && e[i].a.color == blue) {
				      e[i].draw();
				    }
				    else if((aspektWorld.edgesOn&16)==16 && e[i].a.color == gray) {
				      e[i].draw();
				    }
				}
			}
			else {
				while(--i >= 0) {
				    if((aspektWorld.edgesOn&1)==1 && e[i].a.color == red && e[i].b.color == red) {
				      e[i].draw();
				    }
				    else if((aspektWorld.edgesOn&2)==2 && e[i].a.color == yellow && e[i].b.color == yellow) {
				      e[i].draw();
				    }
				    else if((aspektWorld.edgesOn&4)==4 && e[i].a.color == green && e[i].b.color == green) {
				      e[i].draw();
				    }
				    else if((aspektWorld.edgesOn&8)==8 && e[i].a.color == blue && e[i].b.color == blue) {
				      e[i].draw();
				    }
				    else if((aspektWorld.edgesOn&16)==16 && e[i].a.color == gray && e[i].b.color == gray) {
				      e[i].draw();
				    }
				}
			}
		}
		else {
                    while(--i >= 0)
			e[i].draw();
		}
		i = v.length;
		if(aspektWorld.verticesOn==1)
			while(--i >= 0)
				v[i].draw();
	}

}

import processing.core.PVector;


class Line {
	/**
	 * 
	 */
	private final AspektWorld aspektWorld;
	public PVector another, posa;
	public PVector yetanother;
	public int lc;

	public Line( AspektWorld aspektWorld, float ax, float ay, float az, float tx, float ty, float tz ) {
		this.aspektWorld = aspektWorld;
		another = new PVector(ax,ay,az);
		yetanother = new PVector(tx,ty,tz);
		lc = this.aspektWorld.color (180,200,200,200);
	}  

	public void draw() {
		this.aspektWorld.stroke(lc);
		this.aspektWorld.strokeWeight(aspektWorld.MIDDLE_STROKE);
		this.aspektWorld.line(posa.x+another.x,posa.y+another.y,posa.z+another.z,posa.x+yetanother.x,posa.y+yetanother.y,posa.z+yetanother.z);  
	}
}

import java.util.Arrays;
import java.util.HashSet;

public class Matrix {
	public int rows, columns, shortest, symmetric;
	public int i, j;
	public static final int NONE = -1;
	public static final int ZERO = 0;
	public static final int IDENTITY = 1;
	public static final int COUNTERIDENTITY = 2;
	public static final int ONES = 3;
        public static final int COLLOCATION = 5;
	public static final int entry_changed = 1;
	public static final int det_changed = 2;
	public static final int eigen_changed = 4;
	public static final int norms_changed = 8;
	public static final int lu_changed = 16;
	public static final int inverse_changed = 32;
	public static final int minmax_changed = 64;
	public static final int all_changed = 127;
	public static final float ZEROF = 0.0f;
	public static final float IDENTITYF = 1.0f;
	public float eij;
	public float entry[][];
        public HashSet<Integer> hentry[][];
	public float trace;
	public float det;
	public float[] minmax;
	public float eigenvector[][];
	public float eigenvalue[];
	public float krylovs[][];
	public float norms[];
	public float lu_entry[][];
	public float inverse_entry[][];
	public int updates;
	public int eigens_calculated;
	public int norms_calculated;

	public Matrix(int TYPE, int rows_num, int columns_num) {
		rows = rows_num;
		columns = columns_num;
		symmetric = rows == columns ? 1 : 0;
		shortest = rows < columns ? rows : columns;
		updates = all_changed;
		switch (TYPE) {
		case (NONE):
		case (ZERO):
			make_zero();
		break;
		case (IDENTITY):
			make_identity();
		break;
		case (COUNTERIDENTITY):
			make_counterIdentity();
		break;
		case (ONES):
			make_ones();
                case ( COLLOCATION ):
                        make_collocation();
		break;
		}
		eigenvector = new float[rows][columns];
		eigenvalue = new float[rows];
	}

	public void show() {
		int i = rows;
		while(--i >= 0)
			show_v(entry[i]);
		System.out.println(rows + "," + columns);
	}

	public static void show_v (float[] v) {
		int i = v.length;
		while(--i >= 0)
			System.out.print(v[i] + " ");
		System.out.println();
	}

	public float[] minmax() {
		if (minmax!=null && (updates & minmax_changed) == 0)
			return minmax;
		minmax = new float[2];
		if (rows == 0 || columns == 0)
			return minmax;
		minmax[0] = minmax[1] = entry[0][0];
		i = rows;
		while (--i >= 0) {
			j = columns;
			while (--j >= 0) {
				eij = entry[i][j];
				if (eij > minmax[1])
					minmax[1] = eij;
				else if (eij < minmax[0])
					minmax[0] = eij;
			}
		}
		updates &= (all_changed ^ minmax_changed);
		return minmax;
	}

	public float[] minmaxRow(int row) {
		float[] minmaxRow = new float[2];
		if (rows < row || columns == 0)
			return minmaxRow;
		minmaxRow[0] = minmaxRow[1] = entry[row][0];
		j = columns;
		while (--j >= 0) {
			eij = entry[row][j];
			if (eij > minmaxRow[1])
				minmaxRow[1] = eij;
			else if (eij < minmaxRow[0])
				minmaxRow[0] = eij;
		}
		return minmaxRow;
	}

	public float trace() {
		i = shortest;
		trace = 0.0f;
		while (--i >= 0)
			trace += entry[i][i];
		return trace;
	}

        public float traces() {
        	i = shortest;
		trace = 0.0f;
		while (--i >= 0)
			trace += entry[i][i]*entry[i][i];
		return trace;
	  
        }

	public void make_zero() {
		entry = new float[rows][columns];
		updates &= (all_changed ^ entry_changed);
	}

        public void make_collocation() {
	        hentry = new HashSet[rows][columns];
		updates &= (all_changed ^ entry_changed);
	}

	public void make_identity() {
		entry = new float[rows][columns];
		i = shortest;
		while (--i >= 0)
			entry[i][i] = IDENTITYF;
		updates &= (all_changed ^ entry_changed);
	}

	public void make_counterIdentity() {
		entry = new float[rows][columns];
		i = shortest;
		while (--i >= 0)
			entry[i][shortest - i - 1] = IDENTITYF;
		updates &= (all_changed ^ entry_changed);
	}

	public void make_ones() {
		entry = new float[rows][columns];
		i = rows;
		while (--i >= 0) {
			Arrays.fill(entry[i], IDENTITYF);
		}
	}

	public static void copy(float[][] source, float[][] dest) {
		int si = dest.length;
		int sj = dest[0].length;
		while (--si >= 0) {
			System.arraycopy(source[si], 0, dest[si], 0, sj);
		}
	}

	public Matrix(float[][] entries) {
		rows = entries.length;
		if (rows != 0)
			columns = entries[0].length;
		else
			columns = 0;
		symmetric = rows == columns ? 1 : 0;
		shortest = rows < columns ? rows : columns;
		entry = new float[rows][columns];
		Matrix.copy(entries,entry);
		eigenvector = new float[rows][columns];
		eigenvalue = new float[rows];
		updates = all_changed;
	}

	public Matrix(Matrix expired) {
		rows = expired.rows;
		columns = expired.columns;
		entry = new float[rows][columns];
		eigenvector = new float[rows][columns];
		eigenvalue = new float[rows];
		Matrix.copy(expired.entry, entry);
		Matrix.copy(expired.eigenvector, eigenvector);
		System.arraycopy(expired.eigenvalue, 0, eigenvalue, 0,
				eigenvalue.length);
		shortest = expired.shortest;
	}

	public Matrix transpose() {
		int j = columns;
		int i = rows;
		Matrix transpose = new Matrix(Matrix.NONE, j, i);
		while (--j >= 0) {
			while (--i >= 0) {
				transpose.entry[j][i] = entry[i][j];
			}
			i = rows;
		}
		return transpose;
	}

}

import java.util.ArrayList;
import java.util.HashMap;
import java.util.HashSet;
import java.util.Iterator;

import processing.core.PApplet;
import processing.core.PConstants;
import processing.core.PVector;

public class Visionary {
	public final static int VERTICESONSCREEN = 333;
	public final static int NOMOD = 0;
	public final static int MOD = 1;
	public int[] gridsz;
	public int mc = 0;
	public Matrix lrf;
	public Matrix gridv;
	public int etype = 0;
	public int offset = 0;
	public float conMetParam = 1.0f;
	public float conMetStep = 0.5f;
	public float[] conMetMinMax;

	public final AspektWorld aspektWorld;
	public Graph g1;

	public Visionary(AspektWorld aw) {
		aspektWorld = aw;
		gridsz = new int[] { 70,140,210,280,350,420,490,560,630,700 };
		lrf = new Matrix(Matrix.IDENTITY,VERTICESONSCREEN, VERTICESONSCREEN);
	}

	public Graph randomGraph() {
		Graph g = new Graph(aspektWorld);
		ArrayList<Graph.Vertex> knownVertices = new ArrayList<Graph.Vertex>();
		int i = VERTICESONSCREEN;
		while(--i >= 0) {
			Graph.Vertex v = g.new Vertex();
			g.addVertex(v);
			knownVertices.add(v);
			int k,j;
			k = VERTICESONSCREEN-i;
			j = (int)Math.log(aspektWorld.random(k));
			while(--j >= 0) {
				Graph.Vertex randomNeighbour = knownVertices.get((int)(aspektWorld.random(k)));

				Graph.Edge newEdge = g.new Edge(v, randomNeighbour);
				if(j == 3) {
					Graph.Edge newEdge2 = g.new Edge(randomNeighbour,v);
					g.addEdge(newEdge2);
				}
				g.addEdge(newEdge);

			}
		}
		g.crystallizeGraph();
		g1 = g;
		return g;
	}


	public Graph randomizeLayoutColorSize( ) {
		int i = g1.v.length;
		while(--i >= 0) {
			Graph.Vertex v = g1.v[i];
			v.color = aspektWorld.color(aspektWorld.random(360),200,200,200);
			v.radius = 53;
			float x = aspektWorld.random(4*aspektWorld.width)-2*aspektWorld.width;
			float z= aspektWorld.random(4*aspektWorld.width)-2*aspektWorld.width;
			float y = aspektWorld.random(4*aspektWorld.width)-2*aspektWorld.width;
			v.pos.set(x,-y-300,z);
		}
		i = g1.e.length;
		while(--i >= 0) {
			Graph.Edge e = g1.e[i];
		}
		return g1;
	}

	public void squaresTrianglesCirclesLayout() {
		int i = g1.v.length;
		while(--i >= 0) {
			g1.v[i].pos.set(0.0f,0.0f,0.0f);
		}
		ArrayList tcs = findCycles(4);
		ArrayList fcs = findCycles(3);
		ArrayList wcs = findCycles(2);
		ArrayList ocs = findCycles(1);
		embedRegular(tcs, 4);
		embedRegular(fcs, 3);
		embedRegular(wcs, 2);
		embedRegular(ocs, 1);
		embedRest();
	}

	public void embedRest() {
		int j = g1.v.length;
		int xcount = 0;
		int ycount = 0;
		int zcount = 0;
		while(--j >= 0) {
			Graph.Vertex v = g1.v[j];
			if(v.pos.y == 0.0f) {
				v.pos.set(xcount,ycount,zcount);
				xcount-=150;
				v.color = aspektWorld.color(0,0,200,100);
			}
			if(xcount < -1000) {
				zcount -= 150;
				xcount =0;
			}
			if(zcount < -1000) {
				ycount -= 150;
				zcount = 0;

			}


		}

	}

	public ArrayList findCycles(int cnum) {
		Matrix a = g1.matrix[Graph.MatrixTypeIndex.ADJACENCY];
		Matrix h = g1.matrix[Graph.MatrixTypeIndex.COLLOCATION];
		int rows = a.rows;
		ArrayList cs = new ArrayList(a.rows*a.rows/2);
		int columns = a.columns;
		while(--rows >= 0) {
			while(--columns >= 0) {
				if(rows==columns && cnum != 1)
					continue;
				float aij = a.entry[rows][columns];
				if(aij == 1.0f) {
					switch(cnum) {
					case(1):
						cs.add(new int[] { rows} );
					break;
					case(2):
						if(a.entry[columns][rows] == 1.0f) {
							cs.add(new int[] { rows, columns } );
						}
					break;
					case(3):
						HashSet threes = h.hentry[rows][columns];
					int tsz = threes.size();
					Iterator it = threes.iterator();
					while(--tsz >= 0) {
						int third = (Integer)it.next();
						if(third != rows && third != columns) {
							int[] triangle = new int[] { rows, columns, third };
							cs.add(triangle);  
						}
					}
					break;
					case(4):
						int j = a.columns;
					while(--j >= 0) {
						if(j == columns || j == rows)
							continue;
						aij = a.entry[columns][j];
						if(aij == 1.0f) {
							HashSet fours = h.hentry[rows][j];
							int fsz = fours.size();
							it = fours.iterator();
							while(--fsz >= 0) {
								int fourth = (Integer)it.next();
								if(fourth != rows && fourth != columns && fourth != j) {
									int[] square = new int[] { rows, columns, j, fourth };
									cs.add(square);  
								}
							}
						}  
					}
					break;  
					}
				}       
			}
			columns = a.columns;
		}
		return cs;
	}

	public void embedRegular(ArrayList cycled, int size) {
		int j = cycled.size();
		int sz = 1;
		while(--j >= 0)
			sz = embedRegular((int[])cycled.get(j), size, sz);
	}

	public int embedRegular(int[] cycle, int size, int count) {
		int j = size;
		float ang = 2*PConstants.PI/size;
		float twist = PConstants.PI/size;
		float radius = size*200;
		int vw = 200;
		int vh = 168;
		PVector centroid = new PVector(vw*(10),-count*vh,vw*(10*size)-vw*20);
		int vc = g1.v[cycle[0]].color;
		while(--j >= 0) {
			Graph.Vertex v = g1.v[cycle[j]];
			if(v.pos.y != 0.0f)
				return count;
		}
		j = size;
		while(--j >= 0) {
			Graph.Vertex v = g1.v[cycle[j]];
			switch(size) {
			case(1):
				v.color = aspektWorld.color(0,200,200,200);
			break;
			case(2):
				v.color = aspektWorld.color(60,200,200,200);
			break;
			case(3):
				v.color = aspektWorld.color(119,255,250,200);
			break;
			case(4):
				v.color = aspektWorld.color(200,200,200,200);
			break;
			
			}
			v.pos.set(centroid.x+(radius*PApplet.sin(twist+j*ang)),centroid.y, centroid.z+(radius*PApplet.cos(twist+j*ang)));
			//v.pos.set(centroid);
		}
		return count+1;
	}

}