From a3762a99e6731cc8c633967c85a24f2de6e728d1 Mon Sep 17 00:00:00 2001 From: Daniel Weschke Date: Sat, 15 Mar 2014 10:11:53 +0100 Subject: [PATCH] Add some matrix classes in math\matrix package and standard libraries --- src/math/matrix/CholeskyDecomposition.java | 201 ++++ src/math/matrix/Diagonal.java | 175 ++++ src/math/matrix/Diagonal2D.java | 47 + src/math/matrix/EigenvalueDecomposition.java | 961 ++++++++++++++++++ src/math/matrix/LUDecomposition.java | 282 +++++ src/math/matrix/Matrix2D.java | 334 ++++++ src/math/matrix/Matrix3D.java | 193 ++++ src/math/matrix/QRDecomposition.java | 197 ++++ .../matrix/SingularValueDecomposition.java | 467 +++++++++ src/math/matrix/TensorII.java | 13 + src/math/matrix/TensorII2D.java | 58 ++ src/math/matrix/Vector.java | 899 ++++++++++++++++ src/math/matrix/Vector2D.java | 80 ++ src/math/matrix/Vector3D.java | 120 +++ src/stdlib/StdArrayIO.java | 256 +++++ src/stdlib/StdIn.java | 210 ++++ src/stdlib/StdOut.java | 230 +++++ 17 files changed, 4723 insertions(+) create mode 100644 src/math/matrix/CholeskyDecomposition.java create mode 100644 src/math/matrix/Diagonal.java create mode 100644 src/math/matrix/Diagonal2D.java create mode 100644 src/math/matrix/EigenvalueDecomposition.java create mode 100644 src/math/matrix/LUDecomposition.java create mode 100644 src/math/matrix/Matrix2D.java create mode 100644 src/math/matrix/Matrix3D.java create mode 100644 src/math/matrix/QRDecomposition.java create mode 100644 src/math/matrix/SingularValueDecomposition.java create mode 100644 src/math/matrix/TensorII.java create mode 100644 src/math/matrix/TensorII2D.java create mode 100644 src/math/matrix/Vector.java create mode 100644 src/math/matrix/Vector2D.java create mode 100644 src/math/matrix/Vector3D.java create mode 100644 src/stdlib/StdArrayIO.java create mode 100644 src/stdlib/StdIn.java create mode 100644 src/stdlib/StdOut.java diff --git a/src/math/matrix/CholeskyDecomposition.java b/src/math/matrix/CholeskyDecomposition.java new file mode 100644 index 0000000..367e911 --- /dev/null +++ b/src/math/matrix/CholeskyDecomposition.java @@ -0,0 +1,201 @@ +package math.matrix; + +import java.io.Serializable; + + /** Cholesky Decomposition. +

+ For a symmetric, positive definite matrix A, the Cholesky decomposition + is an lower triangular matrix L so that A = L*L'. +

+ If the matrix is not symmetric or positive definite, the constructor + returns a partial decomposition and sets an internal flag that may + be queried by the isSPD() method. + */ + +public class CholeskyDecomposition implements Serializable { + +/* ------------------------ + Class variables + * ------------------------ */ + + /** Array for internal storage of decomposition. + @serial internal array storage. + */ + private double[][] L; + + /** Row and column dimension (square matrix). + @serial matrix dimension. + */ + private int n; + + /** Symmetric and positive definite flag. + @serial is symmetric and positive definite flag. + */ + private boolean isspd; + +/* ------------------------ + Constructor + * ------------------------ */ + + /** Cholesky algorithm for symmetric and positive definite matrix. + Structure to access L and isspd flag. + @param Arg Square, symmetric matrix. + */ + + public CholeskyDecomposition (Matrix Arg) { + // Initialize. + double[][] A = Arg.getD(); + n = Arg.getM(); + L = new double[n][n]; + isspd = (Arg.getN() == n); + // Main loop. + for (int j = 0; j < n; j++) { + double[] Lrowj = L[j]; + double d = 0.0; + for (int k = 0; k < j; k++) { + double[] Lrowk = L[k]; + double s = 0.0; + for (int i = 0; i < k; i++) { + s += Lrowk[i]*Lrowj[i]; + } + Lrowj[k] = s = (A[j][k] - s)/L[k][k]; + d = d + s*s; + isspd = isspd & (A[k][j] == A[j][k]); + } + d = A[j][j] - d; + isspd = isspd & (d > 0.0); + L[j][j] = Math.sqrt(Math.max(d,0.0)); + for (int k = j+1; k < n; k++) { + L[j][k] = 0.0; + } + } + } + +/* ------------------------ + Temporary, experimental code. + * ------------------------ *\ + + \** Right Triangular Cholesky Decomposition. +

+ For a symmetric, positive definite matrix A, the Right Cholesky + decomposition is an upper triangular matrix R so that A = R'*R. + This constructor computes R with the Fortran inspired column oriented + algorithm used in LINPACK and MATLAB. In Java, we suspect a row oriented, + lower triangular decomposition is faster. We have temporarily included + this constructor here until timing experiments confirm this suspicion. + *\ + + \** Array for internal storage of right triangular decomposition. **\ + private transient double[][] R; + + \** Cholesky algorithm for symmetric and positive definite matrix. + @param A Square, symmetric matrix. + @param rightflag Actual value ignored. + @return Structure to access R and isspd flag. + *\ + + public CholeskyDecomposition (Matrix Arg, int rightflag) { + // Initialize. + double[][] A = Arg.getArray(); + n = Arg.getColumnDimension(); + R = new double[n][n]; + isspd = (Arg.getColumnDimension() == n); + // Main loop. + for (int j = 0; j < n; j++) { + double d = 0.0; + for (int k = 0; k < j; k++) { + double s = A[k][j]; + for (int i = 0; i < k; i++) { + s = s - R[i][k]*R[i][j]; + } + R[k][j] = s = s/R[k][k]; + d = d + s*s; + isspd = isspd & (A[k][j] == A[j][k]); + } + d = A[j][j] - d; + isspd = isspd & (d > 0.0); + R[j][j] = Math.sqrt(Math.max(d,0.0)); + for (int k = j+1; k < n; k++) { + R[k][j] = 0.0; + } + } + } + + \** Return upper triangular factor. + @return R + *\ + + public Matrix getR () { + return new Matrix(R,n,n); + } + +\* ------------------------ + End of temporary code. + * ------------------------ */ + +/* ------------------------ + Public Methods + * ------------------------ */ + + /** Is the matrix symmetric and positive definite? + @return true if A is symmetric and positive definite. + */ + + public boolean isSPD () { + return isspd; + } + + /** Return triangular factor. + @return L + */ + + public Matrix getL () { + return new Matrix(L); + } + + /** Solve A*X = B + @param B A Matrix with as many rows as A and any number of columns. + @return X so that L*L'*X = B + @exception IllegalArgumentException Matrix row dimensions must agree. + @exception RuntimeException Matrix is not symmetric positive definite. + */ + + public Matrix solve (Matrix B) { + if (B.getM() != n) { + throw new IllegalArgumentException("Matrix row dimensions must agree."); + } + if (!isspd) { + throw new RuntimeException("Matrix is not symmetric positive definite."); + } + + // Copy right hand side. + double[][] X = B.getCopy(); + int nx = B.getN(); + + // Solve L*Y = B; + for (int k = 0; k < n; k++) { + for (int j = 0; j < nx; j++) { + for (int i = 0; i < k ; i++) { + X[k][j] -= X[i][j]*L[k][i]; + } + X[k][j] /= L[k][k]; + } + } + + // Solve L'*X = Y; + for (int k = n-1; k >= 0; k--) { + for (int j = 0; j < nx; j++) { + for (int i = k+1; i < n ; i++) { + X[k][j] -= X[i][j]*L[i][k]; + } + X[k][j] /= L[k][k]; + } + } + + + return new Matrix(X); + } + private static final long serialVersionUID = 1; + +} + diff --git a/src/math/matrix/Diagonal.java b/src/math/matrix/Diagonal.java new file mode 100644 index 0000000..68fc447 --- /dev/null +++ b/src/math/matrix/Diagonal.java @@ -0,0 +1,175 @@ +package math.matrix; + +import exception.IllegalDimensionException; + +public class Diagonal extends Matrix{ + + /** + * UID + */ + private static final long serialVersionUID = 3206696282474150081L; + + public Diagonal(int n){ + this(n,n); + } + + public Diagonal(int m, int n){ + this.m = m; + this.n = n; + data = new double[m>n?m:n]; + } + + public Diagonal(double... d){ + this(d.length); + int i; + for(i=0; in?m:n); i++) + set(i, d.get(i, i)); + } + + /** + * Generate an n-by-n identity matrix. + * An n-by-n matrix with ones on the diagonal and zeros elsewhere. + * @param n rows/columns + * @return n-by-n identity matrix + */ + public static Diagonal identity(int n){ + int i; + Diagonal c = new Diagonal(n); + for(i=0; iC = s A + */ + public Diagonal times(double s){ + int i; + Diagonal c = new Diagonal(m,n); + for(i=0; iD = s D + */ + public Diagonal timesEquals(double s){ + int i; + for(i=0; iC = A • b = A b + * @throws IllegalDimensionException Illegal vector dimensions. + */ + public Vector times(Vector b) throws IllegalDimensionException{ + int i; + if(n != b.n()) throw new IllegalDimensionException("Illegal vector dimensions."); + Vector c = new Vector(m); + for(i=0; in?m:n]; + } + + /** + * Generate an n-by-n identity matrix. + * An m-by-n matrix with ones on the diagonal and zeros elsewhere. + * @param n rows/columns + * @return n-by-n identity matrix + */ + public static Diagonal2D identity(int n){ + int i; + Diagonal2D c = new Diagonal2D(n); + for(i=0; i + If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is + diagonal and the eigenvector matrix V is orthogonal. + I.e. A = V.times(D.times(V.transpose())) and + V.times(V.transpose()) equals the identity matrix. +

+ If A is not symmetric, then the eigenvalue matrix D is block diagonal + with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, + lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The + columns of V represent the eigenvectors in the sense that A*V = V*D, + i.e. A.times(V) equals V.times(D). The matrix V may be badly + conditioned, or even singular, so the validity of the equation + A = V*D*inverse(V) depends upon V.cond(). +**/ + +public class EigenvalueDecomposition implements Serializable { + + /** + * UID + */ + private static final long serialVersionUID = -4489049767346210616L; + +/* ------------------------ + Class variables + * ------------------------ */ + + /** Row and column dimension (square matrix). + @serial matrix dimension. + */ + private int n; + + /** Symmetry flag. + @serial internal symmetry flag. + */ + private boolean issymmetric; + + /** Arrays for internal storage of eigenvalues. + @serial internal storage of eigenvalues. + */ + private double[] d, e; + + /** Array for internal storage of eigenvectors. + @serial internal storage of eigenvectors. + */ + private double[][] V; + + /** Array for internal storage of nonsymmetric Hessenberg form. + @serial internal storage of nonsymmetric Hessenberg form. + */ + private double[][] H; + + /** Working storage for nonsymmetric algorithm. + @serial working storage for nonsymmetric algorithm. + */ + private double[] ort; + +/* ------------------------ + Private Methods + * ------------------------ */ + + // Symmetric Householder reduction to tridiagonal form. + + private void tred2 () { + + // This is derived from the Algol procedures tred2 by + // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for + // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutine in EISPACK. + + for (int j = 0; j < n; j++) { + d[j] = V[n-1][j]; + } + + // Householder reduction to tridiagonal form. + + for (int i = n-1; i > 0; i--) { + + // Scale to avoid under/overflow. + + double scale = 0.0; + double h = 0.0; + for (int k = 0; k < i; k++) { + scale = scale + Math.abs(d[k]); + } + if (scale == 0.0) { + e[i] = d[i-1]; + for (int j = 0; j < i; j++) { + d[j] = V[i-1][j]; + V[i][j] = 0.0; + V[j][i] = 0.0; + } + } else { + + // Generate Householder vector. + + for (int k = 0; k < i; k++) { + d[k] /= scale; + h += d[k] * d[k]; + } + double f = d[i-1]; + double g = Math.sqrt(h); + if (f > 0) { + g = -g; + } + e[i] = scale * g; + h = h - f * g; + d[i-1] = f - g; + for (int j = 0; j < i; j++) { + e[j] = 0.0; + } + + // Apply similarity transformation to remaining columns. + + for (int j = 0; j < i; j++) { + f = d[j]; + V[j][i] = f; + g = e[j] + V[j][j] * f; + for (int k = j+1; k <= i-1; k++) { + g += V[k][j] * d[k]; + e[k] += V[k][j] * f; + } + e[j] = g; + } + f = 0.0; + for (int j = 0; j < i; j++) { + e[j] /= h; + f += e[j] * d[j]; + } + double hh = f / (h + h); + for (int j = 0; j < i; j++) { + e[j] -= hh * d[j]; + } + for (int j = 0; j < i; j++) { + f = d[j]; + g = e[j]; + for (int k = j; k <= i-1; k++) { + V[k][j] -= (f * e[k] + g * d[k]); + } + d[j] = V[i-1][j]; + V[i][j] = 0.0; + } + } + d[i] = h; + } + + // Accumulate transformations. + + for (int i = 0; i < n-1; i++) { + V[n-1][i] = V[i][i]; + V[i][i] = 1.0; + double h = d[i+1]; + if (h != 0.0) { + for (int k = 0; k <= i; k++) { + d[k] = V[k][i+1] / h; + } + for (int j = 0; j <= i; j++) { + double g = 0.0; + for (int k = 0; k <= i; k++) { + g += V[k][i+1] * V[k][j]; + } + for (int k = 0; k <= i; k++) { + V[k][j] -= g * d[k]; + } + } + } + for (int k = 0; k <= i; k++) { + V[k][i+1] = 0.0; + } + } + for (int j = 0; j < n; j++) { + d[j] = V[n-1][j]; + V[n-1][j] = 0.0; + } + V[n-1][n-1] = 1.0; + e[0] = 0.0; + } + + // Symmetric tridiagonal QL algorithm. + + private void tql2 () { + + // This is derived from the Algol procedures tql2, by + // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for + // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutine in EISPACK. + + for (int i = 1; i < n; i++) { + e[i-1] = e[i]; + } + e[n-1] = 0.0; + + double f = 0.0; + double tst1 = 0.0; + double eps = Math.pow(2.0,-52.0); + for (int l = 0; l < n; l++) { + + // Find small subdiagonal element + + tst1 = Math.max(tst1,Math.abs(d[l]) + Math.abs(e[l])); + int m = l; + while (m < n) { + if (Math.abs(e[m]) <= eps*tst1) { + break; + } + m++; + } + + // If m == l, d[l] is an eigenvalue, + // otherwise, iterate. + + if (m > l) { + int iter = 0; + do { + iter = iter + 1; // (Could check iteration count here.) + + // Compute implicit shift + + double g = d[l]; + double p = (d[l+1] - g) / (2.0 * e[l]); + double r = Maths.hypot(p,1.0); + if (p < 0) { + r = -r; + } + d[l] = e[l] / (p + r); + d[l+1] = e[l] * (p + r); + double dl1 = d[l+1]; + double h = g - d[l]; + for (int i = l+2; i < n; i++) { + d[i] -= h; + } + f = f + h; + + // Implicit QL transformation. + + p = d[m]; + double c = 1.0; + double c2 = c; + double c3 = c; + double el1 = e[l+1]; + double s = 0.0; + double s2 = 0.0; + for (int i = m-1; i >= l; i--) { + c3 = c2; + c2 = c; + s2 = s; + g = c * e[i]; + h = c * p; + r = Maths.hypot(p,e[i]); + e[i+1] = s * r; + s = e[i] / r; + c = p / r; + p = c * d[i] - s * g; + d[i+1] = h + s * (c * g + s * d[i]); + + // Accumulate transformation. + + for (int k = 0; k < n; k++) { + h = V[k][i+1]; + V[k][i+1] = s * V[k][i] + c * h; + V[k][i] = c * V[k][i] - s * h; + } + } + p = -s * s2 * c3 * el1 * e[l] / dl1; + e[l] = s * p; + d[l] = c * p; + + // Check for convergence. + + } while (Math.abs(e[l]) > eps*tst1); + } + d[l] = d[l] + f; + e[l] = 0.0; + } + + // Sort eigenvalues and corresponding vectors. + + for (int i = 0; i < n-1; i++) { + int k = i; + double p = d[i]; + for (int j = i+1; j < n; j++) { + if (d[j] < p) { + k = j; + p = d[j]; + } + } + if (k != i) { + d[k] = d[i]; + d[i] = p; + for (int j = 0; j < n; j++) { + p = V[j][i]; + V[j][i] = V[j][k]; + V[j][k] = p; + } + } + } + } + + // Nonsymmetric reduction to Hessenberg form. + + private void orthes () { + + // This is derived from the Algol procedures orthes and ortran, + // by Martin and Wilkinson, Handbook for Auto. Comp., + // Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutines in EISPACK. + + int low = 0; + int high = n-1; + + for (int m = low+1; m <= high-1; m++) { + + // Scale column. + + double scale = 0.0; + for (int i = m; i <= high; i++) { + scale = scale + Math.abs(H[i][m-1]); + } + if (scale != 0.0) { + + // Compute Householder transformation. + + double h = 0.0; + for (int i = high; i >= m; i--) { + ort[i] = H[i][m-1]/scale; + h += ort[i] * ort[i]; + } + double g = Math.sqrt(h); + if (ort[m] > 0) { + g = -g; + } + h = h - ort[m] * g; + ort[m] = ort[m] - g; + + // Apply Householder similarity transformation + // H = (I-u*u'/h)*H*(I-u*u')/h) + + for (int j = m; j < n; j++) { + double f = 0.0; + for (int i = high; i >= m; i--) { + f += ort[i]*H[i][j]; + } + f = f/h; + for (int i = m; i <= high; i++) { + H[i][j] -= f*ort[i]; + } + } + + for (int i = 0; i <= high; i++) { + double f = 0.0; + for (int j = high; j >= m; j--) { + f += ort[j]*H[i][j]; + } + f = f/h; + for (int j = m; j <= high; j++) { + H[i][j] -= f*ort[j]; + } + } + ort[m] = scale*ort[m]; + H[m][m-1] = scale*g; + } + } + + // Accumulate transformations (Algol's ortran). + + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + V[i][j] = (i == j ? 1.0 : 0.0); + } + } + + for (int m = high-1; m >= low+1; m--) { + if (H[m][m-1] != 0.0) { + for (int i = m+1; i <= high; i++) { + ort[i] = H[i][m-1]; + } + for (int j = m; j <= high; j++) { + double g = 0.0; + for (int i = m; i <= high; i++) { + g += ort[i] * V[i][j]; + } + // Double division avoids possible underflow + g = (g / ort[m]) / H[m][m-1]; + for (int i = m; i <= high; i++) { + V[i][j] += g * ort[i]; + } + } + } + } + } + + + // Complex scalar division. + + private transient double cdivr, cdivi; + private void cdiv(double xr, double xi, double yr, double yi) { + double r,d; + if (Math.abs(yr) > Math.abs(yi)) { + r = yi/yr; + d = yr + r*yi; + cdivr = (xr + r*xi)/d; + cdivi = (xi - r*xr)/d; + } else { + r = yr/yi; + d = yi + r*yr; + cdivr = (r*xr + xi)/d; + cdivi = (r*xi - xr)/d; + } + } + + + // Nonsymmetric reduction from Hessenberg to real Schur form. + + private void hqr2 () { + + // This is derived from the Algol procedure hqr2, + // by Martin and Wilkinson, Handbook for Auto. Comp., + // Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutine in EISPACK. + + // Initialize + + int nn = this.n; + int n = nn-1; + int low = 0; + int high = nn-1; + double eps = Math.pow(2.0,-52.0); + double exshift = 0.0; + double p=0,q=0,r=0,s=0,z=0,t,w,x,y; + + // Store roots isolated by balanc and compute matrix norm + + double norm = 0.0; + for (int i = 0; i < nn; i++) { + if (i < low | i > high) { + d[i] = H[i][i]; + e[i] = 0.0; + } + for (int j = Math.max(i-1,0); j < nn; j++) { + norm = norm + Math.abs(H[i][j]); + } + } + + // Outer loop over eigenvalue index + + int iter = 0; + while (n >= low) { + + // Look for single small sub-diagonal element + + int l = n; + while (l > low) { + s = Math.abs(H[l-1][l-1]) + Math.abs(H[l][l]); + if (s == 0.0) { + s = norm; + } + if (Math.abs(H[l][l-1]) < eps * s) { + break; + } + l--; + } + + // Check for convergence + // One root found + + if (l == n) { + H[n][n] = H[n][n] + exshift; + d[n] = H[n][n]; + e[n] = 0.0; + n--; + iter = 0; + + // Two roots found + + } else if (l == n-1) { + w = H[n][n-1] * H[n-1][n]; + p = (H[n-1][n-1] - H[n][n]) / 2.0; + q = p * p + w; + z = Math.sqrt(Math.abs(q)); + H[n][n] = H[n][n] + exshift; + H[n-1][n-1] = H[n-1][n-1] + exshift; + x = H[n][n]; + + // Real pair + + if (q >= 0) { + if (p >= 0) { + z = p + z; + } else { + z = p - z; + } + d[n-1] = x + z; + d[n] = d[n-1]; + if (z != 0.0) { + d[n] = x - w / z; + } + e[n-1] = 0.0; + e[n] = 0.0; + x = H[n][n-1]; + s = Math.abs(x) + Math.abs(z); + p = x / s; + q = z / s; + r = Math.sqrt(p * p+q * q); + p = p / r; + q = q / r; + + // Row modification + + for (int j = n-1; j < nn; j++) { + z = H[n-1][j]; + H[n-1][j] = q * z + p * H[n][j]; + H[n][j] = q * H[n][j] - p * z; + } + + // Column modification + + for (int i = 0; i <= n; i++) { + z = H[i][n-1]; + H[i][n-1] = q * z + p * H[i][n]; + H[i][n] = q * H[i][n] - p * z; + } + + // Accumulate transformations + + for (int i = low; i <= high; i++) { + z = V[i][n-1]; + V[i][n-1] = q * z + p * V[i][n]; + V[i][n] = q * V[i][n] - p * z; + } + + // Complex pair + + } else { + d[n-1] = x + p; + d[n] = x + p; + e[n-1] = z; + e[n] = -z; + } + n = n - 2; + iter = 0; + + // No convergence yet + + } else { + + // Form shift + + x = H[n][n]; + y = 0.0; + w = 0.0; + if (l < n) { + y = H[n-1][n-1]; + w = H[n][n-1] * H[n-1][n]; + } + + // Wilkinson's original ad hoc shift + + if (iter == 10) { + exshift += x; + for (int i = low; i <= n; i++) { + H[i][i] -= x; + } + s = Math.abs(H[n][n-1]) + Math.abs(H[n-1][n-2]); + x = y = 0.75 * s; + w = -0.4375 * s * s; + } + + // MATLAB's new ad hoc shift + + if (iter == 30) { + s = (y - x) / 2.0; + s = s * s + w; + if (s > 0) { + s = Math.sqrt(s); + if (y < x) { + s = -s; + } + s = x - w / ((y - x) / 2.0 + s); + for (int i = low; i <= n; i++) { + H[i][i] -= s; + } + exshift += s; + x = y = w = 0.964; + } + } + + iter = iter + 1; // (Could check iteration count here.) + + // Look for two consecutive small sub-diagonal elements + + int m = n-2; + while (m >= l) { + z = H[m][m]; + r = x - z; + s = y - z; + p = (r * s - w) / H[m+1][m] + H[m][m+1]; + q = H[m+1][m+1] - z - r - s; + r = H[m+2][m+1]; + s = Math.abs(p) + Math.abs(q) + Math.abs(r); + p = p / s; + q = q / s; + r = r / s; + if (m == l) { + break; + } + if (Math.abs(H[m][m-1]) * (Math.abs(q) + Math.abs(r)) < + eps * (Math.abs(p) * (Math.abs(H[m-1][m-1]) + Math.abs(z) + + Math.abs(H[m+1][m+1])))) { + break; + } + m--; + } + + for (int i = m+2; i <= n; i++) { + H[i][i-2] = 0.0; + if (i > m+2) { + H[i][i-3] = 0.0; + } + } + + // Double QR step involving rows l:n and columns m:n + + + for (int k = m; k <= n-1; k++) { + boolean notlast = (k != n-1); + if (k != m) { + p = H[k][k-1]; + q = H[k+1][k-1]; + r = (notlast ? H[k+2][k-1] : 0.0); + x = Math.abs(p) + Math.abs(q) + Math.abs(r); + if (x == 0.0) { + continue; + } + p = p / x; + q = q / x; + r = r / x; + } + + s = Math.sqrt(p * p + q * q + r * r); + if (p < 0) { + s = -s; + } + if (s != 0) { + if (k != m) { + H[k][k-1] = -s * x; + } else if (l != m) { + H[k][k-1] = -H[k][k-1]; + } + p = p + s; + x = p / s; + y = q / s; + z = r / s; + q = q / p; + r = r / p; + + // Row modification + + for (int j = k; j < nn; j++) { + p = H[k][j] + q * H[k+1][j]; + if (notlast) { + p = p + r * H[k+2][j]; + H[k+2][j] = H[k+2][j] - p * z; + } + H[k][j] = H[k][j] - p * x; + H[k+1][j] = H[k+1][j] - p * y; + } + + // Column modification + + for (int i = 0; i <= Math.min(n,k+3); i++) { + p = x * H[i][k] + y * H[i][k+1]; + if (notlast) { + p = p + z * H[i][k+2]; + H[i][k+2] = H[i][k+2] - p * r; + } + H[i][k] = H[i][k] - p; + H[i][k+1] = H[i][k+1] - p * q; + } + + // Accumulate transformations + + for (int i = low; i <= high; i++) { + p = x * V[i][k] + y * V[i][k+1]; + if (notlast) { + p = p + z * V[i][k+2]; + V[i][k+2] = V[i][k+2] - p * r; + } + V[i][k] = V[i][k] - p; + V[i][k+1] = V[i][k+1] - p * q; + } + } // (s != 0) + } // k loop + } // check convergence + } // while (n >= low) + + // Backsubstitute to find vectors of upper triangular form + + if (norm == 0.0) { + return; + } + + for (n = nn-1; n >= 0; n--) { + p = d[n]; + q = e[n]; + + // Real vector + + if (q == 0) { + int l = n; + H[n][n] = 1.0; + for (int i = n-1; i >= 0; i--) { + w = H[i][i] - p; + r = 0.0; + for (int j = l; j <= n; j++) { + r = r + H[i][j] * H[j][n]; + } + if (e[i] < 0.0) { + z = w; + s = r; + } else { + l = i; + if (e[i] == 0.0) { + if (w != 0.0) { + H[i][n] = -r / w; + } else { + H[i][n] = -r / (eps * norm); + } + + // Solve real equations + + } else { + x = H[i][i+1]; + y = H[i+1][i]; + q = (d[i] - p) * (d[i] - p) + e[i] * e[i]; + t = (x * s - z * r) / q; + H[i][n] = t; + if (Math.abs(x) > Math.abs(z)) { + H[i+1][n] = (-r - w * t) / x; + } else { + H[i+1][n] = (-s - y * t) / z; + } + } + + // Overflow control + + t = Math.abs(H[i][n]); + if ((eps * t) * t > 1) { + for (int j = i; j <= n; j++) { + H[j][n] = H[j][n] / t; + } + } + } + } + + // Complex vector + + } else if (q < 0) { + int l = n-1; + + // Last vector component imaginary so matrix is triangular + + if (Math.abs(H[n][n-1]) > Math.abs(H[n-1][n])) { + H[n-1][n-1] = q / H[n][n-1]; + H[n-1][n] = -(H[n][n] - p) / H[n][n-1]; + } else { + cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q); + H[n-1][n-1] = cdivr; + H[n-1][n] = cdivi; + } + H[n][n-1] = 0.0; + H[n][n] = 1.0; + for (int i = n-2; i >= 0; i--) { + double ra,sa,vr,vi; + ra = 0.0; + sa = 0.0; + for (int j = l; j <= n; j++) { + ra = ra + H[i][j] * H[j][n-1]; + sa = sa + H[i][j] * H[j][n]; + } + w = H[i][i] - p; + + if (e[i] < 0.0) { + z = w; + r = ra; + s = sa; + } else { + l = i; + if (e[i] == 0) { + cdiv(-ra,-sa,w,q); + H[i][n-1] = cdivr; + H[i][n] = cdivi; + } else { + + // Solve complex equations + + x = H[i][i+1]; + y = H[i+1][i]; + vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; + vi = (d[i] - p) * 2.0 * q; + if (vr == 0.0 & vi == 0.0) { + vr = eps * norm * (Math.abs(w) + Math.abs(q) + + Math.abs(x) + Math.abs(y) + Math.abs(z)); + } + cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi); + H[i][n-1] = cdivr; + H[i][n] = cdivi; + if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) { + H[i+1][n-1] = (-ra - w * H[i][n-1] + q * H[i][n]) / x; + H[i+1][n] = (-sa - w * H[i][n] - q * H[i][n-1]) / x; + } else { + cdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q); + H[i+1][n-1] = cdivr; + H[i+1][n] = cdivi; + } + } + + // Overflow control + + t = Math.max(Math.abs(H[i][n-1]),Math.abs(H[i][n])); + if ((eps * t) * t > 1) { + for (int j = i; j <= n; j++) { + H[j][n-1] = H[j][n-1] / t; + H[j][n] = H[j][n] / t; + } + } + } + } + } + } + + // Vectors of isolated roots + + for (int i = 0; i < nn; i++) { + if (i < low | i > high) { + for (int j = i; j < nn; j++) { + V[i][j] = H[i][j]; + } + } + } + + // Back transformation to get eigenvectors of original matrix + + for (int j = nn-1; j >= low; j--) { + for (int i = low; i <= high; i++) { + z = 0.0; + for (int k = low; k <= Math.min(j,high); k++) { + z = z + V[i][k] * H[k][j]; + } + V[i][j] = z; + } + } + } + + +/* ------------------------ + Constructor + * ------------------------ */ + + /** Check for symmetry, then construct the eigenvalue decomposition + Structure to access D and V. + @param Arg Square matrix + */ + + public EigenvalueDecomposition (Matrix Arg) { + double[][] A = Arg.getD(); + n = Arg.getN(); + V = new double[n][n]; + d = new double[n]; + e = new double[n]; + + issymmetric = true; + for (int j = 0; (j < n) & issymmetric; j++) { + for (int i = 0; (i < n) & issymmetric; i++) { + issymmetric = (A[i][j] == A[j][i]); + } + } + + if (issymmetric) { + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + V[i][j] = A[i][j]; + } + } + + // Tridiagonalize. + tred2(); + + // Diagonalize. + tql2(); + + } else { + H = new double[n][n]; + ort = new double[n]; + + for (int j = 0; j < n; j++) { + for (int i = 0; i < n; i++) { + H[i][j] = A[i][j]; + } + } + + // Reduce to Hessenberg form. + orthes(); + + // Reduce Hessenberg to real Schur form. + hqr2(); + } + } + +/* ------------------------ + Public Methods + * ------------------------ */ + + /** + * Return the eigenvector matrix + * @return V + */ + public Matrix getV () { + return new Matrix(V); + } + + /** + * Return the real parts of the eigenvalues + * @return real(diag(D)) + */ + public double[] getRealEigenvalues(){ + return d; + } + + /** + * Return the imaginary parts of the eigenvalues + * @return imag(diag(D)) + */ + public double[] getImagEigenvalues(){ + return e; + } + + /** + * Return the block diagonal eigenvalue matrix + * @return D + */ + public Matrix getD(){ + int i, j; + Matrix D = new Matrix(n,n); + for(i=0; i 0) + D.set(i, i+1, e[i]); + else if (e[i] < 0) + D.set(i, i-1, e[i]); + } + return D; + } +} diff --git a/src/math/matrix/LUDecomposition.java b/src/math/matrix/LUDecomposition.java new file mode 100644 index 0000000..bf277ec --- /dev/null +++ b/src/math/matrix/LUDecomposition.java @@ -0,0 +1,282 @@ +package math.matrix; + +import java.io.Serializable; + +/** + * LU Decomposition. + *

+ * For an m-by-n matrix A with m ≥ n, the LU decomposition is an m-by-n + * unit lower triangular matrix L, an n-by-n upper triangular matrix U, + * and a permutation vector piv of length m so that A(piv,:) = L*U. + * If m < n, then L is m-by-m and U is m-by-n. + *

+ * The LU decompostion with pivoting always exists, even if the matrix is + * singular, so the constructor will never fail. The primary use of the + * LU decomposition is in the solution of square systems of simultaneous + * linear equations. This will fail if isNonsingular() returns false. + */ +public class LUDecomposition implements Serializable { + + /** + * UID + */ + private static final long serialVersionUID = -3852210156107961377L; + + /** + * Array for internal storage of decomposition. + * @serial internal array storage + */ + private double[][] LU; + + /** + * @serial row dimension + */ + private int m; + + /** + * @serial column dimension + */ + private int n; + + /** + * @serial pivot sign + */ + private int pivsign; + + /** + * Internal storage of pivot vector. + * @serial pivot vector. + */ + private int[] piv; + + /** + * LU Decomposition + * Structure to access L, U and piv. + * @param A Rectangular matrix + */ + public LUDecomposition(Matrix A){ + int i, j, k; + + // Use a "left-looking", dot-product, Crout/Doolittle algorithm. + LU = A.getCopy(); + m = A.getM(); + n = A.getN(); + piv = new int[m]; + for(i=0; i Math.abs(LUcolj[p])) + p = i; + if(p!=j){ + for(k=0; k + This constructor computes L and U with the "daxpy"-based elimination + algorithm used in LINPACK and MATLAB. In Java, we suspect the dot-product, + Crout algorithm will be faster. We have temporarily included this + constructor until timing experiments confirm this suspicion. +

+ Structure to access L, U and piv. + @param A Rectangular matrix + @param linpackflag Use Gaussian elimination. Actual value ignored. + */ + public LUDecomposition(Matrix A, int linpackflag){ + int i, j, k; + // Initialize. + LU = A.getCopy(); + m = A.getM(); + n = A.getN(); + piv = new int[m]; + for(i=0; i Math.abs(LU[p][k])) + p = i; + // Exchange if necessary. + if(p != k) { + for(j=0; jj) + L[i][j] = LU[i][j]; + else if(i==j) + L[i][j] = 1.0; + else + L[i][j] = 0.0; + } + } + return X; + } + + /** + * Return upper triangular factor + * @return U + */ + public Matrix getU(){ + int i, j; + Matrix X = new Matrix(n,n); + double[][] U = X.getD(); + for(i=0; i=0; k--) + for(j=0; ji + * @throws NoSquareException + * @throws IllegalDimensionException + * @throws ComplexException + */ + public Vector2D eigenvaluesRe(int analytic) throws NoSquareException, IllegalDimensionException, ComplexException{ + // TODO: also Imaginary + if(m!=n && m!=2) throw new IllegalDimensionException("Illegal matrix dimensions."); + // A = [ a, b; c, d] + // det| a-l, b; c, d-l| = l^2 - tr(A)*l + det (A) + Polynomial p = new Polynomial(1,-trace(),det(analytic)); + return new Vector2D(p.rootsRe()); + } + + /** + * @return eigenvectors xi in matrix Λ + * @throws IllegalDimensionException + * @throws NoSquareException + * @throws ComplexException + */ + public Matrix2D eigenvectors(int analytic) throws NoSquareException, IllegalDimensionException, ComplexException{ + Vector lambda = eigenvaluesRe(analytic); + Matrix2D S = new Matrix2D(); + if(isR2()){ + Matrix2D Al; + Vector2D x = new Vector2D(); + int i; + double a1, a2; + for(i=0;i<2;i++){ + Al = (Matrix2D) this.minus(Matrix.identity(2).times(lambda.get(i))); + a2 = Al.get(1, 0); + a1 = -Al.get(1, 1); // to the other side: a1*x1 = -a2*x2 + if(a1 == 0 && a2 == 0){ + a2 = Al.get(0, 0); + a1 = -Al.get(0, 1); // to the other side: a1*x1 = -a2*x2 + if(a2 > 0){ + x.set(0, a1); + x.set(1, a2); + } else { + x.set(0, -a1); + x.set(1, -a2); + } + } else + if(a2 > 0){ + x.set(0, a1); + x.set(1, a2); + } else { + x.set(0, -a1); + x.set(1, -a2); + } + S.setN(i, x.normalize()); + } + } + return S; + } + + public void plot() throws IllegalDimensionException{ + if(!isR2()) throw new IllegalDimensionException(); + StdDraw sd = new StdDraw(); + double maxM0 = Maths.sumAbs(getM(0).getArray()); + double maxM1 = Maths.sumAbs(getM(1).getArray()); + double maxN0 = Maths.sumAbs(getN(0).getArray()); + double maxN1 = Maths.sumAbs(getN(1).getArray()); + double max = Maths.max(maxM0,maxM1,maxN0,maxN1); + sd.setScale(-max, max); + StdDraw.polygon( + new double[]{0,get(0,0),get(0,0)+get(0,1),get(0,1)}, + new double[]{0,get(1,0),get(1,0)+get(1,1),get(1,1)}); + StdDraw.filledPolygon( + new double[]{0,get(0,0),get(0,0)+get(0,1),get(0,1)}, + new double[]{0,get(1,0),get(1,0)+get(1,1),get(1,1)}); + } + + public static void main(String[] args) { + Matrix2D a = new Matrix2D(0,1,1,0); + try { + Matrix b = new Matrix2D(new double[][]{{0,-1},{1,0}}).times(a); + System.out.println(b); +// System.out.println((Matrix2D) b); + + Matrix2D P = new Matrix2D(new double[][]{ { 1 }, { 1 } }); + System.out.println(P.setName("P")); + System.out.println("rot(P,45°) = \n" + P.rotate(45)); + + Vector2D v = new Vector2D(1,1); + Vector2D w = new Vector2D(0,1); + System.out.println(Diagonal2D.identity(v.n())); + P = (Matrix2D) Matrix2D.mirror(v); + System.out.println("Pw = "); + P.times(w).println(); + + Matrix A = new Matrix2D(new double[][]{ { 2, 1 }, { 0, 1 } }); + System.out.println(A.setName("A")); + Matrix B = A.transpose().times(A); + System.out.println("A'A = \n" + B); + Matrix X = B.eigenvectors(); + System.out.println("X = \n" + X); + System.out.println("X'A'A = \n" + X.transpose().times(B)); + Matrix L2 = X.transpose().times(B).times(X); + System.out.println("L^2 = X'A'AX = \n" + L2); + Matrix L = L2.times(new Matrix(new double[][]{ { 1/Math.sqrt(L2.get(0, 0)), 0 }, { 0, 1/Math.sqrt(L2.get(1, 1)) } })); + System.out.println("L = \n" + L); + Matrix S = X.times(L).times(X); + System.out.println("S = XLX' = \n" + S); + Matrix2D Q = (Matrix2D) A.times(S.inverse()); + System.out.println("Q = AS^-1 = \n" + Q); + Q.plot(); + } catch (IllegalDimensionException | SingularException e) { + e.printStackTrace(); + } + } + +} diff --git a/src/math/matrix/Matrix3D.java b/src/math/matrix/Matrix3D.java new file mode 100644 index 0000000..2b02f95 --- /dev/null +++ b/src/math/matrix/Matrix3D.java @@ -0,0 +1,193 @@ +package math.matrix; + +import math.Maths; +import math.equation.Polynomial; +import exception.ComplexException; +import exception.IllegalDimensionException; +import exception.NoSquareException; + + +public class Matrix3D extends Matrix{ + /** + * UID + */ + private static final long serialVersionUID = 2505785127150028971L; + + /** + * Create an null matrix. + */ + public Matrix3D(){ + super(); + } + + /** + * Create an 3-by-n matrix with zeros. + * @param n number of columns + */ + public Matrix3D(int n){ + super(3,n); + } + + /** + * Create an 3-by-3 matrix based on components. + * @param a11 + * @param a12 + * @param a13 + * @param a21 + * @param a22 + * @param a23 + * @param a31 + * @param a32 + * @param a33 + */ + public Matrix3D(double a11, double a12, double a13, + double a21, double a22, double a23, + double a31, double a32, double a33){ + super(3,3); + set(0,0,a11); + set(0,1,a12); + set(0,2,a13); + set(1,0,a21); + set(1,1,a22); + set(1,2,a23); + set(2,0,a31); + set(2,1,a32); + set(2,2,a33); + } + + /** + * Create an 3-by-3 matrix based on a 2d array. + * @param a 3d double array + */ + public Matrix3D(double[][] a){ + super(a); + } + + public Matrix3D(Matrix B){ + super(B); + } + + /** + * 3d rotation around unit vector u by ϑ degrees counterclockwise or + * clockwise with negative degrees for a right-hand system. + * @param ϑ in degrees + * @param u unit vector + * @return 3d rotated matrix + * @throws IllegalDimensionException + */ + public Matrix3D rotate(double ϑ, double[] u) throws IllegalDimensionException{ + if(!isR3() && u.length == 3) throw new IllegalDimensionException("Illegal matrix dimensions."); + // TODO: verify + double u1 = u[0], u2 = u[1], u3 = u[2]; + return (Matrix3D) new Matrix3D(new double[][]{ + {u1*u1*(1-Maths.cosd(ϑ))+Maths.cosd(ϑ), + u1*u2*(1-Maths.cosd(ϑ))-u3*Maths.sind(ϑ), + u1*u3*(1-Maths.cosd(ϑ))+u2*Maths.sind(ϑ)}, + {u2*u1*(1-Maths.cosd(ϑ))+u3*Maths.sind(ϑ), + u2*u2*(1-Maths.cosd(ϑ))+Maths.cosd(ϑ), + u2*u3*(1-Maths.cosd(ϑ))-u1*Maths.sind(ϑ)}, + {u3*u1*(1-Maths.cosd(ϑ))-u2*Maths.sind(ϑ), + u3*u2*(1-Maths.cosd(ϑ))+u1*Maths.sind(ϑ), + u3*u3*(1-Maths.cosd(ϑ))+Maths.cosd(ϑ)}}).times(this); + } + + /** + * Rotation about the x axis by φ degrees counterclockwise for a right-hand system. + * @param φ in degrees + * @return 3d rotated matrix + * @throws IllegalDimensionException + */ + public Matrix3D rotateX(double φ) throws IllegalDimensionException{ + if(!isR3()) throw new IllegalDimensionException("Illegal matrix dimensions."); + return (Matrix3D) new Matrix3D(new double[][]{ + {1,0,0}, + {0,Maths.cosd(φ),Maths.sind(φ)}, + {0,-Maths.sind(φ),Maths.cosd(φ)}}).times(this); + } + + /** + * Rotation about the y axis by ϑ degrees counterclockwise for a right-hand system. + * @param ϑ in degrees + * @return 3d rotated matrix + * @throws IllegalDimensionException + */ + public Matrix3D rotateY(double ϑ) throws IllegalDimensionException{ + if(!isR3()) throw new IllegalDimensionException("Illegal matrix dimensions."); + return (Matrix3D) new Matrix3D(new double[][]{ + {Maths.cosd(ϑ),0,-Maths.sind(ϑ)}, + {0,1,0}, + {Maths.sind(ϑ),0,Maths.cosd(ϑ)}}).times(this); + } + + /** + * Rotation about the z axis by ψ degrees counterclockwise for a right-hand system. + * @param ψ in degrees + * @return 3d rotated matrix + * @throws IllegalDimensionException + */ + public Matrix3D rotateZ(double ψ) throws IllegalDimensionException{ + if(!isR3()) throw new IllegalDimensionException("Illegal matrix dimensions."); + return (Matrix3D) new Matrix3D(new double[][]{ + {Maths.cosd(ψ),Maths.sind(ψ),0}, + {-Maths.sind(ψ),Maths.cosd(ψ),0}, + {0,0,1}}).times(this); + } + + /** + * 3d rotation by φ around x axis, ϑ around y axis and ψ degrees around z axis counterclockwise or + * clockwise with negative degrees for a right-hand system. + * @param φ in degrees + * @param ϑ in degrees + * @param ψ in degrees + * @return 3d rotated matrix + * @throws IllegalDimensionException + */ + public Matrix3D rotateXYZ(double φ, double ϑ, double ψ) throws IllegalDimensionException{ + if(!isR3()) throw new IllegalDimensionException("Illegal matrix dimensions."); + return (Matrix3D) rotateX(φ).times(rotateY(ϑ)).times(rotateZ(ψ)).times(this); + } + + /** + * Projection onto the x-y plane + * @return matrix projected onto the x-y plane + * @throws IllegalDimensionException + */ + public Matrix3D projectionToXY() throws IllegalDimensionException{ + if(!isR3()) throw new IllegalDimensionException("Illegal matrix dimensions."); + return (Matrix3D) new Matrix(new double[][]{{1,0,0},{0,1,0},{0,0,0}}).times(this); + } + + /** + * Real eigenvalue calculated with the root formula of second and third degree. + * @return eigenvalues λi + * @throws NoSquareException + * @throws IllegalDimensionException + * @throws ComplexException + */ + public Vector3D eigenvaluesRe(int analytic) throws NoSquareException, IllegalDimensionException, ComplexException{ + if(m!=n && m!=3) throw new IllegalDimensionException("Illegal matrix dimensions."); + Vector3D lambda = new Vector3D(m); + // matrix diagonal? + double a12 = get(0, 1); + double a13 = get(0, 2); + double a23 = get(1, 2); + double d = a12*a12 + a13*a13 + a23*a23; + if (d == 0){ // A is diagonal. + lambda = new Vector3D(get(0,0),get(1,1),get(2,2)); + lambda.sort(); // don't have to but ... + } else { + Polynomial p = new Polynomial(1, + -trace(), + -((this.times(this)).trace()-trace()*trace())/2, + -det(analytic)); + lambda = new Vector3D(p.rootsRe()); + } + return lambda; + } + + public static void main(String[] args) { + Matrix3D a = new Matrix3D(0,0,1,0,1,0,1,0,0); + System.out.println(a); + } + +} diff --git a/src/math/matrix/QRDecomposition.java b/src/math/matrix/QRDecomposition.java new file mode 100644 index 0000000..4efa362 --- /dev/null +++ b/src/math/matrix/QRDecomposition.java @@ -0,0 +1,197 @@ +package math.matrix; + +import java.io.Serializable; + +import math.Maths; + +/** QR Decomposition. + *

+ * For an m-by-n matrix A with m ≥ n, the QR decomposition is an m-by-n + * orthogonal matrix Q and an n-by-n upper triangular matrix R so that + * A = Q*R. + *

+ * The QR decompostion always exists, even if the matrix does not have + * full rank, so the constructor will never fail. The primary use of the + * QR decomposition is in the least squares solution of nonsquare systems + * of simultaneous linear equations. This will fail if isFullRank() + * returns false. + */ +public class QRDecomposition implements Serializable { + + /** + * UID + */ + private static final long serialVersionUID = 6913030878857169502L; + +/** + * Array for internal storage of decomposition. + * @serial internal array storage. + */ + private double[][] QR; + + /** + * Row dimensions. + * @serial row dimension. + */ + private int m; + + /** + * Column dimensions. + * @serial column dimension. + */ + private int n; + + /** + * Array for internal storage of diagonal of R. + * @serial diagonal of R. + */ + private double[] Rdiag; + + /** + * QR Decomposition, computed by Householder reflections. + * Structure to access R and the Householder vectors and compute Q. + * @param A rectangular matrix + */ + public QRDecomposition(Matrix A){ + int i, j, k; + // Initialize. + QR = A.getCopy(); + m = A.getM(); + n = A.getN(); + Rdiag = new double[n]; + // Main loop. + for(k=0; k= j) + H.set(i,j, QR[i][j]); + else + H.set(i,j, 0.0); + } + return H; + } + + /** + * Return the upper triangular factor + * @return R + */ + public Matrix getR(){ + int i, j; + Matrix R = new Matrix(n,n); + for(i=0; i=0; k--) { + for(i=0; i=0; k--){ + for(j=0; j + * For an m-by-n matrix A with m ≥ n, the singular value decomposition is + * an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and + * an n-by-n orthogonal matrix V so that A = U*S*V'. + *

+ * The singular values, σk = Skk, are ordered so that + * σ1 ≥ σ2 ≥ ... ≥ σn. + *

+ * The singular value decompostion always exists, so the constructor will + * never fail. The matrix condition number and the effective numerical + * rank can be computed from this decomposition. + */ +public class SingularValueDecomposition implements Serializable { + + /** + * UID + */ + private static final long serialVersionUID = -4816488865370991490L; + +/** + * Array for internal storage of U. + * @serial internal storage of U. + */ + private double[][] U; + + /** + * Array for internal storage of V. + * @serial internal storage of V. + */ + private double[][] V; + + /** + * Array for internal storage of singular values. + * @serial internal storage of singular values. + */ + private double[] s; + + /** + * Row dimensions. + * @serial row dimension. + */ + private int m; + + /** + * Column dimensions. + * @serial column dimension. + */ + private int n; + +/* ------------------------ + Constructor + * ------------------------ */ + + /** + * Construct the singular value decomposition Structure to access U, S and V. + * @param Arg rectangular matrix + */ + public SingularValueDecomposition(Matrix Arg){ + int i, j, k; + // Derived from LINPACK code. + // Initialize. + double[][] A = Arg.getCopy(); + m = Arg.getM(); + n = Arg.getN(); + /* TODO: Apparently the failing cases are only a proper subset of (m= n"); } + */ + int nu = Math.min(m,n); + s = new double [Math.min(m+1,n)]; + U = new double [m][nu]; + V = new double [n][n]; + double[] e = new double [n]; + double[] work = new double [m]; + boolean wantu = true; + boolean wantv = true; + // Reduce A to bidiagonal form, storing the diagonal elements + // in s and the super-diagonal elements in e. + int nct = Math.min(m-1,n); + int nrt = Math.max(0,Math.min(n-2,m)); + for(k=0; k=0; k--){ + if(s[k] != 0.0) { + for(j=k+1; j=0; k--){ + if((k < nrt) & (e[k] != 0.0)) + for(j=k+1; j 0){ + int kase; + // Here is where a test for too many iterations would go. + // + // This section of the program inspects for + // negligible elements in the s and e arrays. On + // completion the variables kase and k are set as follows. + // kase = 1 if s(p) and e[k-1] are negligible and k

=-1; k--){ + if(k == -1) break; + if(Math.abs(e[k]) <= + tiny + eps*(Math.abs(s[k]) + Math.abs(s[k+1]))){ + e[k] = 0.0; + break; + } + } + if(k == p-2) + kase = 4; + else { + int ks; + for(ks=p-1; ks>=k; ks--){ + if (ks == k) break; + double t = (ks != p ? Math.abs(e[ks]) : 0.) + + (ks != k+1 ? Math.abs(e[ks-1]) : 0.); + if(Math.abs(s[ks]) <= tiny + eps*t){ + s[ks] = 0.0; + break; + } + } + if(ks == k) + kase = 3; + else if(ks == p-1) + kase = 1; + else { + kase = 2; + k = ks; + } + } + k++; + // Perform the task indicated by kase. + switch (kase) { + // Deflate negligible s(p). + case 1: { + double f = e[p-2]; + e[p-2] = 0.0; + for(j=p-2; j>=k; j--){ + double t = Maths.hypot(s[j],f); + double cs = s[j]/t; + double sn = f/t; + s[j] = t; + if(j != k){ + f = -sn*e[j-1]; + e[j-1] = cs*e[j-1]; + } + if(wantv) + for(i=0; i= s[k+1]) + break; + double t = s[k]; + s[k] = s[k+1]; + s[k+1] = t; + if(wantv && (k < n-1)){ + for(i=0; i tol) + r++; + return r; + } +} diff --git a/src/math/matrix/TensorII.java b/src/math/matrix/TensorII.java new file mode 100644 index 0000000..5fe3a6e --- /dev/null +++ b/src/math/matrix/TensorII.java @@ -0,0 +1,13 @@ +package math.matrix; + +public class TensorII extends Matrix3D { + + /** + * UID + */ + private static final long serialVersionUID = -5664678448598889794L; + + public static void main(String[] args) { + } + +} diff --git a/src/math/matrix/TensorII2D.java b/src/math/matrix/TensorII2D.java new file mode 100644 index 0000000..63b3cbd --- /dev/null +++ b/src/math/matrix/TensorII2D.java @@ -0,0 +1,58 @@ +package math.matrix; + +import math.Maths; +import exception.IllegalDimensionException; + +/** + * Tensor 2x2 + * @author Daniel + * + */ +public class TensorII2D extends Matrix2D { + + /** + * UID + */ + private static final long serialVersionUID = -3195354821730519791L; + + public TensorII2D() { + super(2); + } + + public TensorII2D(double a11, double a12, double a21, double a22) { + super(a11, a12, a21, a22); + } + + public TensorII2D(double[][] a) { + super(a); + } + + public TensorII2D(Matrix B) { + super(B); + } + + /** + * @param ϑ in degrees + * @return 2d rotation matrix + */ + public static TensorII2D rotation(double ϑ){ + return new TensorII2D(new double[][]{ + {Maths.cosd(ϑ),-Maths.sind(ϑ)}, + {Maths.sind(ϑ), Maths.cosd(ϑ)}}); + } + + public TensorII2D rotate(double ϑ){ + try { + return (TensorII2D) rotation(ϑ).times(this).times(rotation(ϑ).transpose()); + } catch (IllegalDimensionException e) { + e.printStackTrace(); + return null; + } + } + + public static void main(String[] args) { + TensorII2D tau = new TensorII2D(1, -1, -1, 1); + double ϑ = 45; + tau.rotate(ϑ).println(); + } +} diff --git a/src/math/matrix/Vector.java b/src/math/matrix/Vector.java new file mode 100644 index 0000000..5d56c17 --- /dev/null +++ b/src/math/matrix/Vector.java @@ -0,0 +1,899 @@ +package math.matrix; + +import java.io.Serializable; +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; + +import math.Maths; +import math.geometry.Point; +import thisandthat.WObject; +import exception.IllegalDimensionException; + +/** + * Vector of real numbers, v∈ℝ. + * @author Daniel Weschke + */ +public class Vector extends WObject implements Cloneable, Serializable{ + + /** + * UID + */ + private static final long serialVersionUID = 7909642082289701909L; + + protected List vector; +// protected double[] vector; + + public Vector(){ + } + + /** + * Create zero vector + * @param n + */ + public Vector(int n){ + int i; +// vector = new double[n]; + vector = new ArrayList(); + for(i=0; i data){ + this(data.size()); + int i; + int n = data.size(); + for(i=0; i get(){ + return vector; + } + + /** + * @param i index + * @return corresponding coordinate + */ + public double get(int i){ + return vector.get(i); + } + + /** + * @return vector array (integers) + */ + public int[] getInts(){ + int i; + int[] result = new int[n()]; + for(i=0; i0?get(0):Double.NaN; + } + + /** + * get y (second element) + * @return y value + */ + public double y(){ + return n()>1?get(1):Double.NaN; + } + + /** + * get z (third element) + * @return z value + */ + public double z(){ + return n()>2?get(2):Double.NaN; + } + + /** + * length of vector + * @return dimension + */ + public int n(){ + return length(); + } + + /** + * length of vector + * @return dimension + */ + public int length(){ +// return vector.length; + return vector.size(); + } + + /** + * Create vector with zeros + * @param n size + * @return zero vector + */ + public static Vector zeros(int n){ + return new Vector(n); + } + + /** + * Create vector with ones + * @param n size + * @return one vector + */ + public static Vector ones(int n){ + return fill(n, 1); + } + + /** + * Create vector with given scalar + * @param n size + * @param s scalar + * @return filled vector + */ + public static Vector fill(int n, double s){ + Vector a = new Vector(n); + int i; + for(i=0; i0 ? 1 : -1; + return fill(start, increment, end); + } + + /** + * Create vector and fill it from a start value to an end value with n equidistant steps + * @param start + * @param end + * @param n steps + * @return filled vector [start,(end-start)*1/n,(end-start)*2/n,...,end] + */ + public static Vector fillN(double start, double end, int n){ + double increment = (end-start)/n; + return fill(start, increment, end); + } + + /** + * Create vector and fill it from a start value and increments it to an end value + * @param start + * @param increment + * @param end + * @return filled vector [start,start+increment,start+2*increment,...,end] + */ + public static Vector fill(double start, double increment, double end){ + double delta = end-start; + int n = (int)(delta/increment)+1; + if((delta>=0&&increment<0)||(delta<0&&increment>=0)||increment==0) return new Vector(); + Vector a = new Vector(n); + int i; + for(i=0; iC = A • B = ATB + * @throws IllegalDimensionException Inner matrix dimensions must agree. + */ + public Vector times(Matrix B) throws IllegalDimensionException{ + int i,j; + if(n() != B.getM()) throw new IllegalDimensionException("Inner matrix dimensions must agree."); + double[] c = new double[n()]; + for(i=0; iTx. + * A • B = a1b1 + a2b2 + + * … + anbn; + * a scalar quantity + * @return scalar + * @throws IllegalDimensionException + */ + public double dot() throws IllegalDimensionException{ + return dot(this); + } + + /** + * Dot / inner / scalar product xTy. + * A • B = a1b1 + a2b2 + + * … + anbn; + * a scalar quantity + * @param a vector + * @return scalar + * @throws IllegalDimensionException + */ + public double dot(Vector a) throws IllegalDimensionException{ + if(n() != a.n()) throw new IllegalDimensionException("Illegal vector dimension."); + int i; + double c = 0; + for(i=0; ii = aibi + * @param b vector + * @return vector + */ + public Vector timesE(Vector b){ + int i; + Vector c = create(n()); + for(i=0; ii = ai/bi + * @param b vector + * @return vector + */ + public Vector overE(Vector b){ + int i; + Vector c = create(n()); + for(i=0; ia/||a|| + * @return unit vector + */ + public Vector normalize(){ + Vector c = times(1.0/norm()); + this.vector = c.vector; + return this; + } + + /** + * @return the corresponding unit vector + */ + public Vector direction() { + if(norm() == 0.0) throw new RuntimeException("Zero-vector has no direction"); + return times(1.0 / norm()); + } + + /** + * Element-wise cosine of argument in radians. + * @return the cosine for each element of the matrix + */ + public Vector cos() { + int i; + int n = n(); + Vector c = new Vector(n); + for(i=0; i1|, |v2|, …, |vn| } + * @return vector with absolute values + */ + public Vector abs(){ + int i; + Vector c = create(n()); + for(i=0; i=0 ? get(i) : -get(i)); + return c; + } + + /** + * Reflection against a plane due to a vector + * @return reflected vector + * @throws IllegalDimensionException + */ + public Vector flip(Vector v) throws IllegalDimensionException{ + // P = I - alpha v v'; alpha = 2 / (v'v) + return Matrix.mirror(v).times(this); + } + + /** + * Maximum value of the vector + * @return maximum value in vector, -∞ if no such value. + */ + public double max(){ + return Maths.max(getArray()); + } + + /** + * Minimum value of the vector + * @return minimum value in vector, +∞ if no such value. + */ + public double min(){ + return Maths.min(getArray()); + } + + /** + * Find value in vector. + * @param a value + * @return index + */ + public int indexOf(double a){ + int i; + for(i=0; i= 0; i--, j++) + tmp.set(j, get(i)); + vector = tmp.vector; + return this; + } + + /** + * Zero vector? + * @return boolean + */ + public boolean isZero(){ + int i; + boolean result = true; + for(i=0; i 0){ + result = false; + break; + } + } + return result; + } + + public boolean isNormalized(){ + return Math.abs(norm()-1)2. + * n=2 + * @return true or false + */ + public boolean isR2(){ + return (n()==2); + } + + /** + * Vector in euclidean plane ℝ3. + * n=3 + * @return true or false + */ + public boolean isR3(){ + return (n()==3); + } + + /** + * + * @param a int array + * @param b int array + * @return true if two integer arrays have same length and + * all corresponding pairs of integers are equal + */ + public static boolean equals(int[] a, int[] b){ + if(a.length != b.length) return false; // same length? + int i; + for(i=0; i2b3 - a3b2, + * a3b1 - a1b3, + * a1b2 - a2b1)T; + * a vector quantity + * @param a vector + * @return vector + * @throws IllegalDimensionException + */ + public Vector3D cross(Vector3D a){ + Vector3D c = new Vector3D(); + c.set(0, get(1)*a.get(2) - get(2)*a.get(1)); + c.set(1, get(2)*a.get(0) - get(0)*a.get(2)); + c.set(2, get(0)*a.get(1) - get(1)*a.get(0)); + return c; + } + + /** + * Scalar triple product. + * A • (B x C); a scalar quantity + * @param b vector + * @param c vector + * @return scalar + * @throws IllegalDimensionException + */ + public double scalTrip(Vector3D b, Vector3D c){ + try { + return dot(b.cross(c)); + } catch (IllegalDimensionException e) { + e.printStackTrace(); + return Double.NaN; + } + } + + /** + * Vector triple product. + * A x (B x C); a vector quantity + * @param b vector + * @param c vector + * @return vector + * @throws IllegalDimensionException + */ + public Vector3D vecTrip(Vector3D b, Vector3D c) { + return this.cross(b.cross(c)); + } + + /** + * cosine between this and given vector + * @param a vector + * @return cosθ + */ + public double cos(Vector3D a){ + try { + return dot(a)/(norm()*a.norm()); + } catch (IllegalDimensionException e) { + e.printStackTrace(); + return Double.NaN; + } + } + + /** + * sine between this and given vector + * @param a vector + * @return sinθ + */ + public double sin(Vector3D a){ + return cross(a).norm()/(norm()*a.norm()); + } + + public static void main(String[] args) { + Vector3D u = new Vector3D(3,3,0); + Vector3D v = new Vector3D(0,2,2); + System.out.println("cos = " + u.cos(v) + " : " + (u.cos(v)==0.5)); + System.out.println("|w| = " + u.cross(v).norm() + " : " + + ((u.cross(v).norm())==(6*Math.sqrt(3)))); + } + +} diff --git a/src/stdlib/StdArrayIO.java b/src/stdlib/StdArrayIO.java new file mode 100644 index 0000000..7fa4f1e --- /dev/null +++ b/src/stdlib/StdArrayIO.java @@ -0,0 +1,256 @@ +package stdlib; + +/************************************************************************* + * Compilation: javac StdArrayIO.java + * Execution: java StdArrayIO < input.txt + * + * A library for reading in 1D and 2D arrays of integers, doubles, + * and booleans from standard input and printing them out to + * standard output. + * + * % more tinyDouble1D.txt + * 4 + * .000 .246 .222 -.032 + * + * % more tinyDouble2D.txt + * 4 3 + * .000 .270 .000 + * .246 .224 -.036 + * .222 .176 .0893 + * -.032 .739 .270 + * + * % more tinyBoolean2D.txt + * 4 3 + * 1 1 0 + * 0 0 0 + * 0 1 1 + * 1 1 1 + * + * % cat tinyDouble1D.txt tinyDouble2D.txt tinyBoolean2D.txt | java StdArrayIO + * 4 + * 0.00000 0.24600 0.22200 -0.03200 + * + * 4 3 + * 0.00000 0.27000 0.00000 + * 0.24600 0.22400 -0.03600 + * 0.22200 0.17600 0.08930 + * 0.03200 0.73900 0.27000 + * + * 4 3 + * 1 1 0 + * 0 0 0 + * 0 1 1 + * 1 1 1 + * + *************************************************************************/ + + +/** + * Standard array IO. This class provides methods for reading + * in 1D and 2D arrays from standard input and printing out to + * standard output. + *

+ * For additional documentation, see + * Section 2.2 of + * Introduction to Programming in Java: An Interdisciplinary Approach + * by Robert Sedgewick and Kevin Wayne. + * + * @author Robert Sedgewick + * @author Kevin Wayne + */ +public class StdArrayIO { + + // it doesn't make sense to instantiate this class + private StdArrayIO() { } + + /** + * Read in and return an array of doubles from standard input. + */ + public static double[] readDouble1D() { + int N = StdIn.readInt(); + double[] a = new double[N]; + for (int i = 0; i < N; i++) { + a[i] = StdIn.readDouble(); + } + return a; + } + + /** + * Print an array of doubles to standard output. + */ + public static void print(double[] a) { + int N = a.length; + StdOut.println(N); + for (int i = 0; i < N; i++) { + StdOut.printf("%9.5f ", a[i]); + } + StdOut.println(); + } + + + /** + * Read in and return an M-by-N array of doubles from standard input. + */ + public static double[][] readDouble2D() { + int M = StdIn.readInt(); + int N = StdIn.readInt(); + double[][] a = new double[M][N]; + for (int i = 0; i < M; i++) { + for (int j = 0; j < N; j++) { + a[i][j] = StdIn.readDouble(); + } + } + return a; + } + + /** + * Print the M-by-N array of doubles to standard output. + */ + public static void print(double[][] a) { + int M = a.length; + int N = a[0].length; + StdOut.println(M + " " + N); + for (int i = 0; i < M; i++) { + for (int j = 0; j < N; j++) { + StdOut.printf("%9.5f ", a[i][j]); + } + StdOut.println(); + } + } + + + /** + * Read in and return an array of ints from standard input. + */ + public static int[] readInt1D() { + int N = StdIn.readInt(); + int[] a = new int[N]; + for (int i = 0; i < N; i++) { + a[i] = StdIn.readInt(); + } + return a; + } + + /** + * Print an array of ints to standard output. + */ + public static void print(int[] a) { + int N = a.length; + StdOut.println(N); + for (int i = 0; i < N; i++) { + StdOut.printf("%9d ", a[i]); + } + StdOut.println(); + } + + + /** + * Read in and return an M-by-N array of ints from standard input. + */ + public static int[][] readInt2D() { + int M = StdIn.readInt(); + int N = StdIn.readInt(); + int[][] a = new int[M][N]; + for (int i = 0; i < M; i++) { + for (int j = 0; j < N; j++) { + a[i][j] = StdIn.readInt(); + } + } + return a; + } + + /** + * Print the M-by-N array of ints to standard output. + */ + public static void print(int[][] a) { + int M = a.length; + int N = a[0].length; + StdOut.println(M + " " + N); + for (int i = 0; i < M; i++) { + for (int j = 0; j < N; j++) { + StdOut.printf("%9d ", a[i][j]); + } + StdOut.println(); + } + } + + + /** + * Read in and return an array of booleans from standard input. + */ + public static boolean[] readBoolean1D() { + int N = StdIn.readInt(); + boolean[] a = new boolean[N]; + for (int i = 0; i < N; i++) { + a[i] = StdIn.readBoolean(); + } + return a; + } + + /** + * Print an array of booleans to standard output. + */ + public static void print(boolean[] a) { + int N = a.length; + StdOut.println(N); + for (int i = 0; i < N; i++) { + if (a[i]) StdOut.print("1 "); + else StdOut.print("0 "); + } + StdOut.println(); + } + + /** + * Read in and return an M-by-N array of booleans from standard input. + */ + public static boolean[][] readBoolean2D() { + int M = StdIn.readInt(); + int N = StdIn.readInt(); + boolean[][] a = new boolean[M][N]; + for (int i = 0; i < M; i++) { + for (int j = 0; j < N; j++) { + a[i][j] = StdIn.readBoolean(); + } + } + return a; + } + + /** + * Print the M-by-N array of booleans to standard output. + */ + public static void print(boolean[][] a) { + int M = a.length; + int N = a[0].length; + StdOut.println(M + " " + N); + for (int i = 0; i < M; i++) { + for (int j = 0; j < N; j++) { + if (a[i][j]) StdOut.print("1 "); + else StdOut.print("0 "); + } + StdOut.println(); + } + } + + + /** + * Test client. + */ + public static void main(String[] args) { + + // read and print an array of doubles + double[] a = StdArrayIO.readDouble1D(); + StdArrayIO.print(a); + StdOut.println(); + + // read and print a matrix of doubles + double[][] b = StdArrayIO.readDouble2D(); + StdArrayIO.print(b); + StdOut.println(); + + // read and print a matrix of doubles + boolean[][] d = StdArrayIO.readBoolean2D(); + StdArrayIO.print(d); + StdOut.println(); + } + +} \ No newline at end of file diff --git a/src/stdlib/StdIn.java b/src/stdlib/StdIn.java new file mode 100644 index 0000000..baa6026 --- /dev/null +++ b/src/stdlib/StdIn.java @@ -0,0 +1,210 @@ +package stdlib; + +/************************************************************************* + * Compilation: javac StdIn.java + * Execution: java StdIn + * + * Reads in data of various types from standard input. + * + *************************************************************************/ + +import java.io.BufferedInputStream; +import java.util.Locale; +import java.util.Scanner; + +/** + * Standard input. This class provides methods for reading strings + * and numbers from standard input. + *

+ * The Locale used is: language = English, country = US. This is consistent + * with the formatting conventions with Java floating-point literals, + * command-line arguments (via Double.parseDouble()) + * and standard output (via System.out.print()). It ensures that + * standard input works with the input files used in the textbook. + *

+ * For additional documentation, see Section 1.5 of + * Introduction to Programming in Java: An Interdisciplinary Approach by Robert Sedgewick and Kevin Wayne. + */ +public final class StdIn { + + // assume Unicode UTF-8 encoding + private static String charsetName = "UTF-8"; + + // assume language = English, country = US for consistency with System.out. + private static Locale usLocale = new Locale("en", "US"); + + // the scanner object + private static Scanner scanner = new Scanner(new BufferedInputStream(System.in), charsetName); + + // static initializer + static { scanner.useLocale(usLocale); } + + // singleton pattern - can't instantiate + private StdIn() { } + + + /** + * Is there only whitespace left on standard input? + */ + public static boolean isEmpty() { + return !scanner.hasNext(); + } + + /** + * Return next string from standard input + */ + public static String readString() { + return scanner.next(); + } + + /** + * Return next int from standard input + */ + public static int readInt() { + return scanner.nextInt(); + } + + /** + * Return next double from standard input + */ + public static double readDouble() { + return scanner.nextDouble(); + } + + /** + * Return next float from standard input + */ + public static float readFloat() { + return scanner.nextFloat(); + } + + /** + * Return next short from standard input + */ + public static short readShort() { + return scanner.nextShort(); + } + + /** + * Return next long from standard input + */ + public static long readLong() { + return scanner.nextLong(); + } + + /** + * Return next byte from standard input + */ + public static byte readByte() { + return scanner.nextByte(); + } + + /** + * Return next boolean from standard input, allowing "true" or "1" for true, + * and "false" or "0" for false + */ + public static boolean readBoolean() { + String s = readString(); + if (s.equalsIgnoreCase("true")) return true; + if (s.equalsIgnoreCase("false")) return false; + if (s.equals("1")) return true; + if (s.equals("0")) return false; + throw new java.util.InputMismatchException(); + } + + /** + * Does standard input have a next line? + */ + public static boolean hasNextLine() { + return scanner.hasNextLine(); + } + + /** + * Return rest of line from standard input + */ + public static String readLine() { + return scanner.nextLine(); + } + + /** + * Return next char from standard input + */ + // a complete hack and inefficient - email me if you have a better + public static char readChar() { + // (?s) for DOTALL mode so . matches a line termination character + // 1 says look only one character ahead + // consider precompiling the pattern + String s = scanner.findWithinHorizon("(?s).", 1); + return s.charAt(0); + } + + /** + * Return rest of input from standard input + */ + public static String readAll() { + if (!scanner.hasNextLine()) return null; + + // reference: http://weblogs.java.net/blog/pat/archive/2004/10/stupid_scanner_1.html + return scanner.useDelimiter("\\A").next(); + } + + /** + * Read rest of input as array of ints + */ + public static int[] readInts() { + String[] fields = readAll().trim().split("\\s+"); + int[] vals = new int[fields.length]; + for (int i = 0; i < fields.length; i++) + vals[i] = Integer.parseInt(fields[i]); + return vals; + } + + /** + * Read rest of input as array of doubles + */ + public static double[] readDoubles() { + String[] fields = readAll().trim().split("\\s+"); + double[] vals = new double[fields.length]; + for (int i = 0; i < fields.length; i++) + vals[i] = Double.parseDouble(fields[i]); + return vals; + } + + /** + * Read rest of input as array of strings + */ + public static String[] readStrings() { + String[] fields = readAll().trim().split("\\s+"); + return fields; + } + + + + /** + * Unit test + */ + public static void main(String[] args) { + + System.out.println("Type a string: "); + String s = StdIn.readString(); + System.out.println("Your string was: " + s); + System.out.println(); + + System.out.println("Type an int: "); + int a = StdIn.readInt(); + System.out.println("Your int was: " + a); + System.out.println(); + + System.out.println("Type a boolean: "); + boolean b = StdIn.readBoolean(); + System.out.println("Your boolean was: " + b); + System.out.println(); + + System.out.println("Type a double: "); + double c = StdIn.readDouble(); + System.out.println("Your double was: " + c); + System.out.println(); + + } + +} diff --git a/src/stdlib/StdOut.java b/src/stdlib/StdOut.java new file mode 100644 index 0000000..aa611b8 --- /dev/null +++ b/src/stdlib/StdOut.java @@ -0,0 +1,230 @@ +package stdlib; + +/************************************************************************* + * Compilation: javac StdOut.java + * Execution: java StdOut + * + * Writes data of various types to standard output. + * + *************************************************************************/ + +import java.io.OutputStreamWriter; +import java.io.PrintWriter; +import java.io.UnsupportedEncodingException; +import java.util.Locale; + +/** + * Standard output. This class provides methods for writing strings + * and numbers to standard output. + *

+ * For additional documentation, see Section 1.5 of + * Introduction to Programming in Java: An Interdisciplinary Approach by Robert Sedgewick and Kevin Wayne. + */ +public final class StdOut { + + // force Unicode UTF-8 encoding; otherwise it's system dependent + private static final String UTF8 = "UTF-8"; + + // assume language = English, country = US for consistency with StdIn + private static final Locale US_LOCALE = new Locale("en", "US"); + + // send output here + private static PrintWriter out; + + // this is called before invoking any methods + static { + try { + out = new PrintWriter(new OutputStreamWriter(System.out, UTF8), true); + } + catch (UnsupportedEncodingException e) { System.out.println(e); } + } + + // singleton pattern - can't instantiate + private StdOut() { } + + // close the output stream (not required) + /** + * Close standard output. + */ + public static void close() { + out.close(); + } + + /** + * Terminate the current line by printing the line separator string. + */ + public static void println() { + out.println(); + } + + /** + * Print an object to standard output and then terminate the line. + */ + public static void println(Object x) { + out.println(x); + } + + /** + * Print a boolean to standard output and then terminate the line. + */ + public static void println(boolean x) { + out.println(x); + } + + /** + * Print a char to standard output and then terminate the line. + */ + public static void println(char x) { + out.println(x); + } + + /** + * Print a double to standard output and then terminate the line. + */ + public static void println(double x) { + out.println(x); + } + + /** + * Print a float to standard output and then terminate the line. + */ + public static void println(float x) { + out.println(x); + } + + /** + * Print an int to standard output and then terminate the line. + */ + public static void println(int x) { + out.println(x); + } + + /** + * Print a long to standard output and then terminate the line. + */ + public static void println(long x) { + out.println(x); + } + + /** + * Print a short to standard output and then terminate the line. + */ + public static void println(short x) { + out.println(x); + } + + /** + * Print a byte to standard output and then terminate the line. + */ + public static void println(byte x) { + out.println(x); + } + + /** + * Flush standard output. + */ + public static void print() { + out.flush(); + } + + /** + * Print an Object to standard output and flush standard output. + */ + public static void print(Object x) { + out.print(x); + out.flush(); + } + + /** + * Print a boolean to standard output and flush standard output. + */ + public static void print(boolean x) { + out.print(x); + out.flush(); + } + + /** + * Print a char to standard output and flush standard output. + */ + public static void print(char x) { + out.print(x); + out.flush(); + } + + /** + * Print a double to standard output and flush standard output. + */ + public static void print(double x) { + out.print(x); + out.flush(); + } + + /** + * Print a float to standard output and flush standard output. + */ + public static void print(float x) { + out.print(x); + out.flush(); + } + + /** + * Print an int to standard output and flush standard output. + */ + public static void print(int x) { + out.print(x); + out.flush(); + } + + /** + * Print a long to standard output and flush standard output. + */ + public static void print(long x) { + out.print(x); + out.flush(); + } + + /** + * Print a short to standard output and flush standard output. + */ + public static void print(short x) { + out.print(x); + out.flush(); + } + + /** + * Print a byte to standard output and flush standard output. + */ + public static void print(byte x) { + out.print(x); + out.flush(); + } + + /** + * Print a formatted string to standard output using the specified + * format string and arguments, and flush standard output. + */ + public static void printf(String format, Object... args) { + out.printf(US_LOCALE, format, args); + out.flush(); + } + + /** + * Print a formatted string to standard output using the specified + * locale, format string, and arguments, and flush standard output. + */ + public static void printf(Locale locale, String format, Object... args) { + out.printf(locale, format, args); + out.flush(); + } + + // This method is just here to test the class + public static void main(String[] args) { + + // write to stdout + StdOut.println("Test"); + StdOut.println(17); + StdOut.println(true); + StdOut.printf("%.6f\n", 1.0/7.0); + } + +}