Given a polynomial p(x) = sum{i=0..n} a[i]*x^i, compute its companion matrix, i.e., a matrix whose characteristic polynomial is p. One such matrix is the one that is zero except for the entries m[i+1,i]:=1 and m[i,n-1]:=-a[i]/a[n]. This way, computing the roots of the polynomial can be done by computing the eigenvalues of its matrix.

This function computes a QR-decomposition of a given matrix by means of Givens rotations, i.e., for a given matrix m, it computes an orthogonal matrix q (the column vectors of q are orthogonal unit length vectors) and an upper triangular matrix r such that m=q*r holds. Despite the real size of the matrix m, it will only consider the upper left submatrix of size dimRow x dimCol. Note that the inverse of an orthogonal matrix is its transpose.

This function computes a QR-decomposition of a given matrix by means of Givens rotations, i.e., for a given matrix m, it computes an orthogonal matrix q (the column vectors of q are orthogonal unit length vectors) and an upper triangular matrix r such that m=q*r holds. Despite the real size of the matrix m, it will only consider the upper left submatrix of size dimRow x dimCol. Note that the inverse of an orthogonal matrix is its transpose. See the comments on DecomposeQR to understand the implementation.

This function computes an iteration for approximating the eigenvalues of a given matrix m over real numbers. To this end, we define the QR iteration with double shifts as follows:

• m[0] := m
• m[i+1] := trans(q[i])*m[i]*q[i] for the QR-decomposition q[i]*r[i] := m[i]*m[i] - (s1+s2)*m[i] + s1*s2*I using shifts s1,s2

Note that this is indeed a double shift with the sift values s1 and s2: To see this, define the following values as would be done with two linear shifts with values s1 and s2:

• b := r1*q1+s1*I with QR decomposition q1*r1=m[i]-s1*I
• m[i+1] := r2*q2+s2*I with QR decomposition q2*r2=b-s2*I

It will compute the same m[i+1] as our definition above:

