Add php files

include/ECM/PHPExcel/Shared/JAMA/CholeskyDecomposition.php

@@ -0,0 +1,149 @@
+<?php
+/**
+ *	@package JAMA
+ *
+ *	Cholesky decomposition class
+ *
+ *	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.
+ *
+ *	@author Paul Meagher
+ *	@author Michael Bommarito
+ *	@version 1.2
+ */
+class CholeskyDecomposition {
+
+	/**
+	 *	Decomposition storage
+	 *	@var array
+	 *	@access private
+	 */
+	private $L = array();
+
+	/**
+	 *	Matrix row and column dimension
+	 *	@var int
+	 *	@access private
+	 */
+	private $m;
+
+	/**
+	 *	Symmetric positive definite flag
+	 *	@var boolean
+	 *	@access private
+	 */
+	private $isspd = true;
+
+
+	/**
+	 *	CholeskyDecomposition
+	 *
+	 *	Class constructor - decomposes symmetric positive definite matrix
+	 *	@param mixed Matrix square symmetric positive definite matrix
+	 */
+	public function __construct($A = null) {
+		if ($A instanceof Matrix) {
+			$this->L = $A->getArray();
+			$this->m = $A->getRowDimension();
+
+			for($i = 0; $i < $this->m; ++$i) {
+				for($j = $i; $j < $this->m; ++$j) {
+					for($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
+						$sum -= $this->L[$i][$k] * $this->L[$j][$k];
+					}
+					if ($i == $j) {
+						if ($sum >= 0) {
+							$this->L[$i][$i] = sqrt($sum);
+						} else {
+							$this->isspd = false;
+						}
+					} else {
+						if ($this->L[$i][$i] != 0) {
+							$this->L[$j][$i] = $sum / $this->L[$i][$i];
+						}
+					}
+				}
+
+				for ($k = $i+1; $k < $this->m; ++$k) {
+					$this->L[$i][$k] = 0.0;
+				}
+			}
+		} else {
+			throw new Exception(JAMAError(ArgumentTypeException));
+		}
+	}	//	function __construct()
+
+
+	/**
+	 *	Is the matrix symmetric and positive definite?
+	 *
+	 *	@return boolean
+	 */
+	public function isSPD() {
+		return $this->isspd;
+	}	//	function isSPD()
+
+
+	/**
+	 *	getL
+	 *
+	 *	Return triangular factor.
+	 *	@return Matrix Lower triangular matrix
+	 */
+	public function getL() {
+		return new Matrix($this->L);
+	}	//	function getL()
+
+
+	/**
+	 *	Solve A*X = B
+	 *
+	 *	@param $B Row-equal matrix
+	 *	@return Matrix L * L' * X = B
+	 */
+	public function solve($B = null) {
+		if ($B instanceof Matrix) {
+			if ($B->getRowDimension() == $this->m) {
+				if ($this->isspd) {
+					$X  = $B->getArrayCopy();
+					$nx = $B->getColumnDimension();
+
+					for ($k = 0; $k < $this->m; ++$k) {
+						for ($i = $k + 1; $i < $this->m; ++$i) {
+							for ($j = 0; $j < $nx; ++$j) {
+								$X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
+							}
+						}
+						for ($j = 0; $j < $nx; ++$j) {
+							$X[$k][$j] /= $this->L[$k][$k];
+						}
+					}
+
+					for ($k = $this->m - 1; $k >= 0; --$k) {
+						for ($j = 0; $j < $nx; ++$j) {
+							$X[$k][$j] /= $this->L[$k][$k];
+						}
+						for ($i = 0; $i < $k; ++$i) {
+							for ($j = 0; $j < $nx; ++$j) {
+								$X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
+							}
+						}
+					}
+
+					return new Matrix($X, $this->m, $nx);
+				} else {
+					throw new Exception(JAMAError(MatrixSPDException));
+				}
+			} else {
+				throw new Exception(JAMAError(MatrixDimensionException));
+			}
+		} else {
+			throw new Exception(JAMAError(ArgumentTypeException));
+		}
+	}	//	function solve()
+
+}	//	class CholeskyDecomposition

include/ECM/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php

