Norm

        public T Norm(double p) {
            if (M == 1 && N == 1) {
                if (p < 1) throw new DomainException();
                else return this[1, 1];
            }
            else if (M == 1 || N == 1) return VectorNorm(p);
            else return MatrixNorm(p);
        }

        private T VectorNorm(double p) {
            // shouldn't be  thrown, but either way this is only for vectors. use matrixnorm
            if (M != 1 && N != 1) throw new InvalidDimException();

            T ans = (dynamic)0;

            if (p >= 1 && !double.IsInfinity(p)) {
                if (N == 1)
                    for (int i = 1; i <= M; i++)
                        ans += (dynamic)Pow(Abs((dynamic)this[i, 1]), p);
                if (M == 1)
                    for (int j = 1; j <= N; j++)
                        ans += (dynamic)Pow(Abs((dynamic)this[1, j]), p);
                return Pow((dynamic)ans, 1 / p);
            } // p = positive int
            else if (p == double.PositiveInfinity) {
                T[] terms = new T[(int)max(M, N)];

                if (N == 1)
                    for (int i = 1; i <= M; i++)
                        terms[i - 1] = Abs((dynamic)this[i, 1]);
                if (M == 1)
                    for (int j = 1; j <= N; j++)
                        terms[j - 1] = Abs((dynamic)this[1, j]);
                return max((dynamic)terms);
            }
            else if (p == double.NegativeInfinity) {
                T[] terms = new T[(int)max(M, N)];

                if (N == 1)
                    for (int i = 1; i <= M; i++)
                        terms[i - 1] = Abs((dynamic)this[i, 1]);
                if (M == 1)
                    for (int j = 1; j <= N; j++)
                        terms[j - 1] = Abs((dynamic)this[1, j]);
                return min((dynamic)terms);
            }
            else throw new DomainException(); // p is not +/-inf or a positive number
        }

        private T MatrixNorm(double p) {
            if (p == 1) {
                T[] norms = new T[N];
                for (int j = 1; j <= N; j++)
                    norms[j - 1] = GetColumnVector(j).VectorNorm(1);
                return max((dynamic)norms);
            }
            else if (p == 2) {
                Matrix<T> A2 = this^2;
                T[] eiv = Eigenvalues();
                return Pow(max((dynamic)eiv), 1/2);
            }
            else if (p == double.PositiveInfinity) {
                T[] norms = new T[M];
                for (int i = 1; i <= M; i++)
                    norms[i - 1] = GetRowVector(i).VectorNorm(1);
                return max((dynamic)norms);
            }
            else throw new DomainException(); // p must be 1, 2, or a +inf
        }