• We have m[i+1] := trans(q[i])*m[i]*q[i] with q[i]:=q1*q2: m[i+1] = r2*q2+s2*I = (trans(q2)*(b-s2*I))*q2+s2*I = trans(q2)*b*q2 - s2*trans(q2)*q2+s2*I = trans(q2)*b*q2 = trans(q2)*(r1*q1+s1*I)*q2 = trans(q2)*(trans(q1)*(m[i]-s1*I)*q1+s1*I)*q2 = trans(q2)*trans(q1)*m[i]*q1*q2 - s1*trans(q2)*(trans(q1)*q1*q2 + s1*trans(q2)*q2 = trans(q2)*trans(q1)*m[i]*q1*q2 = trans(q)*m[i]*q
• Second, we have q1*q2*r2*r1 = (m[i]-s2*I)*(m[i]-s1*I): q1*q2*r2*r1 = q1*(q2*r2)*r1 = q1*(b-s2*I)*r1 = (m[i]-s2*I)*q1*r1 = (m[i]-s2*I)*(m[i]-s1*I), where q1*(b-s2*I) = (m[i]-s2*I)*q1 can be seen as follows: q1*(b-s2*I) = (m[i]-s2*I)*q1 iff q1*b-s2*q1 = m[i]*q1-s2*q1 iff q1*b = m[i]*q1 iff q1*(r1*q1+s1*I) = m[i]*q1 iff (q1*r1)*q1+s1*q1 = m[i]*q1 iff (m[i]-s1*I)*q1+s1*q1 = m[i]*q1 iff m[i]*q1-s1*q1+s1*q1 = m[i]*q1 iff m[i]*q1 = m[i]*q1
• Third, note that (m[i]-s2*I)*(m[i]-s1*I) is nothing else but the polynomial m[i]*m[i] - (s1+s2)*m[i] + s1*s2*I used in the first iteration. For complex conjugate numbers s1,s2, both s1+s2 and s1*s2 are real numbers, so that only real number computations are required.
The above steps define therefore a double step QR iteration where for a real valued matrix only computations with real numbers are required even though some of the eigenvalues may be complex numbers. For the shift values, one considers the eigenvalues of the lower right 2x2 submatrix.

This function computes an iteration for approximating the eigenvalues of a given matrix m over complex numbers. To this end, we define the QR iteration with linear shifts as follows:

• m[0] := m
• m[i+1] := r[i]*q[i]+s[i]*I for the QR-decomposition q[i]*r[i] := m[i]-s[i]*I using some shift value s[i] (a complex number).

Since trans(q[i]) is the inverse of q[i], the QR decomposition implies that r[i]=trans(q[i])*(m[i]-s[i]*I) holds, and thus m[i+1] = r[i]*q[i]+s[i]*I = trans(q[i])*(m[i]-s[i]*I)*q[i]+s[i]*I = trans(q[i])*m[i]*q[i]-s[i]*trans(q[i])*q[i]+s[i]*I = trans(q[i])*m[i]*q[i]. Since m[i+1] = trans(q[i])*m[i]*q[i] holds, it follows that m[i] and m[i+1] have the same eigenvalues: Note that m[i+1]*x = lambda*x implies trans(q[i])*m[i]*q[i]*x = lambda*x, thus, m[i]*q[i]*x = lambda*(q[i]*x), and m[i]*x = lambda*x implies trans(q[i])*m[i]*x = lambda*trans(q[i])*x, which in turn implies trans(q[i])*m[i]*q[i]*x' = lambda*x' for x' := trans(q[i])*x. Hence, if the sequence of matrices m[i] will converge, it will converge to a triangular matrix, where we can read its eigenvalues from its diagonal. However, convergence is poor unless we will make use of shifts. Good values for a shift are approximations of eigenvalues, in particular, we make use of the Wilkinson shift.

This function implements the Fourier-Motzkin variable elimination to determine a solution of a linear inequality system.

GaussJordan transforms the given matrix into a diagonal matrix of using in-place operations that add multiples of another row to a row. It is the basis for solving linear equation systems or for inverting matrices.

InverseMatrixRat inverts a given matrix in that it concatenates the unit matrix and applies then the Gauss-Jordan algorithm. The given matrix can be inverted if the left half is finally the unit matrix, and in that case the right half contains the inverse of the given matrix.

Convert a linear inequality system to a GeneralSimplexProblem

convert a matrix of floats to a matrix of complex numbers

This function is to be called by the html page. It expects in the qscoll for the key "MatrixTask" one of the following DecomposeQR, EigenvaluesShift EigenvaluesDoubleShift, GaussJordan, Inverse, LinearEquationSystem and will then call one of the functions of this module with the appropriate html output.

Parse a GeneralSimplexProblem from a string as the following let s = "-oo <= 5 1 <= 5 -oo <= 2 1 <= 4 -oo <= 1/2 1 <= 2 -oo <= 1/5 1 <= 1 1 <= 1/2 1 <= +oo 3/2 <= 1 1 <= +oo"

Parse a string to a linear inequality system.

Parse a matrix from a string where rows are separated by newline symbols or semicolons, and entries within a row are separated by blanks or tabs. The entries of the string are expected as floating point numbers.

Parse a matrix from a string where rows are separated by newline symbols or semicolons, and entries within a row are separated by blanks or tabs. The entries of the string are expected as rational numbers.

This function is to be called by the html page. It will parse a polynomial and will compute its roots via the eigenvalues of its companion matrix

print a GeneralSimplex in ASCII

print a GeneralSimplex in Latex

print a GeneralSimplex in MathML

Print a linear inequality system to output stream ostr. The additional array vars holds the names of the variables in the linear inequality system.

Print a linear inequality system to output stream ostr in Latex. The additional array vars holds the names of the variables in the linear inequality system.

Print a linear inequality system to output stream ostr in MathML. The additional array vars holds the names of the variables in the linear inequality system.

print a matrix to output stream ostr

print a matrix to output stream ostr

print a matrix to output stream ostr

SimplexIntSolve solves the integer linear programming problem given by gs. Essentially, it first calls SimplexRatSolve to determine a rational solution if there is one, and narrows then the constraints to integer bounds until an integer solution is found. It therefore solves an integer linear programming problem which is known to be NP-complete. Note however, that it is possible that this method will not terminate since there are inequation systems with unbounded rational solutions that have no integer solutions. For this reason, one can specify an optional maximal number of recursion steps, if that is reached, the function returns None, otherwise, it will find a solution and

Simplex algorithm as decision procedure: Given a general Simplex problem (see the data type of this module, this function performs the required pivoting steps to determine a variable assignment in the contained eval component that satisfies the linear equation system of the tableau and also the constraints. It returns true if a solution has been found, and false if there is no solution on the rational numbers. The computations are done in-place in the given gs. Note that the Simplex algorithm may have an exponential runtime (which shows up rarely in practise), but it is known that the problem can be decided on the rational numbers in polynomial time (Kharmarkar's algorithm, and Kachian's ellipsoid method).

This function is to be called by a teaching tool's html page. It will parse a general simplex problem, will solve it and print the results to the html page.

Solves a linear equation system A*x=0 or A*x=b, where in the first case the given matrix is just A, and in the second case, the given matrix should be matrix (A|b). The result is a tuple (defVarSet,parVarSet,kern) where defVarSet is the set of variable indices that are defined in terms of the other ones contained in the set of parameters parVarSet, and kern is a matrix whose image is the kernel of the given matrix m. For the homogeneous equation system A*x=0, this means that all linear combinations of the column vectors of kern are solutions, and there is just the trivial solution x=0 iff the kernel is a dimCol x 0 matrix. For the inhomogeneous equation system A*x=b, there are no solutions if the variable associated with the rightmost column is a defined variable. Otherwise, the solutions of the inhomogeneous system are obtained from the solutions of the homogeneous system by setting the last parameter -1 and removing the last rows in the kernel's column vectors.

tensor product

solving a linear equation system with an upper triangular matrix uMat and right hand side vector y

solving a linear equation system with a lower triangular matrix lMat and right hand side vector b

Perform a LU-decomposition of the given matrix which is not modified by this function unlike in function luDecomposeInline.

Perform a LU-decomposition of the given matrix where the matrices L and U are stored in the given matrix (L has an additional diagonal of 1s).

Perform a LUP-decomposition of the given matrix which is not modified by this function unlike in function luDecomposeInline. The function returns the matrices L,U,P.

Perform a LUP-decomposition of the given matrix where the matrices L and U are stored in the given matrix (L has an additional diagonal of 1s) and the permutation matrix is returned as a permutation.

solving a linear equation system A*x=b by a LUP decomposition P*A=L*U where the LUP-decomposition is given

matrix * matrix product

matrix * vector product