@@ -0,0 +1,862 @@
+<?php
+/**
+ *	@package JAMA
+ *
+ *	Class to obtain 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.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().
+ *
+ *	@author  Paul Meagher
+ *	@license PHP v3.0
+ *	@version 1.1
+ */
+class EigenvalueDecomposition {
+
+	/**
+	 *	Row and column dimension (square matrix).
+	 *	@var int
+	 */
+	private $n;
+
+	/**
+	 *	Internal symmetry flag.
+	 *	@var int
+	 */
+	private $issymmetric;
+
+	/**
+	 *	Arrays for internal storage of eigenvalues.
+	 *	@var array
+	 */
+	private $d = array();
+	private $e = array();
+
+	/**
+	 *	Array for internal storage of eigenvectors.
+	 *	@var array
+	 */
+	private $V = array();
+
+	/**
+	*	Array for internal storage of nonsymmetric Hessenberg form.
+	*	@var array
+	*/
+	private $H = array();
+
+	/**
+	*	Working storage for nonsymmetric algorithm.
+	*	@var array
+	*/
+	private $ort;
+
+	/**
+	*	Used for complex scalar division.
+	*	@var float
+	*/
+	private $cdivr;
+	private $cdivi;
+
+
+	/**
+	 *	Symmetric Householder reduction to tridiagonal form.
+	 *
+	 *	@access private
+	 */
+	private function 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.
+		$this->d = $this->V[$this->n-1];
+		// Householder reduction to tridiagonal form.
+		for ($i = $this->n-1; $i > 0; --$i) {
+			$i_ = $i -1;
+			// Scale to avoid under/overflow.
+			$h = $scale = 0.0;
+			$scale += array_sum(array_map(abs, $this->d));
+			if ($scale == 0.0) {
+				$this->e[$i] = $this->d[$i_];
+				$this->d = array_slice($this->V[$i_], 0, $i_);
+				for ($j = 0; $j < $i; ++$j) {
+					$this->V[$j][$i] = $this->V[$i][$j] = 0.0;
+				}
+			} else {
+				// Generate Householder vector.
+				for ($k = 0; $k < $i; ++$k) {
+					$this->d[$k] /= $scale;
+					$h += pow($this->d[$k], 2);
+				}
+				$f = $this->d[$i_];
+				$g = sqrt($h);
+				if ($f > 0) {
+					$g = -$g;
+				}
+				$this->e[$i] = $scale * $g;
+				$h = $h - $f * $g;
+				$this->d[$i_] = $f - $g;
+				for ($j = 0; $j < $i; ++$j) {
+					$this->e[$j] = 0.0;
+				}
+				// Apply similarity transformation to remaining columns.
+				for ($j = 0; $j < $i; ++$j) {
+					$f = $this->d[$j];
+					$this->V[$j][$i] = $f;
+					$g = $this->e[$j] + $this->V[$j][$j] * $f;
+					for ($k = $j+1; $k <= $i_; ++$k) {
+						$g += $this->V[$k][$j] * $this->d[$k];
+						$this->e[$k] += $this->V[$k][$j] * $f;
+					}
+					$this->e[$j] = $g;
+				}
+				$f = 0.0;
+				for ($j = 0; $j < $i; ++$j) {
+					$this->e[$j] /= $h;
+					$f += $this->e[$j] * $this->d[$j];
+				}
+				$hh = $f / (2 * $h);
+				for ($j=0; $j < $i; ++$j) {
+					$this->e[$j] -= $hh * $this->d[$j];
+				}
+				for ($j = 0; $j < $i; ++$j) {
+					$f = $this->d[$j];
+					$g = $this->e[$j];
+					for ($k = $j; $k <= $i_; ++$k) {
+						$this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
+					}
+					$this->d[$j] = $this->V[$i-1][$j];
+					$this->V[$i][$j] = 0.0;
+				}
+			}
+			$this->d[$i] = $h;
+		}
+
+		// Accumulate transformations.
+		for ($i = 0; $i < $this->n-1; ++$i) {
+			$this->V[$this->n-1][$i] = $this->V[$i][$i];
+			$this->V[$i][$i] = 1.0;
+			$h = $this->d[$i+1];
+			if ($h != 0.0) {
+				for ($k = 0; $k <= $i; ++$k) {
+					$this->d[$k] = $this->V[$k][$i+1] / $h;
+				}
+				for ($j = 0; $j <= $i; ++$j) {
+					$g = 0.0;
+					for ($k = 0; $k <= $i; ++$k) {
+						$g += $this->V[$k][$i+1] * $this->V[$k][$j];
+					}
+					for ($k = 0; $k <= $i; ++$k) {
+						$this->V[$k][$j] -= $g * $this->d[$k];
+					}
+				}
+			}
+			for ($k = 0; $k <= $i; ++$k) {
+				$this->V[$k][$i+1] = 0.0;
+			}
+		}
+
+		$this->d = $this->V[$this->n-1];
+		$this->V[$this->n-1] = array_fill(0, $j, 0.0);
+		$this->V[$this->n-1][$this->n-1] = 1.0;
+		$this->e[0] = 0.0;
+	}
+
+
+	/**
+	 *	Symmetric tridiagonal QL algorithm.
+	 *
+	 *	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.
+	 *
+	 *	@access private
+	 */
+	private function tql2() {
+		for ($i = 1; $i < $this->n; ++$i) {
+			$this->e[$i-1] = $this->e[$i];
+		}
+		$this->e[$this->n-1] = 0.0;
+		$f = 0.0;
+		$tst1 = 0.0;
+		$eps  = pow(2.0,-52.0);
+
+		for ($l = 0; $l < $this->n; ++$l) {
+			// Find small subdiagonal element
+			$tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
+			$m = $l;
+			while ($m < $this->n) {
+				if (abs($this->e[$m]) <= $eps * $tst1)
+					break;
+				++$m;
+			}
+			// If m == l, $this->d[l] is an eigenvalue,
+			// otherwise, iterate.
+			if ($m > $l) {
+				$iter = 0;
+				do {
+					// Could check iteration count here.
+					$iter += 1;
+					// Compute implicit shift
+					$g = $this->d[$l];
+					$p = ($this->d[$l+1] - $g) / (2.0 * $this->e[$l]);
+					$r = hypo($p, 1.0);
+					if ($p < 0)
+						$r *= -1;
+					$this->d[$l] = $this->e[$l] / ($p + $r);
+					$this->d[$l+1] = $this->e[$l] * ($p + $r);
+					$dl1 = $this->d[$l+1];
+					$h = $g - $this->d[$l];
+					for ($i = $l + 2; $i < $this->n; ++$i)
+						$this->d[$i] -= $h;
+					$f += $h;
+					// Implicit QL transformation.
+					$p = $this->d[$m];
+					$c = 1.0;
+					$c2 = $c3 = $c;
+					$el1 = $this->e[$l + 1];
+					$s = $s2 = 0.0;
+					for ($i = $m-1; $i >= $l; --$i) {
+						$c3 = $c2;
+						$c2 = $c;
+						$s2 = $s;
+						$g  = $c * $this->e[$i];
+						$h  = $c * $p;
+						$r  = hypo($p, $this->e[$i]);
+						$this->e[$i+1] = $s * $r;
+						$s = $this->e[$i] / $r;
+						$c = $p / $r;
+						$p = $c * $this->d[$i] - $s * $g;
+						$this->d[$i+1] = $h + $s * ($c * $g + $s * $this->d[$i]);
+						// Accumulate transformation.
+						for ($k = 0; $k < $this->n; ++$k) {
+							$h = $this->V[$k][$i+1];
+							$this->V[$k][$i+1] = $s * $this->V[$k][$i] + $c * $h;
+							$this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
+						}
+					}
+					$p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
+					$this->e[$l] = $s * $p;
+					$this->d[$l] = $c * $p;
+				// Check for convergence.
+				} while (abs($this->e[$l]) > $eps * $tst1);
+			}
+			$this->d[$l] = $this->d[$l] + $f;
+			$this->e[$l] = 0.0;
+		}
+
+		// Sort eigenvalues and corresponding vectors.
+		for ($i = 0; $i < $this->n - 1; ++$i) {
+			$k = $i;
+			$p = $this->d[$i];
+			for ($j = $i+1; $j < $this->n; ++$j) {
+				if ($this->d[$j] < $p) {
+					$k = $j;
+					$p = $this->d[$j];
+				}
+			}
+			if ($k != $i) {
+				$this->d[$k] = $this->d[$i];
+				$this->d[$i] = $p;
+				for ($j = 0; $j < $this->n; ++$j) {
+					$p = $this->V[$j][$i];
+					$this->V[$j][$i] = $this->V[$j][$k];
+					$this->V[$j][$k] = $p;
+				}
+			}
+		}
+	}
+
+
+	/**
+	 *	Nonsymmetric reduction to Hessenberg form.
+	 *
+	 *	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.
+	 *
+	 *	@access private
+	 */
+	private function orthes () {
+		$low  = 0;
+		$high = $this->n-1;
+
+		for ($m = $low+1; $m <= $high-1; ++$m) {
+			// Scale column.
+			$scale = 0.0;
+			for ($i = $m; $i <= $high; ++$i) {
+				$scale = $scale + abs($this->H[$i][$m-1]);
+			}
+			if ($scale != 0.0) {
+				// Compute Householder transformation.
+				$h = 0.0;
+				for ($i = $high; $i >= $m; --$i) {
+					$this->ort[$i] = $this->H[$i][$m-1] / $scale;
+					$h += $this->ort[$i] * $this->ort[$i];
+				}
+				$g = sqrt($h);
+				if ($this->ort[$m] > 0) {
+					$g *= -1;
+				}
+				$h -= $this->ort[$m] * $g;
+				$this->ort[$m] -= $g;
+				// Apply Householder similarity transformation
+				// H = (I -u * u' / h) * H * (I -u * u') / h)
+				for ($j = $m; $j < $this->n; ++$j) {
+					$f = 0.0;
+					for ($i = $high; $i >= $m; --$i) {
+						$f += $this->ort[$i] * $this->H[$i][$j];
+					}
+					$f /= $h;
+					for ($i = $m; $i <= $high; ++$i) {
+						$this->H[$i][$j] -= $f * $this->ort[$i];
+					}
+				}
+				for ($i = 0; $i <= $high; ++$i) {
+					$f = 0.0;
+					for ($j = $high; $j >= $m; --$j) {
+						$f += $this->ort[$j] * $this->H[$i][$j];
+					}
+					$f = $f / $h;
+					for ($j = $m; $j <= $high; ++$j) {
+						$this->H[$i][$j] -= $f * $this->ort[$j];
+					}
+				}
+				$this->ort[$m] = $scale * $this->ort[$m];
+				$this->H[$m][$m-1] = $scale * $g;
+			}
+		}
+
+		// Accumulate transformations (Algol's ortran).
+		for ($i = 0; $i < $this->n; ++$i) {
+			for ($j = 0; $j < $this->n; ++$j) {
+				$this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
+			}
+		}
+		for ($m = $high-1; $m >= $low+1; --$m) {
+			if ($this->H[$m][$m-1] != 0.0) {
+				for ($i = $m+1; $i <= $high; ++$i) {
+					$this->ort[$i] = $this->H[$i][$m-1];
+				}
+				for ($j = $m; $j <= $high; ++$j) {
+					$g = 0.0;
+					for ($i = $m; $i <= $high; ++$i) {
+						$g += $this->ort[$i] * $this->V[$i][$j];
+					}
+					// Double division avoids possible underflow
+					$g = ($g / $this->ort[$m]) / $this->H[$m][$m-1];
+					for ($i = $m; $i <= $high; ++$i) {
+						$this->V[$i][$j] += $g * $this->ort[$i];
+					}
+				}
+			}
+		}
+	}
+
+
+	/**
+	 *	Performs complex division.
+	 *
+	 *	@access private
+	 */
+	private function cdiv($xr, $xi, $yr, $yi) {
+		if (abs($yr) > abs($yi)) {
+			$r = $yi / $yr;
+			$d = $yr + $r * $yi;
+			$this->cdivr = ($xr + $r * $xi) / $d;
+			$this->cdivi = ($xi - $r * $xr) / $d;
+		} else {
+			$r = $yr / $yi;
+			$d = $yi + $r * $yr;
+			$this->cdivr = ($r * $xr + $xi) / $d;
+			$this->cdivi = ($r * $xi - $xr) / $d;
+		}
+	}
+
+
+	/**
+	 *	Nonsymmetric reduction from Hessenberg to real Schur form.
+	 *
+	 *	Code 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.
+	 *
+	 *	@access private
+	 */
+	private function hqr2 () {
+		//  Initialize
+		$nn = $this->n;
+		$n  = $nn - 1;
+		$low = 0;
+		$high = $nn - 1;
+		$eps = pow(2.0, -52.0);
+		$exshift = 0.0;
+		$p = $q = $r = $s = $z = 0;
+		// Store roots isolated by balanc and compute matrix norm
+		$norm = 0.0;
+
+		for ($i = 0; $i < $nn; ++$i) {
+			if (($i < $low) OR ($i > $high)) {
+				$this->d[$i] = $this->H[$i][$i];
+				$this->e[$i] = 0.0;
+			}
+			for ($j = max($i-1, 0); $j < $nn; ++$j) {
+				$norm = $norm + abs($this->H[$i][$j]);
+			}
+		}
+
+		// Outer loop over eigenvalue index
+		$iter = 0;
+		while ($n >= $low) {
+			// Look for single small sub-diagonal element
+			$l = $n;
+			while ($l > $low) {
+				$s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]);
+				if ($s == 0.0) {
+					$s = $norm;
+				}
+				if (abs($this->H[$l][$l-1]) < $eps * $s) {
+					break;
+				}
+				--$l;
+			}
+			// Check for convergence
+			// One root found
+			if ($l == $n) {
+				$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
+				$this->d[$n] = $this->H[$n][$n];
+				$this->e[$n] = 0.0;
+				--$n;
+				$iter = 0;
+			// Two roots found
+			} else if ($l == $n-1) {
+				$w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
+				$p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0;
+				$q = $p * $p + $w;
+				$z = sqrt(abs($q));
+				$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
+				$this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift;
+				$x = $this->H[$n][$n];
+				// Real pair
+				if ($q >= 0) {
+					if ($p >= 0) {
+						$z = $p + $z;
+					} else {
+						$z = $p - $z;
+					}
+					$this->d[$n-1] = $x + $z;
+					$this->d[$n] = $this->d[$n-1];
+					if ($z != 0.0) {
+						$this->d[$n] = $x - $w / $z;
+					}
+					$this->e[$n-1] = 0.0;
+					$this->e[$n] = 0.0;
+					$x = $this->H[$n][$n-1];
+					$s = abs($x) + abs($z);
+					$p = $x / $s;
+					$q = $z / $s;
+					$r = sqrt($p * $p + $q * $q);
+					$p = $p / $r;
+					$q = $q / $r;
+					// Row modification
+					for ($j = $n-1; $j < $nn; ++$j) {
+						$z = $this->H[$n-1][$j];
+						$this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j];
+						$this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
+					}
+					// Column modification
+					for ($i = 0; $i <= n; ++$i) {
+						$z = $this->H[$i][$n-1];
+						$this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n];
+						$this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
+					}
+					// Accumulate transformations
+					for ($i = $low; $i <= $high; ++$i) {
+						$z = $this->V[$i][$n-1];
+						$this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n];
+						$this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
+					}
+				// Complex pair
+				} else {
+					$this->d[$n-1] = $x + $p;
+					$this->d[$n] = $x + $p;
+					$this->e[$n-1] = $z;
+					$this->e[$n] = -$z;
+				}
+				$n = $n - 2;
+				$iter = 0;
+			// No convergence yet
+			} else {
+				// Form shift
+				$x = $this->H[$n][$n];
+				$y = 0.0;
+				$w = 0.0;
+				if ($l < $n) {
+					$y = $this->H[$n-1][$n-1];
+					$w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
+				}
+				// Wilkinson's original ad hoc shift
+				if ($iter == 10) {
+					$exshift += $x;
+					for ($i = $low; $i <= $n; ++$i) {
+						$this->H[$i][$i] -= $x;
+					}
+					$s = abs($this->H[$n][$n-1]) + abs($this->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 = sqrt($s);
+						if ($y < $x) {
+							$s = -$s;
+						}
+						$s = $x - $w / (($y - $x) / 2.0 + $s);
+						for ($i = $low; $i <= $n; ++$i) {
+							$this->H[$i][$i] -= $s;
+						}
+						$exshift += $s;
+						$x = $y = $w = 0.964;
+					}
+				}
+				// Could check iteration count here.
+				$iter = $iter + 1;
+				// Look for two consecutive small sub-diagonal elements
+				$m = $n - 2;
+				while ($m >= $l) {
+					$z = $this->H[$m][$m];
+					$r = $x - $z;
+					$s = $y - $z;
+					$p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1];
+					$q = $this->H[$m+1][$m+1] - $z - $r - $s;
+					$r = $this->H[$m+2][$m+1];
+					$s = abs($p) + abs($q) + abs($r);
+					$p = $p / $s;
+					$q = $q / $s;
+					$r = $r / $s;
+					if ($m == $l) {
+						break;
+					}
+					if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) <
+						$eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) {
+						break;
+					}
+					--$m;
+				}
+				for ($i = $m + 2; $i <= $n; ++$i) {
+					$this->H[$i][$i-2] = 0.0;
+					if ($i > $m+2) {
+						$this->H[$i][$i-3] = 0.0;
+					}
+				}
+				// Double QR step involving rows l:n and columns m:n
+				for ($k = $m; $k <= $n-1; ++$k) {
+					$notlast = ($k != $n-1);
+					if ($k != $m) {
+						$p = $this->H[$k][$k-1];
+						$q = $this->H[$k+1][$k-1];
+						$r = ($notlast ? $this->H[$k+2][$k-1] : 0.0);
+						$x = abs($p) + abs($q) + abs($r);
+						if ($x != 0.0) {
+							$p = $p / $x;
+							$q = $q / $x;
+							$r = $r / $x;
+						}
+					}
+					if ($x == 0.0) {
+						break;
+					}
+					$s = sqrt($p * $p + $q * $q + $r * $r);
+					if ($p < 0) {
+						$s = -$s;
+					}
+					if ($s != 0) {
+						if ($k != $m) {
+							$this->H[$k][$k-1] = -$s * $x;
+						} elseif ($l != $m) {
+							$this->H[$k][$k-1] = -$this->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 ($j = $k; $j < $nn; ++$j) {
+							$p = $this->H[$k][$j] + $q * $this->H[$k+1][$j];
+							if ($notlast) {
+								$p = $p + $r * $this->H[$k+2][$j];
+								$this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z;
+							}
+							$this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
+							$this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y;
+						}
+						// Column modification
+						for ($i = 0; $i <= min($n, $k+3); ++$i) {
+							$p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1];
+							if ($notlast) {
+								$p = $p + $z * $this->H[$i][$k+2];
+								$this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r;
+							}
+							$this->H[$i][$k] = $this->H[$i][$k] - $p;
+							$this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q;
+						}
+						// Accumulate transformations
+						for ($i = $low; $i <= $high; ++$i) {
+							$p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1];
+							if ($notlast) {
+								$p = $p + $z * $this->V[$i][$k+2];
+								$this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r;
+							}
+							$this->V[$i][$k] = $this->V[$i][$k] - $p;
+							$this->V[$i][$k+1] = $this->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 = $this->d[$n];
+			$q = $this->e[$n];
+			// Real vector
+			if ($q == 0) {
+				$l = $n;
+				$this->H[$n][$n] = 1.0;
+				for ($i = $n-1; $i >= 0; --$i) {
+					$w = $this->H[$i][$i] - $p;
+					$r = 0.0;
+					for ($j = $l; $j <= $n; ++$j) {
+						$r = $r + $this->H[$i][$j] * $this->H[$j][$n];
+					}
+					if ($this->e[$i] < 0.0) {
+						$z = $w;
+						$s = $r;
+					} else {
+						$l = $i;
+						if ($this->e[$i] == 0.0) {
+							if ($w != 0.0) {
+								$this->H[$i][$n] = -$r / $w;
+							} else {
+								$this->H[$i][$n] = -$r / ($eps * $norm);
+							}
+						// Solve real equations
+						} else {
+							$x = $this->H[$i][$i+1];
+							$y = $this->H[$i+1][$i];
+							$q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
+							$t = ($x * $s - $z * $r) / $q;
+							$this->H[$i][$n] = $t;
+							if (abs($x) > abs($z)) {
+								$this->H[$i+1][$n] = (-$r - $w * $t) / $x;
+							} else {
+								$this->H[$i+1][$n] = (-$s - $y * $t) / $z;
+							}
+						}
+						// Overflow control
+						$t = abs($this->H[$i][$n]);
+						if (($eps * $t) * $t > 1) {
+							for ($j = $i; $j <= $n; ++$j) {
+								$this->H[$j][$n] = $this->H[$j][$n] / $t;
+							}
+						}
+					}
+				}
+			// Complex vector
+			} else if ($q < 0) {
+				$l = $n-1;
+				// Last vector component imaginary so matrix is triangular
+				if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) {
+					$this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1];
+					$this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1];
+				} else {
+					$this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q);
+					$this->H[$n-1][$n-1] = $this->cdivr;
+					$this->H[$n-1][$n]   = $this->cdivi;
+				}
+				$this->H[$n][$n-1] = 0.0;
+				$this->H[$n][$n] = 1.0;
+				for ($i = $n-2; $i >= 0; --$i) {
+					// double ra,sa,vr,vi;
+					$ra = 0.0;
+					$sa = 0.0;
+					for ($j = $l; $j <= $n; ++$j) {
+						$ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1];
+						$sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
+					}
+					$w = $this->H[$i][$i] - $p;
+					if ($this->e[$i] < 0.0) {
+						$z = $w;
+						$r = $ra;
+						$s = $sa;
+					} else {
+						$l = $i;
+						if ($this->e[$i] == 0) {
+							$this->cdiv(-$ra, -$sa, $w, $q);
+							$this->H[$i][$n-1] = $this->cdivr;
+							$this->H[$i][$n]   = $this->cdivi;
+						} else {
+							// Solve complex equations
+							$x = $this->H[$i][$i+1];
+							$y = $this->H[$i+1][$i];
+							$vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
+							$vi = ($this->d[$i] - $p) * 2.0 * $q;
+							if ($vr == 0.0 & $vi == 0.0) {
+								$vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
+							}
+							$this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
+							$this->H[$i][$n-1] = $this->cdivr;
+							$this->H[$i][$n]   = $this->cdivi;
+							if (abs($x) > (abs($z) + abs($q))) {
+								$this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x;
+								$this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x;
+							} else {
+								$this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q);
+								$this->H[$i+1][$n-1] = $this->cdivr;
+								$this->H[$i+1][$n]   = $this->cdivi;
+							}
+						}
+						// Overflow control
+						$t = max(abs($this->H[$i][$n-1]),abs($this->H[$i][$n]));
+						if (($eps * $t) * $t > 1) {
+							for ($j = $i; $j <= $n; ++$j) {
+								$this->H[$j][$n-1] = $this->H[$j][$n-1] / $t;
+								$this->H[$j][$n]   = $this->H[$j][$n] / $t;
+							}
+						}
+					} // end else
+				} // end for
+			} // end else for complex case
+		} // end for
+
+		// Vectors of isolated roots
+		for ($i = 0; $i < $nn; ++$i) {
+			if ($i < $low | $i > $high) {
+				for ($j = $i; $j < $nn; ++$j) {
+					$this->V[$i][$j] = $this->H[$i][$j];
+				}
+			}
+		}
+
+		// Back transformation to get eigenvectors of original matrix
+		for ($j = $nn-1; $j >= $low; --$j) {
+			for ($i = $low; $i <= $high; ++$i) {
+				$z = 0.0;
+				for ($k = $low; $k <= min($j,$high); ++$k) {
+					$z = $z + $this->V[$i][$k] * $this->H[$k][$j];
+				}
+				$this->V[$i][$j] = $z;
+			}
+		}
+	} // end hqr2
+
+
+	/**
+	 *	Constructor: Check for symmetry, then construct the eigenvalue decomposition
+	 *
+	 *	@access public
+	 *	@param A  Square matrix
+	 *	@return Structure to access D and V.
+	 */
+	public function __construct($Arg) {
+		$this->A = $Arg->getArray();
+		$this->n = $Arg->getColumnDimension();
+
+		$issymmetric = true;
+		for ($j = 0; ($j < $this->n) & $issymmetric; ++$j) {
+			for ($i = 0; ($i < $this->n) & $issymmetric; ++$i) {
+				$issymmetric = ($this->A[$i][$j] == $this->A[$j][$i]);
+			}
+		}
+
+		if ($issymmetric) {
+			$this->V = $this->A;
+			// Tridiagonalize.
+			$this->tred2();
+			// Diagonalize.
+			$this->tql2();
+		} else {
+			$this->H = $this->A;
+			$this->ort = array();
+			// Reduce to Hessenberg form.
+			$this->orthes();
+			// Reduce Hessenberg to real Schur form.
+			$this->hqr2();
+		}
+	}
+
+
+	/**
+	 *	Return the eigenvector matrix
+	 *
+	 *	@access public
+	 *	@return V
+	 */
+	public function getV() {
+		return new Matrix($this->V, $this->n, $this->n);
+	}
+
+
+	/**
+	 *	Return the real parts of the eigenvalues
+	 *
+	 *	@access public
+	 *	@return real(diag(D))
+	 */
+	public function getRealEigenvalues() {
+		return $this->d;
+	}
+
+
+	/**
+	 *	Return the imaginary parts of the eigenvalues
+	 *
+	 *	@access public
+	 *	@return imag(diag(D))
+	 */
+	public function getImagEigenvalues() {
+		return $this->e;
+	}
+
+
+	/**
+	 *	Return the block diagonal eigenvalue matrix
+	 *
+	 *	@access public
+	 *	@return D
+	 */
+	public function getD() {
+		for ($i = 0; $i < $this->n; ++$i) {
+			$D[$i] = array_fill(0, $this->n, 0.0);
+			$D[$i][$i] = $this->d[$i];
+			if ($this->e[$i] == 0) {
+				continue;
+			}
+			$o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
+			$D[$i][$o] = $this->e[$i];
+		}
+		return new Matrix($D);
+	}
+
+}	//	class EigenvalueDecomposition

