public function solve(Vector $b) : ColumnVector
    {
        $k = $this->m - 1;

        [$l, $u, $p] = $this->lu();

        $pb = $p->dot($b);

        $y = [$pb[0] / ($l[0][0] ?: self::EPSILON)];

        for ($i = 1; $i < $this->m; $i++) {
            $sigma = 0;

            for ($j = 0; $j <= $i - 1; $j++) {
                $sigma += $l[$i][$j] * $y[$j];
            }

            $y[] = ($pb[$i] - $sigma) / $l[$i][$i];
        }

        $x = [];

        $x = [$k => $y[$k] / ($l[$k][$k] ?: self::EPSILON)];

        for ($i = $this->m - 2; $i >= 0; $i--) {
            $sigma = 0;

            for ($j = $i + 1; $j < $this->m; $j++) {
                $sigma += $u[$i][$j] * $x[$j];
            }

            $x[$i] = ($y[$i] - $sigma) / $u[$i][$i];
        }

        return ColumnVector::quick(array_reverse($x));
    }