12,300,841 members (62,355 online)

55.7K views
87 bookmarked
Posted

# Yet Another Math Parser (YAMP)

, 30 Sep 2012 CPOL
Constructing a fast math parser using Reflection to do numerics like Matlab.
 ```using System; using YAMP; using System.Runtime.InteropServices; namespace YAMP.Numerics { /// /// Eigenvalues and eigenvectors of a real matrix. /// 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.Multiply(D.Multiply(V.Transpose())) and /// V.Multiply(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.Multiply(V) equals V.Multiply(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 Eigenvalues { #region Members /// /// Row and column dimension (square matrix). /// int n; /// /// Symmetry flag. /// bool issymmetric; /// /// Arrays for internal storage of eigenvalues. /// double[] d, e; /// /// Array for internal storage of eigenvectors. /// double[][] V; /// /// Array for internal storage of nonsymmetric Hessenberg form. /// double[][] H; /// /// Working storage for nonsymmetric algorithm. /// double[] ort; #endregion // Class variables #region Private Methods // Symmetric Householder reduction to tridiagonal form. 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. var scale = 0.0; var 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]; } var f = d[i - 1]; var 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]; } var 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; var 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++) { var 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. 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; var f = 0.0; var tst1 = 0.0; var 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])); var m = l; while (m < n) { if (System.Math.Abs(e[m]) <= eps * tst1) break; m++; } // If m == l, d[l] is an eigenvalue, // otherwise, iterate. if (m > l) { var iter = 0; do { iter = iter + 1; // (Could check iteration count here.) // Compute implicit shift var g = d[l]; var p = (d[l + 1] - g) / (2.0 * e[l]); var r = NumericHelpers.Hypot(p, 1.0); if (p < 0) r = -r; d[l] = e[l] / (p + r); d[l + 1] = e[l] * (p + r); var dl1 = d[l + 1]; var h = g - d[l]; for (int i = l + 2; i < n; i++) d[i] -= h; f = f + h; // Implicit QL transformation. p = d[m]; var c = 1.0; var c2 = c; var c3 = c; var el1 = e[l + 1]; var s = 0.0; var 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 = NumericHelpers.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++) { var k = i; var 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. 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. var low = 0; var high = n - 1; for (int m = low + 1; m <= high - 1; m++) { // Scale column. var 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++) { var 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. double cdivr, cdivi; 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. 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 var nn = this.n; var n = nn - 1; var low = 0; var high = nn - 1; var eps = Math.Pow(2.0, -52.0); var 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 var 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 var 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++) { bool 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) { p = p / x; q = q / x; r = r / x; } } if (x == 0.0) break; 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 <= System.Math.Min(j, high); k++) { z = z + V[i][k] * H[k][j]; } V[i][j] = z; } } } #endregion // Private Methods #region Constructor /// /// Check for symmetry, then construct the eigenvalue decomposition /// /// Square matrix /// Structure to access D and V. public Eigenvalues(MatrixValue Arg) { var A = Arg.GetRealArray(); n = Arg.DimensionX; V = new double[n][]; for (int i = 0; i < n; i++) V[i] = new double[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][]; for (int i2 = 0; i2 < n; i2++) H[i2] = new double[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(); } } #endregion // Constructor #region Public Properties /// /// Return the real parts of the eigenvalues /// /// real(diag(D)) virtual public double[] RealEigenvalues { get { return d; } } /// /// Return the imaginary parts of the eigenvalues /// /// imag(diag(D)) virtual public double[] ImagEigenvalues { get { return e; } } /// /// Return the block diagonal eigenvalue matrix /// /// D virtual public MatrixValue D { get { var X = new MatrixValue(n, n); for (int i = 0; i < n; i++) X[i, i] = new ScalarValue(d[i], e[i]); return X; } } #endregion // Public Properties #region Public Methods /// /// Return the eigenvector matrix /// /// V public virtual MatrixValue GetV() { return new MatrixValue(V, n, n); } #endregion // Public Methods } }```

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

## Share

 Architect Germany
Florian lives in Munich, Germany. He started his programming career with Perl. After programming C/C++ for some years he discovered his favorite programming language C#. He did work at Siemens as a programmer until he decided to study Physics.

During his studies he worked as an IT consultant for various companies. After graduating with a PhD in theoretical particle Physics he is working as a senior technical consultant in the field of home automation and IoT.

Florian has been giving lectures in C#, HTML5 with CSS3 and JavaScript, software design, and other topics. He is regularly giving talks at user groups, conferences, and companies. He is actively contributing to open-source projects. Florian is the maintainer of AngleSharp, a completely managed browser engine.