include/ECM/PHPExcel/Shared/JAMA/LUDecomposition.php

@@ -0,0 +1,258 @@
+<?php
+/**
+ *	@package JAMA
+ *
+ *	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.
+ *
+ *	@author Paul Meagher
+ *	@author Bartosz Matosiuk
+ *	@author Michael Bommarito
+ *	@version 1.1
+ *	@license PHP v3.0
+ */
+class PHPExcel_Shared_JAMA_LUDecomposition {
+
+	const MatrixSingularException	= "Can only perform operation on singular matrix.";
+	const MatrixSquareException		= "Mismatched Row dimension";
+
+	/**
+	 *	Decomposition storage
+	 *	@var array
+	 */
+	private $LU = array();
+
+	/**
+	 *	Row dimension.
+	 *	@var int
+	 */
+	private $m;
+
+	/**
+	 *	Column dimension.
+	 *	@var int
+	 */
+	private $n;
+
+	/**
+	 *	Pivot sign.
+	 *	@var int
+	 */
+	private $pivsign;
+
+	/**
+	 *	Internal storage of pivot vector.
+	 *	@var array
+	 */
+	private $piv = array();
+
+
+	/**
+	 *	LU Decomposition constructor.
+	 *
+	 *	@param $A Rectangular matrix
+	 *	@return Structure to access L, U and piv.
+	 */
+	public function __construct($A) {
+		if ($A instanceof PHPExcel_Shared_JAMA_Matrix) {
+			// Use a "left-looking", dot-product, Crout/Doolittle algorithm.
+			$this->LU = $A->getArray();
+			$this->m  = $A->getRowDimension();
+			$this->n  = $A->getColumnDimension();
+			for ($i = 0; $i < $this->m; ++$i) {
+				$this->piv[$i] = $i;
+			}
+			$this->pivsign = 1;
+			$LUrowi = $LUcolj = array();
+
+			// Outer loop.
+			for ($j = 0; $j < $this->n; ++$j) {
+				// Make a copy of the j-th column to localize references.
+				for ($i = 0; $i < $this->m; ++$i) {
+					$LUcolj[$i] = &$this->LU[$i][$j];
+				}
+				// Apply previous transformations.
+				for ($i = 0; $i < $this->m; ++$i) {
+					$LUrowi = $this->LU[$i];
+					// Most of the time is spent in the following dot product.
+					$kmax = min($i,$j);
+					$s = 0.0;
+					for ($k = 0; $k < $kmax; ++$k) {
+						$s += $LUrowi[$k] * $LUcolj[$k];
+					}
+					$LUrowi[$j] = $LUcolj[$i] -= $s;
+				}
+				// Find pivot and exchange if necessary.
+				$p = $j;
+				for ($i = $j+1; $i < $this->m; ++$i) {
+					if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
+						$p = $i;
+					}
+				}
+				if ($p != $j) {
+					for ($k = 0; $k < $this->n; ++$k) {
+						$t = $this->LU[$p][$k];
+						$this->LU[$p][$k] = $this->LU[$j][$k];
+						$this->LU[$j][$k] = $t;
+					}
+					$k = $this->piv[$p];
+					$this->piv[$p] = $this->piv[$j];
+					$this->piv[$j] = $k;
+					$this->pivsign = $this->pivsign * -1;
+				}
+				// Compute multipliers.
+				if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
+					for ($i = $j+1; $i < $this->m; ++$i) {
+						$this->LU[$i][$j] /= $this->LU[$j][$j];
+					}
+				}
+			}
+		} else {
+			throw new Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
+		}
+	}	//	function __construct()
+
+
+	/**
+	 *	Get lower triangular factor.
+	 *
+	 *	@return array Lower triangular factor
+	 */
+	public function getL() {
+		for ($i = 0; $i < $this->m; ++$i) {
+			for ($j = 0; $j < $this->n; ++$j) {
+				if ($i > $j) {
+					$L[$i][$j] = $this->LU[$i][$j];
+				} elseif ($i == $j) {
+					$L[$i][$j] = 1.0;
+				} else {
+					$L[$i][$j] = 0.0;
+				}
+			}
+		}
+		return new PHPExcel_Shared_JAMA_Matrix($L);
+	}	//	function getL()
+
+
+	/**
+	 *	Get upper triangular factor.
+	 *
+	 *	@return array Upper triangular factor
+	 */
+	public function getU() {
+		for ($i = 0; $i < $this->n; ++$i) {
+			for ($j = 0; $j < $this->n; ++$j) {
+				if ($i <= $j) {
+					$U[$i][$j] = $this->LU[$i][$j];
+				} else {
+					$U[$i][$j] = 0.0;
+				}
+			}
+		}
+		return new PHPExcel_Shared_JAMA_Matrix($U);
+	}	//	function getU()
+
+
+	/**
+	 *	Return pivot permutation vector.
+	 *
+	 *	@return array Pivot vector
+	 */
+	public function getPivot() {
+		return $this->piv;
+	}	//	function getPivot()
+
+
+	/**
+	 *	Alias for getPivot
+	 *
+	 *	@see getPivot
+	 */
+	public function getDoublePivot() {
+		return $this->getPivot();
+	}	//	function getDoublePivot()
+
+
+	/**
+	 *	Is the matrix nonsingular?
+	 *
+	 *	@return true if U, and hence A, is nonsingular.
+	 */
+	public function isNonsingular() {
+		for ($j = 0; $j < $this->n; ++$j) {
+			if ($this->LU[$j][$j] == 0) {
+				return false;
+			}
+		}
+		return true;
+	}	//	function isNonsingular()
+
+
+	/**
+	 *	Count determinants
+	 *
+	 *	@return array d matrix deterninat
+	 */
+	public function det() {
+		if ($this->m == $this->n) {
+			$d = $this->pivsign;
+			for ($j = 0; $j < $this->n; ++$j) {
+				$d *= $this->LU[$j][$j];
+			}
+			return $d;
+		} else {
+			throw new Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
+		}
+	}	//	function det()
+
+
+	/**
+	 *	Solve A*X = B
+	 *
+	 *	@param  $B  A Matrix with as many rows as A and any number of columns.
+	 *	@return  X so that L*U*X = B(piv,:)
+	 *	@exception  IllegalArgumentException Matrix row dimensions must agree.
+	 *	@exception  RuntimeException  Matrix is singular.
+	 */
+	public function solve($B) {
+		if ($B->getRowDimension() == $this->m) {
+			if ($this->isNonsingular()) {
+				// Copy right hand side with pivoting
+				$nx = $B->getColumnDimension();
+				$X  = $B->getMatrix($this->piv, 0, $nx-1);
+				// Solve L*Y = B(piv,:)
+				for ($k = 0; $k < $this->n; ++$k) {
+					for ($i = $k+1; $i < $this->n; ++$i) {
+						for ($j = 0; $j < $nx; ++$j) {
+							$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
+						}
+					}
+				}
+				// Solve U*X = Y;
+				for ($k = $this->n-1; $k >= 0; --$k) {
+					for ($j = 0; $j < $nx; ++$j) {
+						$X->A[$k][$j] /= $this->LU[$k][$k];
+					}
+					for ($i = 0; $i < $k; ++$i) {
+						for ($j = 0; $j < $nx; ++$j) {
+							$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
+						}
+					}
+				}
+				return $X;
+			} else {
+				throw new Exception(self::MatrixSingularException);
+			}
+		} else {
+			throw new Exception(self::MatrixSquareException);
+		}
+	}	//	function solve()
+
+}	//	class PHPExcel_Shared_JAMA_LUDecomposition

include/ECM/PHPExcel/Shared/JAMA/QRDecomposition.php

@@ -0,0 +1,234 @@
+<?php
+/**
+ *	@package JAMA
+ *
+ *	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.
+ *
+ *	@author  Paul Meagher
+ *	@license PHP v3.0
+ *	@version 1.1
+ */
+class PHPExcel_Shared_JAMA_QRDecomposition {
+
+	const MatrixRankException	= "Can only perform operation on full-rank matrix.";
+
+	/**
+	 *	Array for internal storage of decomposition.
+	 *	@var array
+	 */
+	private $QR = array();
+
+	/**
+	 *	Row dimension.
+	 *	@var integer
+	 */
+	private $m;
+
+	/**
+	*	Column dimension.
+	*	@var integer
+	*/
+	private $n;
+
+	/**
+	 *	Array for internal storage of diagonal of R.
+	 *	@var  array
+	 */
+	private $Rdiag = array();
+
+
+	/**
+	 *	QR Decomposition computed by Householder reflections.
+	 *
+	 *	@param matrix $A Rectangular matrix
+	 *	@return Structure to access R and the Householder vectors and compute Q.
+	 */
+	public function __construct($A) {
+		if($A instanceof PHPExcel_Shared_JAMA_Matrix) {
+			// Initialize.
+			$this->QR = $A->getArrayCopy();
+			$this->m  = $A->getRowDimension();
+			$this->n  = $A->getColumnDimension();
+			// Main loop.
+			for ($k = 0; $k < $this->n; ++$k) {
+				// Compute 2-norm of k-th column without under/overflow.
+				$nrm = 0.0;
+				for ($i = $k; $i < $this->m; ++$i) {
+					$nrm = hypo($nrm, $this->QR[$i][$k]);
+				}
+				if ($nrm != 0.0) {
+					// Form k-th Householder vector.
+					if ($this->QR[$k][$k] < 0) {
+						$nrm = -$nrm;
+					}
+					for ($i = $k; $i < $this->m; ++$i) {
+						$this->QR[$i][$k] /= $nrm;
+					}
+					$this->QR[$k][$k] += 1.0;
+					// Apply transformation to remaining columns.
+					for ($j = $k+1; $j < $this->n; ++$j) {
+						$s = 0.0;
+						for ($i = $k; $i < $this->m; ++$i) {
+							$s += $this->QR[$i][$k] * $this->QR[$i][$j];
+						}
+						$s = -$s/$this->QR[$k][$k];
+						for ($i = $k; $i < $this->m; ++$i) {
+							$this->QR[$i][$j] += $s * $this->QR[$i][$k];
+						}
+					}
+				}
+				$this->Rdiag[$k] = -$nrm;
+			}
+		} else {
+			throw new Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
+		}
+	}	//	function __construct()
+
+
+	/**
+	 *	Is the matrix full rank?
+	 *
+	 *	@return boolean true if R, and hence A, has full rank, else false.
+	 */
+	public function isFullRank() {
+		for ($j = 0; $j < $this->n; ++$j) {
+			if ($this->Rdiag[$j] == 0) {
+				return false;
+			}
+		}
+		return true;
+	}	//	function isFullRank()
+
+
+	/**
+	 *	Return the Householder vectors
+	 *
+	 *	@return Matrix Lower trapezoidal matrix whose columns define the reflections
+	 */
+	public function getH() {
+		for ($i = 0; $i < $this->m; ++$i) {
+			for ($j = 0; $j < $this->n; ++$j) {
+				if ($i >= $j) {
+					$H[$i][$j] = $this->QR[$i][$j];
+				} else {
+					$H[$i][$j] = 0.0;
+				}
+			}
+		}
+		return new PHPExcel_Shared_JAMA_Matrix($H);
+	}	//	function getH()
+
+
+	/**
+	 *	Return the upper triangular factor
+	 *
+	 *	@return Matrix upper triangular factor
+	 */
+	public function getR() {
+		for ($i = 0; $i < $this->n; ++$i) {
+			for ($j = 0; $j < $this->n; ++$j) {
+				if ($i < $j) {
+					$R[$i][$j] = $this->QR[$i][$j];
+				} elseif ($i == $j) {
+					$R[$i][$j] = $this->Rdiag[$i];
+				} else {
+					$R[$i][$j] = 0.0;
+				}
+			}
+		}
+		return new PHPExcel_Shared_JAMA_Matrix($R);
+	}	//	function getR()
+
+
+	/**
+	 *	Generate and return the (economy-sized) orthogonal factor
+	 *
+	 *	@return Matrix orthogonal factor
+	 */
+	public function getQ() {
+		for ($k = $this->n-1; $k >= 0; --$k) {
+			for ($i = 0; $i < $this->m; ++$i) {
+				$Q[$i][$k] = 0.0;
+			}
+			$Q[$k][$k] = 1.0;
+			for ($j = $k; $j < $this->n; ++$j) {
+				if ($this->QR[$k][$k] != 0) {
+					$s = 0.0;
+					for ($i = $k; $i < $this->m; ++$i) {
+						$s += $this->QR[$i][$k] * $Q[$i][$j];
+					}
+					$s = -$s/$this->QR[$k][$k];
+					for ($i = $k; $i < $this->m; ++$i) {
+						$Q[$i][$j] += $s * $this->QR[$i][$k];
+					}
+				}
+			}
+		}
+		/*
+		for($i = 0; $i < count($Q); ++$i) {
+			for($j = 0; $j < count($Q); ++$j) {
+				if(! isset($Q[$i][$j]) ) {
+					$Q[$i][$j] = 0;
+				}
+			}
+		}
+		*/
+		return new PHPExcel_Shared_JAMA_Matrix($Q);
+	}	//	function getQ()
+
+
+	/**
+	 *	Least squares solution of A*X = B
+	 *
+	 *	@param Matrix $B A Matrix with as many rows as A and any number of columns.
+	 *	@return Matrix Matrix that minimizes the two norm of Q*R*X-B.
+	 */
+	public function solve($B) {
+		if ($B->getRowDimension() == $this->m) {
+			if ($this->isFullRank()) {
+				// Copy right hand side
+				$nx = $B->getColumnDimension();
+				$X  = $B->getArrayCopy();
+				// Compute Y = transpose(Q)*B
+				for ($k = 0; $k < $this->n; ++$k) {
+					for ($j = 0; $j < $nx; ++$j) {
+						$s = 0.0;
+						for ($i = $k; $i < $this->m; ++$i) {
+							$s += $this->QR[$i][$k] * $X[$i][$j];
+						}
+						$s = -$s/$this->QR[$k][$k];
+						for ($i = $k; $i < $this->m; ++$i) {
+							$X[$i][$j] += $s * $this->QR[$i][$k];
+						}
+					}
+				}
+				// Solve R*X = Y;
+				for ($k = $this->n-1; $k >= 0; --$k) {
+					for ($j = 0; $j < $nx; ++$j) {
+						$X[$k][$j] /= $this->Rdiag[$k];
+					}
+					for ($i = 0; $i < $k; ++$i) {
+						for ($j = 0; $j < $nx; ++$j) {
+							$X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];
+						}
+					}
+				}
+				$X = new PHPExcel_Shared_JAMA_Matrix($X);
+				return ($X->getMatrix(0, $this->n-1, 0, $nx));
+			} else {
+				throw new Exception(self::MatrixRankException);
+			}
+		} else {
+			throw new Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
+		}
+	}	//	function solve()
+
+}	//	PHPExcel_Shared_JAMA_class QRDecomposition

