8889841cPK     Õ,[³ˆ\µï  ï    utils/Maths.phpnu „[µü¤        <?php

/**
 *    Pythagorean Theorem:.
 *
 *    a = 3
 *    b = 4
 *    r = sqrt(square(a) + square(b))
 *    r = 5
 *
 *    r = sqrt(a^2 + b^2) without under/overflow.
 *
 * @param mixed $a
 * @param mixed $b
 *
 * @return float
 */
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;
}
PK     Õ,[ùU/D  D    LUDecomposition.phpnu „[µü¤        <?php

namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;

use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;

/**
 *    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
 */
class LUDecomposition
{
    const MATRIX_SINGULAR_EXCEPTION = 'Can only perform operation on singular matrix.';
    const MATRIX_SQUARE_EXCEPTION = 'Mismatched Row dimension';

    /**
     * Decomposition storage.
     *
     * @var array
     */
    private $LU = [];

    /**
     * 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 = [];

    /**
     * LU Decomposition constructor.
     *
     * @param Matrix $A Rectangular matrix
     */
    public function __construct($A)
    {
        if ($A instanceof 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 = [];

            // 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 CalculationException(Matrix::ARGUMENT_TYPE_EXCEPTION);
        }
    }

    //    function __construct()

    /**
     * Get lower triangular factor.
     *
     * @return Matrix Lower triangular factor
     */
    public function getL()
    {
        $L = [];
        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 Matrix($L);
    }

    //    function getL()

    /**
     * Get upper triangular factor.
     *
     * @return Matrix Upper triangular factor
     */
    public function getU()
    {
        $U = [];
        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 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 bool 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 float
     */
    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;
        }

        throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
    }

    //    function det()

    /**
     * Solve A*X = B.
     *
     * @param Matrix $B a Matrix with as many rows as A and any number of columns
     *
     * @return Matrix X so that L*U*X = B(piv,:)
     */
    public function solve(Matrix $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;
            }

            throw new CalculationException(self::MATRIX_SINGULAR_EXCEPTION);
        }

        throw new CalculationException(self::MATRIX_SQUARE_EXCEPTION);
    }
}
PK     Õ,[^99  9    QRDecomposition.phpnu „[µü¤        <?php

namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;

use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;

/**
 *    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
 *
 * @version 1.1
 */
class QRDecomposition
{
    const MATRIX_RANK_EXCEPTION = 'Can only perform operation on full-rank matrix.';

    /**
     * Array for internal storage of decomposition.
     *
     * @var array
     */
    private $QR = [];

    /**
     * Row dimension.
     *
     * @var int
     */
    private $m;

    /**
     * Column dimension.
     *
     * @var int
     */
    private $n;

    /**
     * Array for internal storage of diagonal of R.
     *
     * @var array
     */
    private $Rdiag = [];

    /**
     * QR Decomposition computed by Householder reflections.
     *
     * @param Matrix $A Rectangular matrix
     */
    public function __construct(Matrix $A)
    {
        // Initialize.
        $this->QR = $A->getArray();
        $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;
        }
    }

    //    function __construct()

    /**
     *    Is the matrix full rank?
     *
     * @return bool 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()
    {
        $H = [];
        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 Matrix($H);
    }

    //    function getH()

    /**
     * Return the upper triangular factor.
     *
     * @return Matrix upper triangular factor
     */
    public function getR()
    {
        $R = [];
        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 Matrix($R);
    }

    //    function getR()

    /**
     * Generate and return the (economy-sized) orthogonal factor.
     *
     * @return Matrix orthogonal factor
     */
    public function getQ()
    {
        $Q = [];
        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];
                    }
                }
            }
        }

        return new 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(Matrix $B)
    {
        if ($B->getRowDimension() == $this->m) {
            if ($this->isFullRank()) {
                // Copy right hand side
                $nx = $B->getColumnDimension();
                $X = $B->getArray();
                // 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 Matrix($X);

                return $X->getMatrix(0, $this->n - 1, 0, $nx);
            }

            throw new CalculationException(self::MATRIX_RANK_EXCEPTION);
        }

        throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
    }
}
PK     Õ,[Møô  ô    CHANGELOG.TXTnu „[µü¤        Mar 1, 2005 11:15 AST by PM

+ For consistency, renamed Math.php to Maths.java, utils to util, 
  tests to test, docs to doc - 

+ Removed conditional logic from top of Matrix class.

+ Switched to using hypo function in Maths.php for all php-hypot calls.
  NOTE TO SELF: Need to make sure that all decompositions have been 
  switched over to using the bundled hypo.

Feb 25, 2005 at 10:00 AST by PM

+ Recommend using simpler Error.php instead of JAMA_Error.php but 
  can be persuaded otherwise.

PK     Õ,[úsR<“  “    EigenvalueDecomposition.phpnu „[µü¤        <?php

namespace PhpOffice\PhpSpreadsheet\Shared\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
 *
 * @version 1.1
 */
class EigenvalueDecomposition
{
    /**
     * Row and column dimension (square matrix).
     *
     * @var int
     */
    private $n;

    /**
     * Arrays for internal storage of eigenvalues.
     *
     * @var array
     */
    private $d = [];

    private $e = [];

    /**
     * Array for internal storage of eigenvectors.
     *
     * @var array
     */
    private $V = [];

    /**
     * Array for internal storage of nonsymmetric Hessenberg form.
     *
     * @var array
     */
    private $H = [];

    /**
     * Working storage for nonsymmetric algorithm.
     *
     * @var array
     */
    private $ort;

    /**
     * Used for complex scalar division.
     *
     * @var float
     */
    private $cdivr;

    private $cdivi;

    /**
     * @var array
     */
    private $A;

    /**
     * Symmetric Householder reduction to tridiagonal form.
     */
    private function tred2(): void
    {
        //  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];
        $j = 0;
        // 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 += $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.
     */
    private function tql2(): void
    {
        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 = 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;
                    // 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.
     */
    private function orthes(): void
    {
        $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.
     *
     * @param mixed $xr
     * @param mixed $xi
     * @param mixed $yr
     * @param mixed $yi
     */
    private function cdiv($xr, $xi, $yr, $yi): void
    {
        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.
     */
    private function hqr2(): void
    {
        //  Initialize
        $nn = $this->n;
        $n = $nn - 1;
        $low = 0;
        $high = $nn - 1;
        $eps = 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) || ($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
            } elseif ($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
                        $iMax = min($n, $k + 3);
                        for ($i = 0; $i <= $iMax; ++$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
            } elseif ($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;
                $kMax = min($j, $high);
                for ($k = $low; $k <= $kMax; ++$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.
     *
     * @param Matrix $Arg A Square matrix
     */
    public function __construct(Matrix $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 = [];
            // Reduce to Hessenberg form.
            $this->orthes();
            // Reduce Hessenberg to real Schur form.
            $this->hqr2();
        }
    }

    /**
     * Return the eigenvector matrix.
     *
     * @return Matrix V
     */
    public function getV()
    {
        return new Matrix($this->V, $this->n, $this->n);
    }

    /**
     * Return the real parts of the eigenvalues.
     *
     * @return array real(diag(D))
     */
    public function getRealEigenvalues()
    {
        return $this->d;
    }

    /**
     * Return the imaginary parts of the eigenvalues.
     *
     * @return array imag(diag(D))
     */
    public function getImagEigenvalues()
    {
        return $this->e;
    }

    /**
     * Return the block diagonal eigenvalue matrix.
     *
     * @return Matrix D
     */
    public function getD()
    {
        $D = [];
        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);
    }
}
PK     Õ,[|Úp      CholeskyDecomposition.phpnu „[µü¤        <?php

namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;

use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;

/**
 *    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
     */
    private $L = [];

    /**
     * Matrix row and column dimension.
     *
     * @var int
     */
    private $m;

    /**
     * Symmetric positive definite flag.
     *
     * @var bool
     */
    private $isspd = true;

    /**
     * CholeskyDecomposition.
     *
     *    Class constructor - decomposes symmetric positive definite matrix
     *
     * @param Matrix $A Matrix square symmetric positive definite matrix
     */
    public function __construct(Matrix $A)
    {
        $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;
            }
        }
    }

    /**
     *    Is the matrix symmetric and positive definite?
     *
     * @return bool
     */
    public function isSPD()
    {
        return $this->isspd;
    }

    /**
     * getL.
     *
     * Return triangular factor.
     *
     * @return Matrix Lower triangular matrix
     */
    public function getL()
    {
        return new Matrix($this->L);
    }

    /**
     * Solve A*X = B.
     *
     * @param Matrix $B Row-equal matrix
     *
     * @return Matrix L * L' * X = B
     */
    public function solve(Matrix $B)
    {
        if ($B->getRowDimension() == $this->m) {
            if ($this->isspd) {
                $X = $B->getArray();
                $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);
            }

            throw new CalculationException(Matrix::MATRIX_SPD_EXCEPTION);
        }

        throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
    }
}
PK     Õ,[ÙRïŒŽ  ŒŽ  
  Matrix.phpnu „[µü¤        <?php

namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;

use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
use PhpOffice\PhpSpreadsheet\Calculation\Functions;
use PhpOffice\PhpSpreadsheet\Shared\StringHelper;

/**
 * Matrix class.
 *
 * @author Paul Meagher
 * @author Michael Bommarito
 * @author Lukasz Karapuda
 * @author Bartek Matosiuk
 *
 * @version 1.8
 *
 * @see https://math.nist.gov/javanumerics/jama/
 */
class Matrix
{
    const POLYMORPHIC_ARGUMENT_EXCEPTION = 'Invalid argument pattern for polymorphic function.';
    const ARGUMENT_TYPE_EXCEPTION = 'Invalid argument type.';
    const ARGUMENT_BOUNDS_EXCEPTION = 'Invalid argument range.';
    const MATRIX_DIMENSION_EXCEPTION = 'Matrix dimensions are not equal.';
    const ARRAY_LENGTH_EXCEPTION = 'Array length must be a multiple of m.';
    const MATRIX_SPD_EXCEPTION = 'Can only perform operation on symmetric positive definite matrix.';

    /**
     * Matrix storage.
     *
     * @var array
     */
    public $A = [];

    /**
     * Matrix row dimension.
     *
     * @var int
     */
    private $m;

    /**
     * Matrix column dimension.
     *
     * @var int
     */
    private $n;

    /**
     * Polymorphic constructor.
     *
     * As PHP has no support for polymorphic constructors, we use tricks to make our own sort of polymorphism using func_num_args, func_get_arg, and gettype. In essence, we're just implementing a simple RTTI filter and calling the appropriate constructor.
     */
    public function __construct(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                //Rectangular matrix - m x n initialized from 2D array
                case 'array':
                    $this->m = count($args[0]);
                    $this->n = count($args[0][0]);
                    $this->A = $args[0];

                    break;
                //Square matrix - n x n
                case 'integer':
                    $this->m = $args[0];
                    $this->n = $args[0];
                    $this->A = array_fill(0, $this->m, array_fill(0, $this->n, 0));

                    break;
                //Rectangular matrix - m x n
                case 'integer,integer':
                    $this->m = $args[0];
                    $this->n = $args[1];
                    $this->A = array_fill(0, $this->m, array_fill(0, $this->n, 0));

                    break;
                //Rectangular matrix - m x n initialized from packed array
                case 'array,integer':
                    $this->m = $args[1];
                    if ($this->m != 0) {
                        $this->n = count($args[0]) / $this->m;
                    } else {
                        $this->n = 0;
                    }
                    if (($this->m * $this->n) == count($args[0])) {
                        for ($i = 0; $i < $this->m; ++$i) {
                            for ($j = 0; $j < $this->n; ++$j) {
                                $this->A[$i][$j] = $args[0][$i + $j * $this->m];
                            }
                        }
                    } else {
                        throw new CalculationException(self::ARRAY_LENGTH_EXCEPTION);
                    }

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
        } else {
            throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
        }
    }

    /**
     * getArray.
     *
     * @return array Matrix array
     */
    public function getArray()
    {
        return $this->A;
    }

    /**
     * getRowDimension.
     *
     * @return int Row dimension
     */
    public function getRowDimension()
    {
        return $this->m;
    }

    /**
     * getColumnDimension.
     *
     * @return int Column dimension
     */
    public function getColumnDimension()
    {
        return $this->n;
    }

    /**
     * get.
     *
     * Get the i,j-th element of the matrix.
     *
     * @param int $i Row position
     * @param int $j Column position
     *
     * @return float|int
     */
    public function get($i = null, $j = null)
    {
        return $this->A[$i][$j];
    }

    /**
     * getMatrix.
     *
     *    Get a submatrix
     *
     * @return Matrix Submatrix
     */
    public function getMatrix(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                //A($i0...; $j0...)
                case 'integer,integer':
                    [$i0, $j0] = $args;
                    if ($i0 >= 0) {
                        $m = $this->m - $i0;
                    } else {
                        throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
                    }
                    if ($j0 >= 0) {
                        $n = $this->n - $j0;
                    } else {
                        throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
                    }
                    $R = new self($m, $n);
                    for ($i = $i0; $i < $this->m; ++$i) {
                        for ($j = $j0; $j < $this->n; ++$j) {
                            $R->set($i, $j, $this->A[$i][$j]);
                        }
                    }

                    return $R;

                    break;
                //A($i0...$iF; $j0...$jF)
                case 'integer,integer,integer,integer':
                    [$i0, $iF, $j0, $jF] = $args;
                    if (($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0)) {
                        $m = $iF - $i0;
                    } else {
                        throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
                    }
                    if (($jF > $j0) && ($this->n >= $jF) && ($j0 >= 0)) {
                        $n = $jF - $j0;
                    } else {
                        throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
                    }
                    $R = new self($m + 1, $n + 1);
                    for ($i = $i0; $i <= $iF; ++$i) {
                        for ($j = $j0; $j <= $jF; ++$j) {
                            $R->set($i - $i0, $j - $j0, $this->A[$i][$j]);
                        }
                    }

                    return $R;

                    break;
                //$R = array of row indices; $C = array of column indices
                case 'array,array':
                    [$RL, $CL] = $args;
                    if (count($RL) > 0) {
                        $m = count($RL);
                    } else {
                        throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
                    }
                    if (count($CL) > 0) {
                        $n = count($CL);
                    } else {
                        throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
                    }
                    $R = new self($m, $n);
                    for ($i = 0; $i < $m; ++$i) {
                        for ($j = 0; $j < $n; ++$j) {
                            $R->set($i, $j, $this->A[$RL[$i]][$CL[$j]]);
                        }
                    }

                    return $R;

                    break;
                //A($i0...$iF); $CL = array of column indices
                case 'integer,integer,array':
                    [$i0, $iF, $CL] = $args;
                    if (($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0)) {
                        $m = $iF - $i0;
                    } else {
                        throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
                    }
                    if (count($CL) > 0) {
                        $n = count($CL);
                    } else {
                        throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
                    }
                    $R = new self($m, $n);
                    for ($i = $i0; $i < $iF; ++$i) {
                        for ($j = 0; $j < $n; ++$j) {
                            $R->set($i - $i0, $j, $this->A[$i][$CL[$j]]);
                        }
                    }

                    return $R;

                    break;
                //$RL = array of row indices
                case 'array,integer,integer':
                    [$RL, $j0, $jF] = $args;
                    if (count($RL) > 0) {
                        $m = count($RL);
                    } else {
                        throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
                    }
                    if (($jF >= $j0) && ($this->n >= $jF) && ($j0 >= 0)) {
                        $n = $jF - $j0;
                    } else {
                        throw new CalculationException(self::ARGUMENT_BOUNDS_EXCEPTION);
                    }
                    $R = new self($m, $n + 1);
                    for ($i = 0; $i < $m; ++$i) {
                        for ($j = $j0; $j <= $jF; ++$j) {
                            $R->set($i, $j - $j0, $this->A[$RL[$i]][$j]);
                        }
                    }

                    return $R;

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
        } else {
            throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
        }
    }

    /**
     * checkMatrixDimensions.
     *
     *    Is matrix B the same size?
     *
     * @param Matrix $B Matrix B
     *
     * @return bool
     */
    public function checkMatrixDimensions($B = null)
    {
        if ($B instanceof self) {
            if (($this->m == $B->getRowDimension()) && ($this->n == $B->getColumnDimension())) {
                return true;
            }

            throw new CalculationException(self::MATRIX_DIMENSION_EXCEPTION);
        }

        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
    }

    //    function checkMatrixDimensions()

    /**
     * set.
     *
     * Set the i,j-th element of the matrix.
     *
     * @param int $i Row position
     * @param int $j Column position
     * @param float|int $c value
     */
    public function set($i = null, $j = null, $c = null): void
    {
        // Optimized set version just has this
        $this->A[$i][$j] = $c;
    }

    //    function set()

    /**
     * identity.
     *
     * Generate an identity matrix.
     *
     * @param int $m Row dimension
     * @param int $n Column dimension
     *
     * @return Matrix Identity matrix
     */
    public function identity($m = null, $n = null)
    {
        return $this->diagonal($m, $n, 1);
    }

    /**
     * diagonal.
     *
     *    Generate a diagonal matrix
     *
     * @param int $m Row dimension
     * @param int $n Column dimension
     * @param mixed $c Diagonal value
     *
     * @return Matrix Diagonal matrix
     */
    public function diagonal($m = null, $n = null, $c = 1)
    {
        $R = new self($m, $n);
        for ($i = 0; $i < $m; ++$i) {
            $R->set($i, $i, $c);
        }

        return $R;
    }

    /**
     * getMatrixByRow.
     *
     *    Get a submatrix by row index/range
     *
     * @param int $i0 Initial row index
     * @param int $iF Final row index
     *
     * @return Matrix Submatrix
     */
    public function getMatrixByRow($i0 = null, $iF = null)
    {
        if (is_int($i0)) {
            if (is_int($iF)) {
                return $this->getMatrix($i0, 0, $iF + 1, $this->n);
            }

            return $this->getMatrix($i0, 0, $i0 + 1, $this->n);
        }

        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
    }

    /**
     * getMatrixByCol.
     *
     *    Get a submatrix by column index/range
     *
     * @param int $j0 Initial column index
     * @param int $jF Final column index
     *
     * @return Matrix Submatrix
     */
    public function getMatrixByCol($j0 = null, $jF = null)
    {
        if (is_int($j0)) {
            if (is_int($jF)) {
                return $this->getMatrix(0, $j0, $this->m, $jF + 1);
            }

            return $this->getMatrix(0, $j0, $this->m, $j0 + 1);
        }

        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
    }

    /**
     * transpose.
     *
     *    Tranpose matrix
     *
     * @return Matrix Transposed matrix
     */
    public function transpose()
    {
        $R = new self($this->n, $this->m);
        for ($i = 0; $i < $this->m; ++$i) {
            for ($j = 0; $j < $this->n; ++$j) {
                $R->set($j, $i, $this->A[$i][$j]);
            }
        }

        return $R;
    }

    //    function transpose()

    /**
     * trace.
     *
     *    Sum of diagonal elements
     *
     * @return float Sum of diagonal elements
     */
    public function trace()
    {
        $s = 0;
        $n = min($this->m, $this->n);
        for ($i = 0; $i < $n; ++$i) {
            $s += $this->A[$i][$i];
        }

        return $s;
    }

    /**
     * plus.
     *
     *    A + B
     *
     * @return Matrix Sum
     */
    public function plus(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $M->set($i, $j, $M->get($i, $j) + $this->A[$i][$j]);
                }
            }

            return $M;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * plusEquals.
     *
     *    A = A + B
     *
     * @return $this
     */
    public function plusEquals(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $validValues = true;
                    $value = $M->get($i, $j);
                    if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
                        $this->A[$i][$j] = trim($this->A[$i][$j], '"');
                        $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
                    }
                    if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
                        $value = trim($value, '"');
                        $validValues &= StringHelper::convertToNumberIfFraction($value);
                    }
                    if ($validValues) {
                        $this->A[$i][$j] += $value;
                    } else {
                        $this->A[$i][$j] = Functions::NAN();
                    }
                }
            }

            return $this;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * minus.
     *
     *    A - B
     *
     * @return Matrix Sum
     */
    public function minus(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $M->set($i, $j, $M->get($i, $j) - $this->A[$i][$j]);
                }
            }

            return $M;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * minusEquals.
     *
     *    A = A - B
     *
     * @return $this
     */
    public function minusEquals(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $validValues = true;
                    $value = $M->get($i, $j);
                    if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
                        $this->A[$i][$j] = trim($this->A[$i][$j], '"');
                        $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
                    }
                    if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
                        $value = trim($value, '"');
                        $validValues &= StringHelper::convertToNumberIfFraction($value);
                    }
                    if ($validValues) {
                        $this->A[$i][$j] -= $value;
                    } else {
                        $this->A[$i][$j] = Functions::NAN();
                    }
                }
            }

            return $this;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * arrayTimes.
     *
     *    Element-by-element multiplication
     *    Cij = Aij * Bij
     *
     * @return Matrix Matrix Cij
     */
    public function arrayTimes(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $M->set($i, $j, $M->get($i, $j) * $this->A[$i][$j]);
                }
            }

            return $M;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * arrayTimesEquals.
     *
     *    Element-by-element multiplication
     *    Aij = Aij * Bij
     *
     * @return $this
     */
    public function arrayTimesEquals(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $validValues = true;
                    $value = $M->get($i, $j);
                    if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
                        $this->A[$i][$j] = trim($this->A[$i][$j], '"');
                        $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
                    }
                    if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
                        $value = trim($value, '"');
                        $validValues &= StringHelper::convertToNumberIfFraction($value);
                    }
                    if ($validValues) {
                        $this->A[$i][$j] *= $value;
                    } else {
                        $this->A[$i][$j] = Functions::NAN();
                    }
                }
            }

            return $this;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * arrayRightDivide.
     *
     *    Element-by-element right division
     *    A / B
     *
     * @return Matrix Division result
     */
    public function arrayRightDivide(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $validValues = true;
                    $value = $M->get($i, $j);
                    if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
                        $this->A[$i][$j] = trim($this->A[$i][$j], '"');
                        $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
                    }
                    if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
                        $value = trim($value, '"');
                        $validValues &= StringHelper::convertToNumberIfFraction($value);
                    }
                    if ($validValues) {
                        if ($value == 0) {
                            //    Trap for Divide by Zero error
                            $M->set($i, $j, '#DIV/0!');
                        } else {
                            $M->set($i, $j, $this->A[$i][$j] / $value);
                        }
                    } else {
                        $M->set($i, $j, Functions::NAN());
                    }
                }
            }

            return $M;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * arrayRightDivideEquals.
     *
     *    Element-by-element right division
     *    Aij = Aij / Bij
     *
     * @return Matrix Matrix Aij
     */
    public function arrayRightDivideEquals(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $this->A[$i][$j] = $this->A[$i][$j] / $M->get($i, $j);
                }
            }

            return $M;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * arrayLeftDivide.
     *
     *    Element-by-element Left division
     *    A / B
     *
     * @return Matrix Division result
     */
    public function arrayLeftDivide(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $M->set($i, $j, $M->get($i, $j) / $this->A[$i][$j]);
                }
            }

            return $M;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * arrayLeftDivideEquals.
     *
     *    Element-by-element Left division
     *    Aij = Aij / Bij
     *
     * @return Matrix Matrix Aij
     */
    public function arrayLeftDivideEquals(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $this->A[$i][$j] = $M->get($i, $j) / $this->A[$i][$j];
                }
            }

            return $M;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * times.
     *
     *    Matrix multiplication
     *
     * @return Matrix Product
     */
    public function times(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $B = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }
                    if ($this->n == $B->m) {
                        $C = new self($this->m, $B->n);
                        for ($j = 0; $j < $B->n; ++$j) {
                            $Bcolj = [];
                            for ($k = 0; $k < $this->n; ++$k) {
                                $Bcolj[$k] = $B->A[$k][$j];
                            }
                            for ($i = 0; $i < $this->m; ++$i) {
                                $Arowi = $this->A[$i];
                                $s = 0;
                                for ($k = 0; $k < $this->n; ++$k) {
                                    $s += $Arowi[$k] * $Bcolj[$k];
                                }
                                $C->A[$i][$j] = $s;
                            }
                        }

                        return $C;
                    }

                    throw new CalculationException(self::MATRIX_DIMENSION_EXCEPTION);
                case 'array':
                    $B = new self($args[0]);
                    if ($this->n == $B->m) {
                        $C = new self($this->m, $B->n);
                        for ($i = 0; $i < $C->m; ++$i) {
                            for ($j = 0; $j < $C->n; ++$j) {
                                $s = '0';
                                for ($k = 0; $k < $C->n; ++$k) {
                                    $s += $this->A[$i][$k] * $B->A[$k][$j];
                                }
                                $C->A[$i][$j] = $s;
                            }
                        }

                        return $C;
                    }

                    throw new CalculationException(self::MATRIX_DIMENSION_EXCEPTION);
                case 'integer':
                    $C = new self($this->A);
                    for ($i = 0; $i < $C->m; ++$i) {
                        for ($j = 0; $j < $C->n; ++$j) {
                            $C->A[$i][$j] *= $args[0];
                        }
                    }

                    return $C;
                case 'double':
                    $C = new self($this->m, $this->n);
                    for ($i = 0; $i < $C->m; ++$i) {
                        for ($j = 0; $j < $C->n; ++$j) {
                            $C->A[$i][$j] = $args[0] * $this->A[$i][$j];
                        }
                    }

                    return $C;
                case 'float':
                    $C = new self($this->A);
                    for ($i = 0; $i < $C->m; ++$i) {
                        for ($j = 0; $j < $C->n; ++$j) {
                            $C->A[$i][$j] *= $args[0];
                        }
                    }

                    return $C;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
            }
        } else {
            throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
        }
    }

    /**
     * power.
     *
     *    A = A ^ B
     *
     * @return $this
     */
    public function power(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $validValues = true;
                    $value = $M->get($i, $j);
                    if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
                        $this->A[$i][$j] = trim($this->A[$i][$j], '"');
                        $validValues &= StringHelper::convertToNumberIfFraction($this->A[$i][$j]);
                    }
                    if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
                        $value = trim($value, '"');
                        $validValues &= StringHelper::convertToNumberIfFraction($value);
                    }
                    if ($validValues) {
                        $this->A[$i][$j] = $this->A[$i][$j] ** $value;
                    } else {
                        $this->A[$i][$j] = Functions::NAN();
                    }
                }
            }

            return $this;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * concat.
     *
     *    A = A & B
     *
     * @return $this
     */
    public function concat(...$args)
    {
        if (count($args) > 0) {
            $match = implode(',', array_map('gettype', $args));

            switch ($match) {
                case 'object':
                    if ($args[0] instanceof self) {
                        $M = $args[0];
                    } else {
                        throw new CalculationException(self::ARGUMENT_TYPE_EXCEPTION);
                    }

                    break;
                case 'array':
                    $M = new self($args[0]);

                    break;
                default:
                    throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);

                    break;
            }
            $this->checkMatrixDimensions($M);
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    $this->A[$i][$j] = trim($this->A[$i][$j], '"') . trim($M->get($i, $j), '"');
                }
            }

            return $this;
        }

        throw new CalculationException(self::POLYMORPHIC_ARGUMENT_EXCEPTION);
    }

    /**
     * Solve A*X = B.
     *
     * @param Matrix $B Right hand side
     *
     * @return Matrix ... Solution if A is square, least squares solution otherwise
     */
    public function solve(self $B)
    {
        if ($this->m == $this->n) {
            $LU = new LUDecomposition($this);

            return $LU->solve($B);
        }
        $QR = new QRDecomposition($this);

        return $QR->solve($B);
    }

    /**
     * Matrix inverse or pseudoinverse.
     *
     * @return Matrix ... Inverse(A) if A is square, pseudoinverse otherwise.
     */
    public function inverse()
    {
        return $this->solve($this->identity($this->m, $this->m));
    }

    /**
     * det.
     *
     *    Calculate determinant
     *
     * @return float Determinant
     */
    public function det()
    {
        $L = new LUDecomposition($this);

        return $L->det();
    }
}
PK     Õ,[|×MbH  bH    SingularValueDecomposition.phpnu „[µü¤        <?php

namespace PhpOffice\PhpSpreadsheet\Shared\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
 *
 *    @version 1.1
 */
class SingularValueDecomposition
{
    /**
     * Internal storage of U.
     *
     * @var array
     */
    private $U = [];

    /**
     * Internal storage of V.
     *
     * @var array
     */
    private $V = [];

    /**
     * Internal storage of singular values.
     *
     * @var array
     */
    private $s = [];

    /**
     * Row dimension.
     *
     * @var int
     */
    private $m;

    /**
     * Column dimension.
     *
     * @var int
     */
    private $n;

    /**
     * Construct the singular value decomposition.
     *
     * Derived from LINPACK code.
     *
     * @param mixed $Arg Rectangular matrix
     */
    public function __construct($Arg)
    {
        // Initialize.
        $A = $Arg->getArray();
        $this->m = $Arg->getRowDimension();
        $this->n = $Arg->getColumnDimension();
        $nu = min($this->m, $this->n);
        $e = [];
        $work = [];
        $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.
        $kMax = max($nct, $nrt);
        for ($k = 0; $k < $kMax; ++$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 && ($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) && ($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) && ($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 = 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;
                } elseif ($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 && ($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 && ($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
    }

    /**
     * Return the left singular vectors.
     *
     * @return Matrix U
     */
    public function getU()
    {
        return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
    }

    /**
     * Return the right singular vectors.
     *
     * @return Matrix V
     */
    public function getV()
    {
        return new Matrix($this->V);
    }

    /**
     * Return the one-dimensional array of singular values.
     *
     * @return array diagonal of S
     */
    public function getSingularValues()
    {
        return $this->s;
    }

    /**
     * Return the diagonal matrix of singular values.
     *
     * @return Matrix S
     */
    public function getS()
    {
        $S = [];
        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.
     *
     * @return float max(S)
     */
    public function norm2()
    {
        return $this->s[0];
    }

    /**
     * Two norm condition number.
     *
     * @return float max(S)/min(S)
     */
    public function cond()
    {
        return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
    }

    /**
     * Effective numerical matrix rank.
     *
     * @return int Number of nonnegligible singular values
     */
    public function rank()
    {
        $eps = 2.0 ** (-52.0);
        $tol = max($this->m, $this->n) * $this->s[0] * $eps;
        $r = 0;
        $iMax = count($this->s);
        for ($i = 0; $i < $iMax; ++$i) {
            if ($this->s[$i] > $tol) {
                ++$r;
            }
        }

        return $r;
    }
}
PK       Õ,[³ˆ\µï  ï                  utils/Maths.phpnu „[µü¤        PK       Õ,[ùU/D  D              .  LUDecomposition.phpnu „[µü¤        PK       Õ,[^99  9              µ   QRDecomposition.phpnu „[µü¤        PK       Õ,[Møô  ô              1<  CHANGELOG.TXTnu „[µü¤        PK       Õ,[úsR<“  “              b>  EigenvalueDecomposition.phpnu „[µü¤        PK       Õ,[|Úp                @¾  CholeskyDecomposition.phpnu „[µü¤        PK       Õ,[ÙRïŒŽ  ŒŽ  
            ŸÍ  Matrix.phpnu „[µü¤        PK       Õ,[|×MbH  bH              e\ SingularValueDecomposition.phpnu „[µü¤        PK      ž  ¥   