include/ECM/PHPExcel/Shared/JAMA/SingularValueDecomposition.php

@@ -0,0 +1,526 @@
+<?php
+/**
+ *	@package JAMA
+ *
+ *	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, sigma[$k] = S[$k][$k], are ordered so that
+ *	sigma[0] >= sigma[1] >= ... >= sigma[n-1].
+ *
+ *	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.
+ *
+ *	@author  Paul Meagher
+ *	@license PHP v3.0
+ *	@version 1.1
+ */
+class SingularValueDecomposition  {
+
+	/**
+	 *	Internal storage of U.
+	 *	@var array
+	 */
+	private $U = array();
+
+	/**
+	 *	Internal storage of V.
+	 *	@var array
+	 */
+	private $V = array();
+
+	/**
+	 *	Internal storage of singular values.
+	 *	@var array
+	 */
+	private $s = array();
+
+	/**
+	 *	Row dimension.
+	 *	@var int
+	 */
+	private $m;
+
+	/**
+	 *	Column dimension.
+	 *	@var int
+	 */
+	private $n;
+
+
+	/**
+	 *	Construct the singular value decomposition
+	 *
+	 *	Derived from LINPACK code.
+	 *
+	 *	@param $A Rectangular matrix
+	 *	@return Structure to access U, S and V.
+	 */
+	public function __construct($Arg) {
+
+		// Initialize.
+		$A = $Arg->getArrayCopy();
+		$this->m = $Arg->getRowDimension();
+		$this->n = $Arg->getColumnDimension();
+		$nu      = min($this->m, $this->n);
+		$e       = array();
+		$work    = array();
+		$wantu   = true;
+		$wantv   = true;
+		$nct = min($this->m - 1, $this->n);
+		$nrt = max(0, min($this->n - 2, $this->m));
+
+		// Reduce A to bidiagonal form, storing the diagonal elements
+		// in s and the super-diagonal elements in e.
+		for ($k = 0; $k < max($nct,$nrt); ++$k) {
+
+			if ($k < $nct) {
+				// Compute the transformation for the k-th column and
+				// place the k-th diagonal in s[$k].
+				// Compute 2-norm of k-th column without under/overflow.
+				$this->s[$k] = 0;
+				for ($i = $k; $i < $this->m; ++$i) {
+					$this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
+				}
+				if ($this->s[$k] != 0.0) {
+					if ($A[$k][$k] < 0.0) {
+						$this->s[$k] = -$this->s[$k];
+					}
+					for ($i = $k; $i < $this->m; ++$i) {
+						$A[$i][$k] /= $this->s[$k];
+					}
+					$A[$k][$k] += 1.0;
+				}
+				$this->s[$k] = -$this->s[$k];
+			}
+
+			for ($j = $k + 1; $j < $this->n; ++$j) {
+				if (($k < $nct) & ($this->s[$k] != 0.0)) {
+					// Apply the transformation.
+					$t = 0;
+					for ($i = $k; $i < $this->m; ++$i) {
+						$t += $A[$i][$k] * $A[$i][$j];
+					}
+					$t = -$t / $A[$k][$k];
+					for ($i = $k; $i < $this->m; ++$i) {
+						$A[$i][$j] += $t * $A[$i][$k];
+					}
+					// Place the k-th row of A into e for the
+					// subsequent calculation of the row transformation.
+					$e[$j] = $A[$k][$j];
+				}
+			}
+
+			if ($wantu AND ($k < $nct)) {
+				// Place the transformation in U for subsequent back
+				// multiplication.
+				for ($i = $k; $i < $this->m; ++$i) {
+					$this->U[$i][$k] = $A[$i][$k];
+				}
+			}
+
+			if ($k < $nrt) {
+				// Compute the k-th row transformation and place the
+				// k-th super-diagonal in e[$k].
+				// Compute 2-norm without under/overflow.
+				$e[$k] = 0;
+				for ($i = $k + 1; $i < $this->n; ++$i) {
+					$e[$k] = hypo($e[$k], $e[$i]);
+				}
+				if ($e[$k] != 0.0) {
+					if ($e[$k+1] < 0.0) {
+						$e[$k] = -$e[$k];
+					}
+					for ($i = $k + 1; $i < $this->n; ++$i) {
+						$e[$i] /= $e[$k];
+					}
+					$e[$k+1] += 1.0;
+				}
+				$e[$k] = -$e[$k];
+				if (($k+1 < $this->m) AND ($e[$k] != 0.0)) {
+					// Apply the transformation.
+					for ($i = $k+1; $i < $this->m; ++$i) {
+						$work[$i] = 0.0;
+					}
+					for ($j = $k+1; $j < $this->n; ++$j) {
+						for ($i = $k+1; $i < $this->m; ++$i) {
+							$work[$i] += $e[$j] * $A[$i][$j];
+						}
+					}
+					for ($j = $k + 1; $j < $this->n; ++$j) {
+						$t = -$e[$j] / $e[$k+1];
+						for ($i = $k + 1; $i < $this->m; ++$i) {
+							$A[$i][$j] += $t * $work[$i];
+						}
+					}
+				}
+				if ($wantv) {
+					// Place the transformation in V for subsequent
+					// back multiplication.
+					for ($i = $k + 1; $i < $this->n; ++$i) {
+						$this->V[$i][$k] = $e[$i];
+					}
+				}
+			}
+		}
+
+		// Set up the final bidiagonal matrix or order p.
+		$p = min($this->n, $this->m + 1);
+		if ($nct < $this->n) {
+			$this->s[$nct] = $A[$nct][$nct];
+		}
+		if ($this->m < $p) {
+			$this->s[$p-1] = 0.0;
+		}
+		if ($nrt + 1 < $p) {
+			$e[$nrt] = $A[$nrt][$p-1];
+		}
+		$e[$p-1] = 0.0;
+		// If required, generate U.
+		if ($wantu) {
+			for ($j = $nct; $j < $nu; ++$j) {
+				for ($i = 0; $i < $this->m; ++$i) {
+					$this->U[$i][$j] = 0.0;
+				}
+				$this->U[$j][$j] = 1.0;
+			}
+			for ($k = $nct - 1; $k >= 0; --$k) {
+				if ($this->s[$k] != 0.0) {
+					for ($j = $k + 1; $j < $nu; ++$j) {
+						$t = 0;
+						for ($i = $k; $i < $this->m; ++$i) {
+							$t += $this->U[$i][$k] * $this->U[$i][$j];
+						}
+						$t = -$t / $this->U[$k][$k];
+						for ($i = $k; $i < $this->m; ++$i) {
+							$this->U[$i][$j] += $t * $this->U[$i][$k];
+						}
+					}
+					for ($i = $k; $i < $this->m; ++$i ) {
+						$this->U[$i][$k] = -$this->U[$i][$k];
+					}
+					$this->U[$k][$k] = 1.0 + $this->U[$k][$k];
+					for ($i = 0; $i < $k - 1; ++$i) {
+						$this->U[$i][$k] = 0.0;
+					}
+				} else {
+					for ($i = 0; $i < $this->m; ++$i) {
+						$this->U[$i][$k] = 0.0;
+					}
+					$this->U[$k][$k] = 1.0;
+				}
+			}
+		}
+
+		// If required, generate V.
+		if ($wantv) {
+			for ($k = $this->n - 1; $k >= 0; --$k) {
+				if (($k < $nrt) AND ($e[$k] != 0.0)) {
+					for ($j = $k + 1; $j < $nu; ++$j) {
+						$t = 0;
+						for ($i = $k + 1; $i < $this->n; ++$i) {
+							$t += $this->V[$i][$k]* $this->V[$i][$j];
+						}
+						$t = -$t / $this->V[$k+1][$k];
+						for ($i = $k + 1; $i < $this->n; ++$i) {
+							$this->V[$i][$j] += $t * $this->V[$i][$k];
+						}
+					}
+				}
+				for ($i = 0; $i < $this->n; ++$i) {
+					$this->V[$i][$k] = 0.0;
+				}
+				$this->V[$k][$k] = 1.0;
+			}
+		}
+
+		// Main iteration loop for the singular values.
+		$pp   = $p - 1;
+		$iter = 0;
+		$eps  = pow(2.0, -52.0);
+
+		while ($p > 0) {
+			// 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<p
+			// kase = 2  if s(k) is negligible and k<p
+			// kase = 3  if e[k-1] is negligible, k<p, and
+			//           s(k), ..., s(p) are not negligible (qr step).
+			// kase = 4  if e(p-1) is negligible (convergence).
+			for ($k = $p - 2; $k >= -1; --$k) {
+				if ($k == -1) {
+					break;
+				}
+				if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) {
+					$e[$k] = 0.0;
+					break;
+				}
+			}
+			if ($k == $p - 2) {
+				$kase = 4;
+			} else {
+				for ($ks = $p - 1; $ks >= $k; --$ks) {
+					if ($ks == $k) {
+						break;
+					}
+					$t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.);
+					if (abs($this->s[$ks]) <= $eps * $t)  {
+						$this->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:
+						$f = $e[$p-2];
+						$e[$p-2] = 0.0;
+						for ($j = $p - 2; $j >= $k; --$j) {
+							$t  = hypo($this->s[$j],$f);
+							$cs = $this->s[$j] / $t;
+							$sn = $f / $t;
+							$this->s[$j] = $t;
+							if ($j != $k) {
+								$f = -$sn * $e[$j-1];
+								$e[$j-1] = $cs * $e[$j-1];
+							}
+							if ($wantv) {
+								for ($i = 0; $i < $this->n; ++$i) {
+									$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1];
+									$this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1];
+									$this->V[$i][$j] = $t;
+								}
+							}
+						}
+						break;
+				// Split at negligible s(k).
+				case 2:
+						$f = $e[$k-1];
+						$e[$k-1] = 0.0;
+						for ($j = $k; $j < $p; ++$j) {
+							$t = hypo($this->s[$j], $f);
+							$cs = $this->s[$j] / $t;
+							$sn = $f / $t;
+							$this->s[$j] = $t;
+							$f = -$sn * $e[$j];
+							$e[$j] = $cs * $e[$j];
+							if ($wantu) {
+								for ($i = 0; $i < $this->m; ++$i) {
+									$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1];
+									$this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1];
+									$this->U[$i][$j] = $t;
+								}
+							}
+						}
+						break;
+				// Perform one qr step.
+				case 3:
+						// Calculate the shift.
+						$scale = max(max(max(max(
+									abs($this->s[$p-1]),abs($this->s[$p-2])),abs($e[$p-2])),
+									abs($this->s[$k])), abs($e[$k]));
+						$sp   = $this->s[$p-1] / $scale;
+						$spm1 = $this->s[$p-2] / $scale;
+						$epm1 = $e[$p-2] / $scale;
+						$sk   = $this->s[$k] / $scale;
+						$ek   = $e[$k] / $scale;
+						$b    = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
+						$c    = ($sp * $epm1) * ($sp * $epm1);
+						$shift = 0.0;
+						if (($b != 0.0) || ($c != 0.0)) {
+							$shift = sqrt($b * $b + $c);
+							if ($b < 0.0) {
+								$shift = -$shift;
+							}
+							$shift = $c / ($b + $shift);
+						}
+						$f = ($sk + $sp) * ($sk - $sp) + $shift;
+						$g = $sk * $ek;
+						// Chase zeros.
+						for ($j = $k; $j < $p-1; ++$j) {
+							$t  = hypo($f,$g);
+							$cs = $f/$t;
+							$sn = $g/$t;
+							if ($j != $k) {
+								$e[$j-1] = $t;
+							}
+							$f = $cs * $this->s[$j] + $sn * $e[$j];
+							$e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
+							$g = $sn * $this->s[$j+1];
+							$this->s[$j+1] = $cs * $this->s[$j+1];
+							if ($wantv) {
+								for ($i = 0; $i < $this->n; ++$i) {
+									$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1];
+									$this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1];
+									$this->V[$i][$j] = $t;
+								}
+							}
+							$t = hypo($f,$g);
+							$cs = $f/$t;
+							$sn = $g/$t;
+							$this->s[$j] = $t;
+							$f = $cs * $e[$j] + $sn * $this->s[$j+1];
+							$this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1];
+							$g = $sn * $e[$j+1];
+							$e[$j+1] = $cs * $e[$j+1];
+							if ($wantu && ($j < $this->m - 1)) {
+								for ($i = 0; $i < $this->m; ++$i) {
+									$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1];
+									$this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1];
+									$this->U[$i][$j] = $t;
+								}
+							}
+						}
+						$e[$p-2] = $f;
+						$iter = $iter + 1;
+						break;
+				// Convergence.
+				case 4:
+						// Make the singular values positive.
+						if ($this->s[$k] <= 0.0) {
+							$this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);
+							if ($wantv) {
+								for ($i = 0; $i <= $pp; ++$i) {
+									$this->V[$i][$k] = -$this->V[$i][$k];
+								}
+							}
+						}
+						// Order the singular values.
+						while ($k < $pp) {
+							if ($this->s[$k] >= $this->s[$k+1]) {
+								break;
+							}
+							$t = $this->s[$k];
+							$this->s[$k] = $this->s[$k+1];
+							$this->s[$k+1] = $t;
+							if ($wantv AND ($k < $this->n - 1)) {
+								for ($i = 0; $i < $this->n; ++$i) {
+									$t = $this->V[$i][$k+1];
+									$this->V[$i][$k+1] = $this->V[$i][$k];
+									$this->V[$i][$k] = $t;
+								}
+							}
+							if ($wantu AND ($k < $this->m-1)) {
+								for ($i = 0; $i < $this->m; ++$i) {
+									$t = $this->U[$i][$k+1];
+									$this->U[$i][$k+1] = $this->U[$i][$k];
+									$this->U[$i][$k] = $t;
+								}
+							}
+							++$k;
+						}
+						$iter = 0;
+						--$p;
+						break;
+			} // end switch
+		} // end while
+
+	} // end constructor
+
+
+	/**
+	 *	Return the left singular vectors
+	 *
+	 *	@access public
+	 *	@return U
+	 */
+	public function getU() {
+		return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
+	}
+
+
+	/**
+	 *	Return the right singular vectors
+	 *
+	 *	@access public
+	 *	@return V
+	 */
+	public function getV() {
+		return new Matrix($this->V);
+	}
+
+
+	/**
+	 *	Return the one-dimensional array of singular values
+	 *
+	 *	@access public
+	 *	@return diagonal of S.
+	 */
+	public function getSingularValues() {
+		return $this->s;
+	}
+
+
+	/**
+	 *	Return the diagonal matrix of singular values
+	 *
+	 *	@access public
+	 *	@return S
+	 */
+	public function getS() {
+		for ($i = 0; $i < $this->n; ++$i) {
+			for ($j = 0; $j < $this->n; ++$j) {
+				$S[$i][$j] = 0.0;
+			}
+			$S[$i][$i] = $this->s[$i];
+		}
+		return new Matrix($S);
+	}
+
+
+	/**
+	 *	Two norm
+	 *
+	 *	@access public
+	 *	@return max(S)
+	 */
+	public function norm2() {
+		return $this->s[0];
+	}
+
+
+	/**
+	 *	Two norm condition number
+	 *
+	 *	@access public
+	 *	@return max(S)/min(S)
+	 */
+	public function cond() {
+		return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
+	}
+
+
+	/**
+	 *	Effective numerical matrix rank
+	 *
+	 *	@access public
+	 *	@return Number of nonnegligible singular values.
+	 */
+	public function rank() {
+		$eps = pow(2.0, -52.0);
+		$tol = max($this->m, $this->n) * $this->s[0] * $eps;
+		$r = 0;
+		for ($i = 0; $i < count($this->s); ++$i) {
+			if ($this->s[$i] > $tol) {
+				++$r;
+			}
+		}
+		return $r;
+	}
+
+}	//	class SingularValueDecomposition

include/ECM/PHPExcel/Shared/JAMA/examples/LMQuadTest.php

@@ -0,0 +1,116 @@
+<?php
+/**
+ * quadratic (p-o)'S'S(p-o)
+ * solve for o, S
+ * S is a single scale factor
+ */
+class LMQuadTest {
+
+	/**
+	 * @param array[] $x
+	 * @param array[] $a
+	 */
+	function val($x, $a) {
+		if (count($a) != 3) die ("Wrong number of elements in array a");
+		if (count($x) != 2) die ("Wrong number of elements in array x");
+
+		$ox = $a[0];
+		$oy = $a[1];
+		$s  = $a[2];
+
+		$sdx = $s * ($x[0] - $ox);
+		$sdy = $s * ($x[1] - $oy);
+
+		return ($sdx * $sdx) + ($sdy * $sdy);
+   }	//	function val()
+
+
+	/**
+	 * z = (p-o)'S'S(p-o)
+	 * dz/dp = 2S'S(p-o)
+	 *
+	 * z = (s*(px-ox))^2 + (s*(py-oy))^2
+	 * dz/dox = -2(s*(px-ox))*s
+	 * dz/ds = 2*s*[(px-ox)^2 + (py-oy)^2]
+	 *
+	 * z = (s*dx)^2 + (s*dy)^2
+	 * dz/ds = 2(s*dx)*dx + 2(s*dy)*dy
+	 *
+	 * @param array[] $x
+	 * @param array[] $a
+	 * @param int $a_k
+	 * @param array[] $a
+	 */
+	function grad($x, $a, $a_k) {
+		if (count($a) != 3) die ("Wrong number of elements in array a");
+		if (count($x) != 2) die ("Wrong number of elements in array x");
+		if ($a_k < 3) die ("a_k=".$a_k);
+
+		$ox = $a[0];
+		$oy = $a[1];
+		$s  = $a[2];
+
+		$dx = ($x[0] - $ox);
+		$dy = ($x[1] - $oy);
+
+		if ($a_k == 0)
+			return -2.*$s*$s*$dx;
+		elseif ($a_k == 1)
+			return -2.*$s*$s*$dy;
+		else
+			return 2.*$s*($dx*$dx + $dy*$dy);
+	}	//	function grad()
+
+
+	/**
+	 * @return array[] $a
+	 */
+	function initial() {
+		$a[0] = 0.05;
+		$a[1] = 0.1;
+		$a[2] = 1.0;
+
+		return $a;
+	}	//	function initial()
+
+
+	/**
+	 * @return Object[] $a
+	 */
+	function testdata() {
+		$npts = 25;
+
+		$a[0] = 0.;
+		$a[1] = 0.;
+		$a[2] = 0.9;
+
+		$i = 0;
+
+		for ($r = -2; $r <= 2; ++$r) {
+			for ($c = -2; $c <= 2; ++$c) {
+				$x[$i][0] = $c;
+				$x[$i][1] = $r;
+				$y[$i] = $this->val($x[$i], $a);
+				print("Quad ".$c.",".$r." -> ".$y[$i]."<br />");
+				$s[$i] = 1.;
+				++$i;
+			}
+		}
+		print("quad x= ");
+
+		$qx = new Matrix($x);
+		$qx->print(10, 2);
+
+		print("quad y= ");
+		$qy = new Matrix($y, $npts);
+		$qy->print(10, 2);
+
+		$o[0] = $x;
+		$o[1] = $a;
+		$o[2] = $y;
+		$o[3] = $s;
+
+		return $o;
+	}	//	function testdata()
+
+}	//	class LMQuadTest

include/ECM/PHPExcel/Shared/JAMA/examples/LagrangeInterpolation.php

@@ -0,0 +1,59 @@
+<?php
+
+require_once "../Matrix.php";
+
+/**
+ * Given n points (x0,y0)...(xn-1,yn-1), the following methid computes
+ * the polynomial factors of the n-1't degree polynomial passing through
+ * the n points.
+ *
+ * Example: Passing in three points (2,3) (1,4) and (3,7) will produce
+ * the results [2.5, -8.5, 10] which means that the points are on the
+ * curve y = 2.5x² - 8.5x + 10.
+ *
+ * @see http://geosoft.no/software/lagrange/LagrangeInterpolation.java.html
+ * @author Jacob Dreyer
+ * @author Paul Meagher (port to PHP and minor changes)
+ *
+ * @param x[] float
+ * @param y[] float
+ */
+class LagrangeInterpolation {
+
+	public function findPolynomialFactors($x, $y) {
+		$n = count($x);
+
+		$data = array();  // double[n][n];
+		$rhs  = array();  // double[n];
+
+		for ($i = 0; $i < $n; ++$i) {
+			$v = 1;
+			for ($j = 0; $j < $n; ++$j) {
+				$data[$i][$n-$j-1] = $v;
+				$v *= $x[$i];
+			}
+			$rhs[$i] = $y[$i];
+		}
+
+		// Solve m * s = b
+		$m = new Matrix($data);
+		$b = new Matrix($rhs, $n);
+
+		$s = $m->solve($b);
+
+		return $s->getRowPackedCopy();
+	}	//	function findPolynomialFactors()
+
+}	//	class LagrangeInterpolation
+
+
+$x = array(2.0, 1.0, 3.0);
+$y = array(3.0, 4.0, 7.0);
+
+$li = new LagrangeInterpolation;
+$f = $li->findPolynomialFactors($x, $y);
+
+
+for ($i = 0; $i < 3; ++$i) {
+	echo $f[$i]."<br />";
+}

include/ECM/PHPExcel/Shared/JAMA/examples/LagrangeInterpolation2.php

@@ -0,0 +1,59 @@
+<?php
+
+require_once "../Matrix.php";
+
+/**
+ * Given n points (x0,y0)...(xn-1,yn-1), the following method computes
+ * the polynomial factors of the n-1't degree polynomial passing through
+ * the n points.
+ *
+ * Example: Passing in three points (2,3) (1,4) and (3,7) will produce
+ * the results [2.5, -8.5, 10] which means that the points are on the
+ * curve y = 2.5x² - 8.5x + 10.
+ *
+ * @see http://geosoft.no/software/lagrange/LagrangeInterpolation.java.html
+ * @see http://source.freehep.org/jcvsweb/ilc/LCSIM/wdview/lcsim/src/org/lcsim/fit/polynomial/PolynomialFitter.java
+ * @author Jacob Dreyer
+ * @author Paul Meagher (port to PHP and minor changes)
+ *
+ * @param x[] float
+ * @param y[] float
+ */
+class LagrangeInterpolation {
+
+	public function findPolynomialFactors($x, $y) {
+		$n = count($x);
+
+		$data = array();  // double[n][n];
+		$rhs  = array();  // double[n];
+
+		for ($i = 0; $i < $n; ++$i) {
+			$v = 1;
+			for ($j = 0; $j < $n; ++$j) {
+				$data[$i][$n-$j-1] = $v;
+				$v *= $x[$i];
+			}
+			$rhs[$i] = $y[$i];
+		}
+
+		// Solve m * s = b
+		$m = new Matrix($data);
+		$b = new Matrix($rhs, $n);
+
+		$s = $m->solve($b);
+
+		return $s->getRowPackedCopy();
+	}	//	function findPolynomialFactors()
+
+}	//	class LagrangeInterpolation
+
+
+$x = array(2.0, 1.0, 3.0);
+$y = array(3.0, 4.0, 7.0);
+
+$li = new LagrangeInterpolation;
+$f = $li->findPolynomialFactors($x, $y);
+
+for ($i = 0; $i < 3; ++$i) {
+	echo $f[$i]."<br />";
+}

include/ECM/PHPExcel/Shared/JAMA/examples/LevenbergMarquardt.php

@@ -0,0 +1,185 @@
+<?php
+
+// Levenberg-Marquardt in PHP
+
+// http://www.idiom.com/~zilla/Computer/Javanumeric/LM.java
+
+class LevenbergMarquardt {
+
+	/**
+	 * Calculate the current sum-squared-error
+	 *
+	 * Chi-squared is the distribution of squared Gaussian errors,
+	 * thus the name.
+	 *
+	 * @param double[][] $x
+	 * @param double[] $a
+	 * @param double[] $y,
+	 * @param double[] $s,
+	 * @param object $f
+	 */
+	function chiSquared($x, $a, $y, $s, $f) {
+		$npts = count($y);
+		$sum = 0.0;
+
+		for ($i = 0; $i < $npts; ++$i) {
+			$d = $y[$i] - $f->val($x[$i], $a);
+			$d = $d / $s[$i];
+			$sum = $sum + ($d*$d);
+		}
+
+		return $sum;
+	}	//	function chiSquared()
+
+
+	/**
+	 * Minimize E = sum {(y[k] - f(x[k],a)) / s[k]}^2
+	 * The individual errors are optionally scaled by s[k].
+	 * Note that LMfunc implements the value and gradient of f(x,a),
+	 * NOT the value and gradient of E with respect to a!
+	 *
+	 * @param x array of domain points, each may be multidimensional
+	 * @param y corresponding array of values
+	 * @param a the parameters/state of the model
+	 * @param vary false to indicate the corresponding a[k] is to be held fixed
+	 * @param s2 sigma^2 for point i
+	 * @param lambda blend between steepest descent (lambda high) and
+	 *	jump to bottom of quadratic (lambda zero).
+	 * 	Start with 0.001.
+	 * @param termepsilon termination accuracy (0.01)
+	 * @param maxiter	stop and return after this many iterations if not done
+	 * @param verbose	set to zero (no prints), 1, 2
+	 *
+	 * @return the new lambda for future iterations.
+	 *  Can use this and maxiter to interleave the LM descent with some other
+	 *  task, setting maxiter to something small.
+	 */
+	function solve($x, $a, $y, $s, $vary, $f, $lambda, $termepsilon, $maxiter, $verbose) {
+		$npts = count($y);
+		$nparm = count($a);
+
+		if ($verbose > 0) {
+			print("solve x[".count($x)."][".count($x[0])."]");
+			print(" a[".count($a)."]");
+			println(" y[".count(length)."]");
+		}
+
+		$e0 = $this->chiSquared($x, $a, $y, $s, $f);
+
+		//double lambda = 0.001;
+		$done = false;
+
+		// g = gradient, H = hessian, d = step to minimum
+		// H d = -g, solve for d
+		$H = array();
+		$g = array();
+
+		//double[] d = new double[nparm];
+
+		$oos2 = array();
+
+		for($i = 0; $i < $npts; ++$i) {
+			$oos2[$i] = 1./($s[$i]*$s[$i]);
+		}
+		$iter = 0;
+		$term = 0;	// termination count test
+
+		do {
+			++$iter;
+
+			// hessian approximation
+			for( $r = 0; $r < $nparm; ++$r) {
+				for( $c = 0; $c < $nparm; ++$c) {
+					for( $i = 0; $i < $npts; ++$i) {
+						if ($i == 0) $H[$r][$c] = 0.;
+						$xi = $x[$i];
+						$H[$r][$c] += ($oos2[$i] * $f->grad($xi, $a, $r) * $f->grad($xi, $a, $c));
+					}  //npts
+				} //c
+			} //r
+
+			// boost diagonal towards gradient descent
+			for( $r = 0; $r < $nparm; ++$r)
+				$H[$r][$r] *= (1. + $lambda);
+
+			// gradient
+			for( $r = 0; $r < $nparm; ++$r) {
+				for( $i = 0; $i < $npts; ++$i) {
+					if ($i == 0) $g[$r] = 0.;
+					$xi = $x[$i];
+					$g[$r] += ($oos2[$i] * ($y[$i]-$f->val($xi,$a)) * $f->grad($xi, $a, $r));
+				}
+			} //npts
+
+			// scale (for consistency with NR, not necessary)
+			if ($false) {
+				for( $r = 0; $r < $nparm; ++$r) {
+					$g[$r] = -0.5 * $g[$r];
+					for( $c = 0; $c < $nparm; ++$c) {
+						$H[$r][$c] *= 0.5;
+					}
+				}
+			}
+
+			// solve H d = -g, evaluate error at new location
+			//double[] d = DoubleMatrix.solve(H, g);
+//			double[] d = (new Matrix(H)).lu().solve(new Matrix(g, nparm)).getRowPackedCopy();
+			//double[] na = DoubleVector.add(a, d);
+//			double[] na = (new Matrix(a, nparm)).plus(new Matrix(d, nparm)).getRowPackedCopy();
+//			double e1 = chiSquared(x, na, y, s, f);
+
+//			if (verbose > 0) {
+//				System.out.println("\n\niteration "+iter+" lambda = "+lambda);
+//				System.out.print("a = ");
+//				(new Matrix(a, nparm)).print(10, 2);
+//				if (verbose > 1) {
+//					System.out.print("H = ");
+//					(new Matrix(H)).print(10, 2);
+//					System.out.print("g = ");
+//					(new Matrix(g, nparm)).print(10, 2);
+//					System.out.print("d = ");
+//					(new Matrix(d, nparm)).print(10, 2);
+//				}
+//				System.out.print("e0 = " + e0 + ": ");
+//				System.out.print("moved from ");
+//				(new Matrix(a, nparm)).print(10, 2);
+//				System.out.print("e1 = " + e1 + ": ");
+//				if (e1 < e0) {
+//					System.out.print("to ");
+//					(new Matrix(na, nparm)).print(10, 2);
+//				} else {
+//					System.out.println("move rejected");
+//				}
+//			}
+
+			// termination test (slightly different than NR)
+//			if (Math.abs(e1-e0) > termepsilon) {
+//				term = 0;
+//			} else {
+//				term++;
+//				if (term == 4) {
+//					System.out.println("terminating after " + iter + " iterations");
+//					done = true;
+//				}
+//			}
+//			if (iter >= maxiter) done = true;
+
+			// in the C++ version, found that changing this to e1 >= e0
+			// was not a good idea.  See comment there.
+			//
+//			if (e1 > e0 || Double.isNaN(e1)) { // new location worse than before
+//				lambda *= 10.;
+//			} else {		// new location better, accept new parameters
+//				lambda *= 0.1;
+//				e0 = e1;
+//				// simply assigning a = na will not get results copied back to caller
+//				for( int i = 0; i < nparm; i++ ) {
+//					if (vary[i]) a[i] = na[i];
+//				}
+//			}
+		} while(!$done);
+
+		return $lambda;
+	}	//	function solve()
+
+}	//	class LevenbergMarquardt

include/ECM/PHPExcel/Shared/JAMA/examples/MagicSquareExample.php

@@ -0,0 +1,182 @@
+<?php
+/**
+* @package JAMA
+*/
+
+require_once "../Matrix.php";
+
+/**
+* Example of use of Matrix Class, featuring magic squares.
+*/
+class MagicSquareExample {
+
+  /**
+  * Generate magic square test matrix.
+  * @param int n dimension of matrix
+  */
+  function magic($n) {
+
+    // Odd order
+
+    if (($n % 2) == 1) {
+      $a = ($n+1)/2;
+      $b = ($n+1);
+      for ($j = 0; $j < $n; ++$j)
+        for ($i = 0; $i < $n; ++$i)
+          $M[$i][$j] = $n*(($i+$j+$a) % $n) + (($i+2*$j+$b) % $n) + 1;
+
+    // Doubly Even Order
+
+    } else if (($n % 4) == 0) {
+      for ($j = 0; $j < $n; ++$j) {
+        for ($i = 0; $i < $n; ++$i) {
+          if ((($i+1)/2)%2 == (($j+1)/2)%2)
+            $M[$i][$j] = $n*$n-$n*$i-$j;
+          else
+            $M[$i][$j] = $n*$i+$j+1;
+        }
+      }
+
+    // Singly Even Order
+
+    } else {
+
+      $p = $n/2;
+      $k = ($n-2)/4;
+      $A = $this->magic($p);
+      $M = array();
+      for ($j = 0; $j < $p; ++$j) {
+        for ($i = 0; $i < $p; ++$i) {
+          $aij = $A->get($i,$j);
+          $M[$i][$j]       = $aij;
+          $M[$i][$j+$p]    = $aij + 2*$p*$p;
+          $M[$i+$p][$j]    = $aij + 3*$p*$p;
+          $M[$i+$p][$j+$p] = $aij + $p*$p;
+        }
+      }
+
+      for ($i = 0; $i < $p; ++$i) {
+        for ($j = 0; $j < $k; ++$j) {
+          $t = $M[$i][$j];
+          $M[$i][$j] = $M[$i+$p][$j];
+          $M[$i+$p][$j] = $t;
+        }
+        for ($j = $n-$k+1; $j < $n; ++$j) {
+          $t = $M[$i][$j];
+          $M[$i][$j] = $M[$i+$p][$j];
+          $M[$i+$p][$j] = $t;
+        }
+      }
+
+      $t = $M[$k][0];  $M[$k][0]  = $M[$k+$p][0];  $M[$k+$p][0]  = $t;
+      $t = $M[$k][$k]; $M[$k][$k] = $M[$k+$p][$k]; $M[$k+$p][$k] = $t;
+
+    }
+
+    return new Matrix($M);
+
+  }
+
+  /**
+  * Simple function to replicate PHP 5 behaviour
+  */
+  function microtime_float() {
+    list($usec, $sec) = explode(" ", microtime());
+    return ((float)$usec + (float)$sec);
+  }
+
+  /**
+  * Tests LU, QR, SVD and symmetric Eig decompositions.
+  *
+  *   n       = order of magic square.
+  *   trace   = diagonal sum, should be the magic sum, (n^3 + n)/2.
+  *   max_eig = maximum eigenvalue of (A + A')/2, should equal trace.
+  *   rank    = linear algebraic rank, should equal n if n is odd,
+  *             be less than n if n is even.
+  *   cond    = L_2 condition number, ratio of singular values.
+  *   lu_res  = test of LU factorization, norm1(L*U-A(p,:))/(n*eps).
+  *   qr_res  = test of QR factorization, norm1(Q*R-A)/(n*eps).
+  */
+  function main() {
+    ?>
+    <p>Test of Matrix Class, using magic squares.</p>
+    <p>See MagicSquareExample.main() for an explanation.</p>
+    <table border='1' cellspacing='0' cellpadding='4'>
+      <tr>
+        <th>n</th>
+        <th>trace</th>
+        <th>max_eig</th>
+        <th>rank</th>
+        <th>cond</th>
+        <th>lu_res</th>
+        <th>qr_res</th>
+      </tr>
+      <?php
+      $start_time = $this->microtime_float();
+      $eps = pow(2.0,-52.0);
+      for ($n = 3; $n <= 6; ++$n) {
+        echo "<tr>";
+
+        echo "<td align='right'>$n</td>";
+
+        $M = $this->magic($n);
+        $t = (int) $M->trace();
+
+        echo "<td align='right'>$t</td>";
+
+        $O = $M->plus($M->transpose());
+        $E = new EigenvalueDecomposition($O->times(0.5));
+        $d = $E->getRealEigenvalues();
+
+        echo "<td align='right'>".$d[$n-1]."</td>";
+
+        $r = $M->rank();
+
+        echo "<td align='right'>".$r."</td>";
+
+        $c = $M->cond();
+
+        if ($c < 1/$eps)
+          echo "<td align='right'>".sprintf("%.3f",$c)."</td>";
+        else
+          echo "<td align='right'>Inf</td>";
+
+        $LU = new LUDecomposition($M);
+        $L = $LU->getL();
+        $U = $LU->getU();
+        $p = $LU->getPivot();
+        // Java version: R = L.times(U).minus(M.getMatrix(p,0,n-1));
+        $S = $L->times($U);
+        $R = $S->minus($M->getMatrix($p,0,$n-1));
+        $res = $R->norm1()/($n*$eps);
+
+        echo "<td align='right'>".sprintf("%.3f",$res)."</td>";
+
+        $QR = new QRDecomposition($M);
+        $Q = $QR->getQ();
+        $R = $QR->getR();
+        $S = $Q->times($R);
+        $R = $S->minus($M);
+        $res = $R->norm1()/($n*$eps);
+
+        echo "<td align='right'>".sprintf("%.3f",$res)."</td>";
+
+        echo "</tr>";
+
+     }
+     echo "<table>";
+     echo "<br />";
+
+     $stop_time = $this->microtime_float();
+     $etime = $stop_time - $start_time;
+
+     echo "<p>Elapsed time is ". sprintf("%.4f",$etime) ." seconds.</p>";
+
+  }
+
+}
+
+$magic = new MagicSquareExample();
+$magic->main();
+
+?>

include/ECM/PHPExcel/Shared/JAMA/examples/benchmark.php

@@ -0,0 +1,263 @@
+<?php
+
+error_reporting(E_ALL);
+
+/**
+ * @package JAMA
+ */
+
+require_once '../Matrix.php';
+require_once 'Stats.php';
+
+
+/**
+ * Example of use of Matrix Class, featuring magic squares.
+ */
+class Benchmark {
+	public $stat;
+
+
+	/**
+	 * Simple function to replicate PHP 5 behaviour
+	 */
+	function microtime_float() {
+		list($usec, $sec) = explode(" ", microtime());
+
+		return ((float)$usec + (float)$sec);
+	}	//	function microtime_float()
+
+
+	function displayStats($times = null) {
+		$this->stat->setData($times);
+		$stats = $this->stat->calcFull();
+
+		echo '<table style="margin-left:32px;">';
+		echo '<tr><td style="text-align:right;"><b>n:</b><td style="text-align:right;">' . $stats['count'] . ' </td></tr>';
+		echo '<tr><td style="text-align:right;"><b>Mean:</b><td style="text-align:right;">' . $stats['mean'] . ' </td></tr>';
+		echo '<tr><td style="text-align:right;"><b>Min.:</b><td style="text-align:right;">' . $stats['min'] . ' </td></tr>';
+		echo '<tr><td style="text-align:right;"><b>Max.:</b><td style="text-align:right;">' . $stats['max'] . ' </td></tr>';
+		echo '<tr><td style="text-align:right;"><b>&sigma;:</b><td style="text-align:right;">' . $stats['stdev'] . ' </td></tr>';
+		echo '<tr><td style="text-align:right;"><b>Variance:</b><td style="text-align:right;">' . $stats['variance'] . ' </td></tr>';
+		echo '<tr><td style="text-align:right;"><b>Range:</b><td style="text-align:right;">' . $stats['range'] . ' </td></tr>';
+		echo '</table>';
+
+		return $stats;
+	}	//	function displayStats()
+
+
+	function runEig($n = 4, $t = 100) {
+		$times = array();
+
+		for ($i = 0; $i < $t; ++$i) {
+			$M = Matrix::random($n, $n);
+			$start_time = $this->microtime_float();
+			$E = new EigenvalueDecomposition($M);
+			$stop_time = $this->microtime_float();
+			$times[] = $stop_time - $start_time;
+		}
+
+		return $times;
+	}	//	function runEig()
+
+
+	function runLU($n = 4, $t = 100) {
+		$times = array();
+
+		for ($i = 0; $i < $t; ++$i) {
+			$M = Matrix::random($n, $n);
+			$start_time = $this->microtime_float();
+			$E = new LUDecomposition($M);
+			$stop_time = $this->microtime_float();
+			$times[] = $stop_time - $start_time;
+		}
+
+		return $times;
+	}	//	function runLU()
+
+
+	function runQR($n = 4, $t = 100) {
+		$times = array();
+
+		for ($i = 0; $i < $t; ++$i) {
+			$M = Matrix::random($n, $n);
+			$start_time = $this->microtime_float();
+			$E = new QRDecomposition($M);
+			$stop_time = $this->microtime_float();
+			$times[] = $stop_time - $start_time;
+		}
+
+		return $times;
+	}	//	function runQR()
+
+
+	function runCholesky($n = 4, $t = 100) {
+		$times = array();
+
+		for ($i = 0; $i < $t; ++$i) {
+			$M = Matrix::random($n, $n);
+			$start_time = $this->microtime_float();
+			$E = new CholeskyDecomposition($M);
+			$stop_time = $this->microtime_float();
+			$times[] = $stop_time - $start_time;
+		}
+
+		return $times;
+	}	//	function runCholesky()
+
+
+	function runSVD($n = 4, $t = 100) {
+		$times = array();
+
+		for ($i = 0; $i < $t; ++$i) {
+			$M = Matrix::random($n, $n);
+			$start_time = $this->microtime_float();
+			$E = new SingularValueDecomposition($M);
+			$stop_time = $this->microtime_float();
+			$times[] = $stop_time - $start_time;
+		}
+
+		return $times;
+	}	//	function runSVD()
+
+
+	function run() {
+		$n = 8;
+		$t = 16;
+		$sum = 0;
+		echo "<b>Cholesky decomposition: $t random {$n}x{$n} matrices</b><br />";
+		$r = $this->displayStats($this->runCholesky($n, $t));
+		$sum += $r['mean'] * $n;
+
+		echo '<hr />';
+
+		echo "<b>Eigenvalue decomposition: $t random {$n}x{$n} matrices</b><br />";
+		$r = $this->displayStats($this->runEig($n, $t));
+		$sum += $r['mean'] * $n;
+
+		echo '<hr />';
+
+		echo "<b>LU decomposition: $t random {$n}x{$n} matrices</b><br />";
+		$r = $this->displayStats($this->runLU($n, $t));
+		$sum += $r['mean'] * $n;
+
+		echo '<hr />';
+
+		echo "<b>QR decomposition: $t random {$n}x{$n} matrices</b><br />";
+		$r = $this->displayStats($this->runQR($n, $t));
+		$sum += $r['mean'] * $n;
+
+		echo '<hr />';
+
+		echo "<b>Singular Value decomposition: $t random {$n}x{$n} matrices</b><br />";
+		$r = $this->displayStats($this->runSVD($n, $t));
+		$sum += $r['mean'] * $n;
+
+		return $sum;
+	}	//	function run()
+
+
+	public function __construct() {
+		$this->stat = new Base();
+	}	//	function Benchmark()
+
+}  // class Benchmark		(end MagicSquareExample)
+
+
+$benchmark = new Benchmark();
+
+switch($_REQUEST['decomposition']) {
+	case 'cholesky':
+		$m = array();
+		for ($i = 2; $i <= 8; $i *= 2) {
+			$t = 32 / $i;
+			echo "<b>Cholesky decomposition: $t random {$i}x{$i} matrices</b><br />";
+			$s = $benchmark->displayStats($benchmark->runCholesky($i, $t));
+			$m[$i] = $s['mean'];
+			echo "<br />";
+		}
+		echo '<pre>';
+		foreach($m as $x => $y) {
+			echo "$x\t" . 1000*$y . "\n";
+		}
+		echo '</pre>';
+		break;
+	case 'eigenvalue':
+		$m = array();
+		for ($i = 2; $i <= 8; $i *= 2) {
+			$t = 32 / $i;
+			echo "<b>Eigenvalue decomposition: $t random {$i}x{$i} matrices</b><br />";
+			$s = $benchmark->displayStats($benchmark->runEig($i, $t));
+			$m[$i] = $s['mean'];
+			echo "<br />";
+		}
+		echo '<pre>';
+		foreach($m as $x => $y) {
+			echo "$x\t" . 1000*$y . "\n";
+		}
+		echo '</pre>';
+		break;
+	case 'lu':
+		$m = array();
+		for ($i = 2; $i <= 8; $i *= 2) {
+			$t = 32 / $i;
+			echo "<b>LU decomposition: $t random {$i}x{$i} matrices</b><br />";
+			$s = $benchmark->displayStats($benchmark->runLU($i, $t));
+			$m[$i] = $s['mean'];
+			echo "<br />";
+		}
+		echo '<pre>';
+		foreach($m as $x => $y) {
+			echo "$x\t" . 1000*$y . "\n";
+		}
+		echo '</pre>';
+		break;
+	case 'qr':
+		$m = array();
+		for ($i = 2; $i <= 8; $i *= 2) {
+			$t = 32 / $i;
+			echo "<b>QR decomposition: $t random {$i}x{$i} matrices</b><br />";
+			$s = $benchmark->displayStats($benchmark->runQR($i, $t));
+			$m[$i] = $s['mean'];
+			echo "<br />";
+		}
+		echo '<pre>';
+		foreach($m as $x => $y) {
+			echo "$x\t" . 1000*$y . "\n";
+		}
+		echo '</pre>';
+		break;
+	case 'svd':
+		$m = array();
+		for($i = 2; $i <= 8; $i *= 2) {
+			$t = 32 / $i;
+			echo "<b>Singular value decomposition: $t random {$i}x{$i} matrices</b><br />";
+			$s = $benchmark->displayStats($benchmark->runSVD($i, $t));
+			$m[$i] = $s['mean'];
+			echo "<br />";
+		}
+		echo '<pre>';
+		foreach($m as $x => $y) {
+			echo "$x\t" . 1000*$y . "\n";
+		}
+		echo '</pre>';
+		break;
+	case 'all':
+		$s = $benchmark->run();
+		print("<br /><b>Total<b>: {$s}s<br />");
+		break;
+	default:
+		?>
+		<ul>
+			<li><a href="benchmark.php?decomposition=all">Complete Benchmark</a>
+				<ul>
+					<li><a href="benchmark.php?decomposition=cholesky">Cholesky</a></li>
+					<li><a href="benchmark.php?decomposition=eigenvalue">Eigenvalue</a></li>
+					<li><a href="benchmark.php?decomposition=lu">LU</a></li>
+					<li><a href="benchmark.php?decomposition=qr">QR</a></li>
+					<li><a href="benchmark.php?decomposition=svd">Singular Value</a></li>
+				</ul>
+			</li>
+		</ul>
+		<?php
+		break;
+}

include/ECM/PHPExcel/Shared/JAMA/examples/polyfit.php

@@ -0,0 +1,73 @@
+<?php
+require_once "../Matrix.php";
+/*
+* @package JAMA
+* @author Michael Bommarito
+* @author Paul Meagher
+* @version 0.1
+*
+* Function to fit an order n polynomial function through
+* a series of x-y data points using least squares.
+*
+* @param $X array x values
+* @param $Y array y values
+* @param $n int order of polynomial to be used for fitting
+* @returns array $coeffs of polynomial coefficients
+* Pre-Conditions: the system is not underdetermined: sizeof($X) > $n+1
+*/
+function polyfit($X, $Y, $n) {
+	for ($i = 0; $i < sizeof($X); ++$i)
+		for ($j = 0; $j <= $n; ++$j)
+			$A[$i][$j] = pow($X[$i], $j);
+	for ($i=0; $i < sizeof($Y); ++$i)
+		$B[$i] = array($Y[$i]);
+	$matrixA = new Matrix($A);
+	$matrixB = new Matrix($B);
+	$C = $matrixA->solve($matrixB);
+	return $C->getMatrix(0, $n, 0, 1);
+}
+
+function printpoly( $C = null ) {
+	for($i = $C->m - 1; $i >= 0; --$i) {
+		$r = $C->get($i, 0);
+		if ( abs($r) <= pow(10, -9) )
+			$r = 0;
+		if ($i == $C->m - 1)
+			echo $r . "x<sup>$i</sup>";
+		else if ($i < $C->m - 1)
+			echo " + " . $r . "x<sup>$i</sup>";
+		else if ($i == 0)
+			echo " + " . $r;
+	}
+}
+
+$X = array(0,1,2,3,4,5);
+$Y = array(4,3,12,67,228, 579);
+$points = new Matrix(array($X, $Y));
+$points->toHTML();
+printpoly(polyfit($X, $Y, 4));
+
+echo '<hr />';
+
+$X = array(0,1,2,3,4,5);
+$Y = array(1,2,5,10,17, 26);
+$points = new Matrix(array($X, $Y));
+$points->toHTML();
+printpoly(polyfit($X, $Y, 2));
+
+echo '<hr />';
+
+$X = array(0,1,2,3,4,5,6);
+$Y = array(-90,-104,-178,-252,-26, 1160, 4446);
+$points = new Matrix(array($X, $Y));
+$points->toHTML();
+printpoly(polyfit($X, $Y, 5));
+
+echo '<hr />';
+
+$X = array(0,1,2,3,4);
+$Y = array(mt_rand(0, 10), mt_rand(40, 80), mt_rand(240, 400), mt_rand(1800, 2215), mt_rand(8000, 9000));
+$points = new Matrix(array($X, $Y));
+$points->toHTML();
+printpoly(polyfit($X, $Y, 3));
+?>

include/ECM/PHPExcel/Shared/JAMA/examples/tile.php

@@ -0,0 +1,78 @@
+<?php
+
+include "../Matrix.php";
+
+/**
+* Tiling of matrix X in [rowWise by colWise] dimension. Tiling
+* creates a larger matrix than the original data X. Example, if
+* X is to be tiled in a [3 x 4] manner, then:
+*
+*     /            \
+*     | X  X  X  X |
+* C = | X  X  X  X |
+*     | X  X  X  X |
+*     \           /
+*
+* @param X Matrix
+* @param rowWise int
+* @param colWise int
+* @return Matrix
+*/
+
+function tile(&$X, $rowWise, $colWise){
+
+  $xArray = $X->getArray();
+  print_r($xArray);
+
+  $countRow    = 0;
+  $countColumn = 0;
+
+  $m = $X->getRowDimension();
+  $n = $X->getColumnDimension();
+
+  if( $rowWise<1 || $colWise<1 ){
+    die("tile : Array index is out-of-bound.");
+  }
+
+  $newRowDim = $m*$rowWise;
+  $newColDim = $n*$colWise;
+
+  $result = array();
+
+  for($i=0 ; $i<$newRowDim; ++$i) {
+
+    $holder = array();
+
+    for($j=0 ; $j<$newColDim ; ++$j) {
+
+      $holder[$j] = $xArray[$countRow][$countColumn++];
+
+      // reset the column-index to zero to avoid reference to out-of-bound index in xArray[][]
+
+      if($countColumn == $n) { $countColumn = 0; }
+
+    } // end for
+
+    ++$countRow;
+
+    // reset the row-index to zero to avoid reference to out-of-bound index in xArray[][]
+
+    if($countRow == $m) { $countRow = 0; }
+
+    $result[$i] = $holder;
+
+  } // end for
+
+  return new Matrix($result);
+
+}
+
+
+$X =array(1,2,3,4,5,6,7,8,9);
+$nRow = 3;
+$nCol = 3;
+$tiled_matrix = tile(new Matrix($X), $nRow, $nCol);
+echo "<pre>";
+print_r($tiled_matrix);
+echo "</pre>";
+?>

include/ECM/PHPExcel/Shared/JAMA/tests/TestMatrix.php

@@ -0,0 +1,415 @@
+<?php
+
+require_once "../Matrix.php";
+
+class TestMatrix {
+
+  function TestMatrix() {
+
+    // define test variables
+
+    $errorCount   = 0;
+    $warningCount = 0;
+    $columnwise   = array(1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.);
+    $rowwise      = array(1.,4.,7.,10.,2.,5.,8.,11.,3.,6.,9.,12.);
+    $avals        = array(array(1.,4.,7.,10.),array(2.,5.,8.,11.),array(3.,6.,9.,12.));
+    $rankdef      = $avals;
+    $tvals        = array(array(1.,2.,3.),array(4.,5.,6.),array(7.,8.,9.),array(10.,11.,12.));
+    $subavals     = array(array(5.,8.,11.),array(6.,9.,12.));
+    $rvals        = array(array(1.,4.,7.),array(2.,5.,8.,11.),array(3.,6.,9.,12.));
+    $pvals        = array(array(1.,1.,1.),array(1.,2.,3.),array(1.,3.,6.));
+    $ivals        = array(array(1.,0.,0.,0.),array(0.,1.,0.,0.),array(0.,0.,1.,0.));
+    $evals        = array(array(0.,1.,0.,0.),array(1.,0.,2.e-7,0.),array(0.,-2.e-7,0.,1.),array(0.,0.,1.,0.));
+    $square       = array(array(166.,188.,210.),array(188.,214.,240.),array(210.,240.,270.));
+    $sqSolution   = array(array(13.),array(15.));
+    $condmat      = array(array(1.,3.),array(7.,9.));
+    $rows         = 3;
+    $cols         = 4;
+    $invalidID    = 5; /* should trigger bad shape for construction with val        */
+    $raggedr      = 0; /* (raggedr,raggedc) should be out of bounds in ragged array */
+    $raggedc      = 4;
+    $validID      = 3; /* leading dimension of intended test Matrices               */
+    $nonconformld = 4; /* leading dimension which is valid, but nonconforming       */
+    $ib           = 1; /* index ranges for sub Matrix                               */
+    $ie           = 2;
+    $jb           = 1;
+    $je           = 3;
+    $rowindexset       = array(1,2);
+    $badrowindexset    = array(1,3);
+    $columnindexset    = array(1,2,3);
+    $badcolumnindexset = array(1,2,4);
+    $columnsummax      = 33.;
+    $rowsummax         = 30.;
+    $sumofdiagonals    = 15;
+    $sumofsquares      = 650;
+
+    /**
+    * Test matrix methods
+    */
+
+    /**
+    * Constructors and constructor-like methods:
+    *
+    *   Matrix(double[], int)
+    *   Matrix(double[][])
+    *   Matrix(int, int)
+    *   Matrix(int, int, double)
+    *   Matrix(int, int, double[][])
+    *   constructWithCopy(double[][])
+    *   random(int,int)
+    *   identity(int)
+    */
+    echo "<p>Testing constructors and constructor-like methods...</p>";
+
+    $A = new Matrix($columnwise, 3);
+    if($A instanceof Matrix) {
+      $this->try_success("Column-packed constructor...");
+    } else
+      $errorCount = $this->try_failure($errorCount, "Column-packed constructor...", "Unable to construct Matrix");
+
+    $T = new Matrix($tvals);
+    if($T instanceof Matrix)
+      $this->try_success("2D array constructor...");
+    else
+      $errorCount = $this->try_failure($errorCount, "2D array constructor...", "Unable to construct Matrix");
+
+    $A = new Matrix($columnwise, $validID);
+    $B = new Matrix($avals);
+    $tmp = $B->get(0,0);
+    $avals[0][0] = 0.0;
+    $C = $B->minus($A);
+    $avals[0][0] = $tmp;
+    $B = Matrix::constructWithCopy($avals);
+    $tmp = $B->get(0,0);
+    $avals[0][0] = 0.0;
+    /** check that constructWithCopy behaves properly **/
+    if ( ( $tmp - $B->get(0,0) ) != 0.0 )
+      $errorCount = $this->try_failure($errorCount,"constructWithCopy... ","copy not effected... data visible outside");
+    else
+      $this->try_success("constructWithCopy... ","");
+
+    $I = new Matrix($ivals);
+    if ( $this->checkMatrices($I,Matrix::identity(3,4)) )
+      $this->try_success("identity... ","");
+    else
+      $errorCount = $this->try_failure($errorCount,"identity... ","identity Matrix not successfully created");
+
+    /**
+    * Access Methods:
+    *
+    *   getColumnDimension()
+    *   getRowDimension()
+    *   getArray()
+    *   getArrayCopy()
+    *   getColumnPackedCopy()
+    *   getRowPackedCopy()
+    *   get(int,int)
+    *   getMatrix(int,int,int,int)
+    *   getMatrix(int,int,int[])
+    *   getMatrix(int[],int,int)
+    *   getMatrix(int[],int[])
+    *   set(int,int,double)
+    *   setMatrix(int,int,int,int,Matrix)
+    *   setMatrix(int,int,int[],Matrix)
+    *   setMatrix(int[],int,int,Matrix)
+    *   setMatrix(int[],int[],Matrix)
+    */
+    print "<p>Testing access methods...</p>";
+
+	$B = new Matrix($avals);
+	if($B->getRowDimension() == $rows)
+	  $this->try_success("getRowDimension...");
+	else
+	  $errorCount = $this->try_failure($errorCount, "getRowDimension...");
+
+	if($B->getColumnDimension() == $cols)
+	  $this->try_success("getColumnDimension...");
+	else
+	  $errorCount = $this->try_failure($errorCount, "getColumnDimension...");
+
+	$barray = $B->getArray();
+	if($this->checkArrays($barray, $avals))
+	  $this->try_success("getArray...");
+	else
+	  $errorCount = $this->try_failure($errorCount, "getArray...");
+
+	$bpacked = $B->getColumnPackedCopy();
+	if($this->checkArrays($bpacked, $columnwise))
+	  $this->try_success("getColumnPackedCopy...");
+	else
+	  $errorCount = $this->try_failure($errorCount, "getColumnPackedCopy...");
+
+	$bpacked = $B->getRowPackedCopy();
+	if($this->checkArrays($bpacked, $rowwise))
+	  $this->try_success("getRowPackedCopy...");
+	else
+	  $errorCount = $this->try_failure($errorCount, "getRowPackedCopy...");
+
+    /**
+    * Array-like methods:
+    *   minus
+    *   minusEquals
+    *   plus
+    *   plusEquals
+    *   arrayLeftDivide
+    *   arrayLeftDivideEquals
+    *   arrayRightDivide
+    *   arrayRightDivideEquals
+    *   arrayTimes
+    *   arrayTimesEquals
+    *   uminus
+    */
+    print "<p>Testing array-like methods...</p>";
+
+    /**
+    * I/O methods:
+    *   read
+    *   print
+    *   serializable:
+    *   writeObject
+    *   readObject
+    */
+    print "<p>Testing I/O methods...</p>";
+
+    /**
+    * Test linear algebra methods
+    */
+    echo "<p>Testing linear algebra methods...<p>";
+
+    $A = new Matrix($columnwise, 3);
+    if( $this->checkMatrices($A->transpose(), $T) )
+      $this->try_success("Transpose check...");
+    else
+      $errorCount = $this->try_failure($errorCount, "Transpose check...", "Matrices are not equal");
+
+    if($this->checkScalars($A->norm1(), $columnsummax))
+      $this->try_success("Maximum column sum...");
+    else
+      $errorCount = $this->try_failure($errorCount, "Maximum column sum...", "Incorrect: " . $A->norm1() . " != " . $columnsummax);
+
+    if($this->checkScalars($A->normInf(), $rowsummax))
+      $this->try_success("Maximum row sum...");
+    else
+      $errorCount = $this->try_failure($errorCount, "Maximum row sum...", "Incorrect: " . $A->normInf() . " != " . $rowsummax );
+
+    if($this->checkScalars($A->normF(), sqrt($sumofsquares)))
+      $this->try_success("Frobenius norm...");
+    else
+      $errorCount = $this->try_failure($errorCount, "Frobenius norm...", "Incorrect:" . $A->normF() . " != " . sqrt($sumofsquares));
+
+    if($this->checkScalars($A->trace(), $sumofdiagonals))
+      $this->try_success("Matrix trace...");
+    else
+      $errorCount = $this->try_failure($errorCount, "Matrix trace...", "Incorrect: " . $A->trace() . " != " . $sumofdiagonals);
+
+    $B = $A->getMatrix(0, $A->getRowDimension(), 0, $A->getRowDimension());
+    if( $B->det() == 0 )
+      $this->try_success("Matrix determinant...");
+    else
+      $errorCount = $this->try_failure($errorCount, "Matrix determinant...", "Incorrect: " . $B->det() . " != " . 0);
+
+    $A = new Matrix($columnwise,3);
+    $SQ = new Matrix($square);
+    if ($this->checkMatrices($SQ, $A->times($A->transpose())))
+      $this->try_success("times(Matrix)...");
+    else {
+      $errorCount = $this->try_failure($errorCount, "times(Matrix)...", "Unable to multiply matrices");
+      $SQ->toHTML();
+      $AT->toHTML();
+    }
+
+    $A = new Matrix($columnwise, 4);
+
+    $QR = $A->qr();
+    $R = $QR->getR();
+    $Q = $QR->getQ();
+    if($this->checkMatrices($A, $Q->times($R)))
+      $this->try_success("QRDecomposition...","");
+    else
+      $errorCount = $this->try_failure($errorCount,"QRDecomposition...","incorrect qr decomposition calculation");
+
+    $A = new Matrix($columnwise, 4);
+    $SVD = $A->svd();
+    $U = $SVD->getU();
+    $S = $SVD->getS();
+    $V = $SVD->getV();
+    if ($this->checkMatrices($A, $U->times($S->times($V->transpose()))))
+      $this->try_success("SingularValueDecomposition...","");
+    else
+      $errorCount = $this->try_failure($errorCount,"SingularValueDecomposition...","incorrect singular value decomposition calculation");
+
+    $n = $A->getColumnDimension();
+    $A = $A->getMatrix(0,$n-1,0,$n-1);
+    $A->set(0,0,0.);
+
+    $LU = $A->lu();
+    $L  = $LU->getL();
+    if ( $this->checkMatrices($A->getMatrix($LU->getPivot(),0,$n-1), $L->times($LU->getU())) )
+      $this->try_success("LUDecomposition...","");
+    else
+      $errorCount = $this->try_failure($errorCount,"LUDecomposition...","incorrect LU decomposition calculation");
+
+    $X = $A->inverse();
+    if ( $this->checkMatrices($A->times($X),Matrix::identity(3,3)) )
+       $this->try_success("inverse()...","");
+     else
+       $errorCount = $this->try_failure($errorCount, "inverse()...","incorrect inverse calculation");
+
+    $DEF = new Matrix($rankdef);
+    if($this->checkScalars($DEF->rank(), min($DEF->getRowDimension(), $DEF->getColumnDimension())-1))
+      $this->try_success("Rank...");
+    else
+      $this->try_failure("Rank...", "incorrect rank calculation");
+
+    $B = new Matrix($condmat);
+    $SVD = $B->svd();
+    $singularvalues = $SVD->getSingularValues();
+    if($this->checkScalars($B->cond(), $singularvalues[0]/$singularvalues[min($B->getRowDimension(), $B->getColumnDimension())-1]))
+      $this->try_success("Condition number...");
+    else
+      $this->try_failure("Condition number...", "incorrect condition number calculation");
+
+    $SUB = new Matrix($subavals);
+    $O   = new Matrix($SUB->getRowDimension(),1,1.0);
+    $SOL = new Matrix($sqSolution);
+    $SQ = $SUB->getMatrix(0,$SUB->getRowDimension()-1,0,$SUB->getRowDimension()-1);
+    if ( $this->checkMatrices($SQ->solve($SOL),$O) )
+      $this->try_success("solve()...","");
+    else
+     $errorCount = $this->try_failure($errorCount,"solve()...","incorrect lu solve calculation");
+
+    $A = new Matrix($pvals);
+    $Chol = $A->chol();
+    $L = $Chol->getL();
+    if ( $this->checkMatrices($A, $L->times($L->transpose())) )
+      $this->try_success("CholeskyDecomposition...","");
+    else
+      $errorCount = $this->try_failure($errorCount,"CholeskyDecomposition...","incorrect Cholesky decomposition calculation");
+
+    $X = $Chol->solve(Matrix::identity(3,3));
+    if ( $this->checkMatrices($A->times($X), Matrix::identity(3,3)) )
+      $this->try_success("CholeskyDecomposition solve()...","");
+    else
+      $errorCount = $this->try_failure($errorCount,"CholeskyDecomposition solve()...","incorrect Choleskydecomposition solve calculation");
+
+    $Eig = $A->eig();
+    $D = $Eig->getD();
+    $V = $Eig->getV();
+    if( $this->checkMatrices($A->times($V),$V->times($D)) )
+      $this->try_success("EigenvalueDecomposition (symmetric)...","");
+    else
+      $errorCount = $this->try_failure($errorCount,"EigenvalueDecomposition (symmetric)...","incorrect symmetric Eigenvalue decomposition calculation");
+
+    $A = new Matrix($evals);
+    $Eig = $A->eig();
+    $D = $Eig->getD();
+    $V = $Eig->getV();
+    if ( $this->checkMatrices($A->times($V),$V->times($D)) )
+      $this->try_success("EigenvalueDecomposition (nonsymmetric)...","");
+    else
+      $errorCount = $this->try_failure($errorCount,"EigenvalueDecomposition (nonsymmetric)...","incorrect nonsymmetric Eigenvalue decomposition calculation");
+
+	print("<b>{$errorCount} total errors</b>.");
+  }
+
+  /**
+  * Print appropriate messages for successful outcome try
+  * @param string $s
+  * @param string $e
+  */
+  function try_success($s, $e = "") {
+    print "> ". $s ."success<br />";
+    if ($e != "")
+      print "> Message: ". $e ."<br />";
+  }
+
+  /**
+  * Print appropriate messages for unsuccessful outcome try
+  * @param int $count
+  * @param string $s
+  * @param string $e
+  * @return int incremented counter
+  */
+  function try_failure($count, $s, $e="") {
+    print "> ". $s ."*** failure ***<br />> Message: ". $e ."<br />";
+    return ++$count;
+  }
+
+  /**
+  * Print appropriate messages for unsuccessful outcome try
+  * @param int $count
+  * @param string $s
+  * @param string $e
+  * @return int incremented counter
+  */
+  function try_warning($count, $s, $e="") {
+    print "> ". $s ."*** warning ***<br />> Message: ". $e ."<br />";
+    return ++$count;
+  }
+
+  /**
+  * Check magnitude of difference of "scalars".
+  * @param float $x
+  * @param float $y
+  */
+  function checkScalars($x, $y) {
+    $eps = pow(2.0,-52.0);
+    if ($x == 0 & abs($y) < 10*$eps) return;
+    if ($y == 0 & abs($x) < 10*$eps) return;
+    if (abs($x-$y) > 10 * $eps * max(abs($x),abs($y)))
+      return false;
+    else
+      return true;
+  }
+
+  /**
+  * Check norm of difference of "vectors".
+  * @param float $x[]
+  * @param float $y[]
+  */
+  function checkVectors($x, $y) {
+    $nx = count($x);
+    $ny = count($y);
+    if ($nx == $ny)
+      for($i=0; $i < $nx; ++$i)
+        $this->checkScalars($x[$i],$y[$i]);
+    else
+      die("Attempt to compare vectors of different lengths");
+  }
+
+  /**
+  * Check norm of difference of "arrays".
+  * @param float $x[][]
+  * @param float $y[][]
+  */
+  function checkArrays($x, $y) {
+    $A = new Matrix($x);
+    $B = new Matrix($y);
+    return $this->checkMatrices($A,$B);
+  }
+
+  /**
+  * Check norm of difference of "matrices".
+  * @param matrix $X
+  * @param matrix $Y
+  */
+  function checkMatrices($X = null, $Y = null) {
+    if( $X == null || $Y == null )
+      return false;
+
+    $eps = pow(2.0,-52.0);
+    if ($X->norm1() == 0. & $Y->norm1() < 10*$eps) return true;
+    if ($Y->norm1() == 0. & $X->norm1() < 10*$eps) return true;
+
+    $A = $X->minus($Y);
+
+    if ($A->norm1() > 1000 * $eps * max($X->norm1(),$Y->norm1()))
+      die("The norm of (X-Y) is too large: ".$A->norm1());
+    else
+      return true;
+  }
+
+}
+
+$test = new TestMatrix;
+?>

include/ECM/PHPExcel/Shared/JAMA/utils/Error.php

@@ -0,0 +1,82 @@
+<?php
+/**
+ *	@package JAMA
+ *
+ *	Error handling
+ *	@author Michael Bommarito
+ *	@version 01292005
+ */
+
+//Language constant
+define('JAMALANG', 'EN');
+
+
+//All errors may be defined by the following format:
+//define('ExceptionName', N);
+//$error['lang'][ExceptionName] = 'Error message';
+$error = array();
+
+/*
+I've used Babelfish and a little poor knowledge of Romance/Germanic languages for the translations here.
+Feel free to correct anything that looks amiss to you.
+*/
+
+define('PolymorphicArgumentException', -1);
+$error['EN'][PolymorphicArgumentException] = "Invalid argument pattern for polymorphic function.";
+$error['FR'][PolymorphicArgumentException] = "Modèle inadmissible d'argument pour la fonction polymorphe.".
+$error['DE'][PolymorphicArgumentException] = "Unzulässiges Argumentmuster für polymorphe Funktion.";
+
+define('ArgumentTypeException', -2);
+$error['EN'][ArgumentTypeException] = "Invalid argument type.";
+$error['FR'][ArgumentTypeException] = "Type inadmissible d'argument.";
+$error['DE'][ArgumentTypeException] = "Unzulässige Argumentart.";
+
+define('ArgumentBoundsException', -3);
+$error['EN'][ArgumentBoundsException] = "Invalid argument range.";
+$error['FR'][ArgumentBoundsException] = "Gamme inadmissible d'argument.";
+$error['DE'][ArgumentBoundsException] = "Unzulässige Argumentstrecke.";
+
+define('MatrixDimensionException', -4);
+$error['EN'][MatrixDimensionException] = "Matrix dimensions are not equal.";
+$error['FR'][MatrixDimensionException] = "Les dimensions de Matrix ne sont pas égales.";
+$error['DE'][MatrixDimensionException] = "Matrixmaße sind nicht gleich.";
+
+define('PrecisionLossException', -5);
+$error['EN'][PrecisionLossException] = "Significant precision loss detected.";
+$error['FR'][PrecisionLossException] = "Perte significative de précision détectée.";
+$error['DE'][PrecisionLossException] = "Bedeutender Präzision Verlust ermittelte.";
+
+define('MatrixSPDException', -6);
+$error['EN'][MatrixSPDException] = "Can only perform operation on symmetric positive definite matrix.";
+$error['FR'][MatrixSPDException] = "Perte significative de précision détectée.";
+$error['DE'][MatrixSPDException] = "Bedeutender Präzision Verlust ermittelte.";
+
+define('MatrixSingularException', -7);
+$error['EN'][MatrixSingularException] = "Can only perform operation on singular matrix.";
+
+define('MatrixRankException', -8);
+$error['EN'][MatrixRankException] = "Can only perform operation on full-rank matrix.";
+
+define('ArrayLengthException', -9);
+$error['EN'][ArrayLengthException] = "Array length must be a multiple of m.";
+
+define('RowLengthException', -10);
+$error['EN'][RowLengthException] = "All rows must have the same length.";
+
+/**
+ *	Custom error handler
+ *	@param int $num Error number
+ */
+function JAMAError($errorNumber = null) {
+	global $error;
+
+	if (isset($errorNumber)) {
+		if (isset($error[JAMALANG][$errorNumber])) {
+			return $error[JAMALANG][$errorNumber];
+		} else {
+			return $error['EN'][$errorNumber];
+		}
+	} else {
+		return ("Invalid argument to JAMAError()");
+	}
+}

include/ECM/PHPExcel/Shared/JAMA/utils/Maths.php

@@ -0,0 +1,43 @@
+<?php
+/**
+ *	@package JAMA
+ *
+ *	Pythagorean Theorem:
+ *
+ *	a = 3
+ *	b = 4
+ *	r = sqrt(square(a) + square(b))
+ *	r = 5
+ *
+ *	r = sqrt(a^2 + b^2) without under/overflow.
+ */
+function hypo($a, $b) {
+	if (abs($a) > abs($b)) {
+		$r = $b / $a;
+		$r = abs($a) * sqrt(1 + $r * $r);
+	} elseif ($b != 0) {
+		$r = $a / $b;
+		$r = abs($b) * sqrt(1 + $r * $r);
+	} else {
+		$r = 0.0;
+	}
+	return $r;
+}	//	function hypo()
+
+
+/**
+ *	Mike Bommarito's version.
+ *	Compute n-dimensional hyotheneuse.
+ *
+function hypot() {
+	$s = 0;
+	foreach (func_get_args() as $d) {
+		if (is_numeric($d)) {
+			$s += pow($d, 2);
+		} else {
+			throw new Exception(JAMAError(ArgumentTypeException));
+		}
+	}
+	return sqrt($s);
+}
+*/