C++程序设计原理与实践习题答案第二十四章第24章习题答案

雨泛回忆殇何潜

已于 2022-04-30 17:06:18 修改

阅读量565

点赞数

分类专栏：《C++程序设计原理与实践》第二版文章标签： c++ 开发语言

于 2022-01-02 20:46:08 首次发布

本文链接：https://blog.csdn.net/qq_36765586/article/details/122280555

版权

《C++程序设计原理与实践》第二版专栏收录该内容

24 篇文章 20 订阅

订阅专栏

第二十四章：数值计算部分习题答案

本章习题所用到的头文件和实现
- Matrix.h
- MatrixIO.h
24.1
24.2
24.4, 24.5, 24.6
24.7
24.9 利用decltype推导模板实参
24.10
24.11
24.12 二维矩阵和一维矩阵的乘法，两个相同维度矩阵的加法

本章习题所用到的头文件和实现

Matrix.h


/*
    warning: this small multidimensional matrix library uses a few features
    not taught in ENGR112 and not explained in elementary textbooks

    (c) Bjarne Stroustrup, Texas A&M University. 

    Use as you like as long as you acknowledge the source.
*/

#ifndef MATRIX_LIB
#define MATRIX_LIB

#include<string>
#include<algorithm>
//#include<iostream>

namespace Numeric_lib {

//-----------------------------------------------------------------------------

struct Matrix_error {
    std::string name;
    Matrix_error(const char* q) :name(q) { }
    Matrix_error(std::string n) :name(n) { }
};

//-----------------------------------------------------------------------------

inline void error(const char* p)
{
    throw Matrix_error(p);
}

//-----------------------------------------------------------------------------

typedef long Index;    // I still dislike unsigned

//-----------------------------------------------------------------------------

// The general Matrix template is simply a prop for its specializations:
template<class T = double, int D = 1> class Matrix {
    // multidimensional matrix class
    // ( ) does multidimensional subscripting
    // [ ] does C style "slicing": gives an N-1 dimensional matrix from an N dimensional one
    // row() is equivalent to [ ]
    // column() is not (yet) implemented because it requires strides.
    // = has copy semantics
    // ( ) and [ ] are range checked
    // slice() to give sub-ranges 
private:
    Matrix();    // this should never be compiled
};

//-----------------------------------------------------------------------------

template<class T = double, int D = 1> class Row ;    // forward declaration

//-----------------------------------------------------------------------------

// function objects for various apply() operations:

template<class T> struct Assign {
    void operator()(T& a, const T& c) { a = c; }
};

template<class T> struct Add_assign {
    void operator()(T& a, const T& c) { a += c; }
};
template<class T> struct Mul_assign {
    void operator()(T& a, const T& c) { a *= c; }
};
template<class T> struct Minus_assign {
    void operator()(T& a, const T& c) { a -= c; }
};
template<class T> struct Div_assign {
    void operator()(T& a, const T& c) { a /= c; }
};
template<class T> struct Mod_assign {
    void operator()(T& a, const T& c) { a %= c; }
};
template<class T> struct Or_assign {
    void operator()(T& a, const T& c) { a |= c; }
};
template<class T> struct Xor_assign {
    void operator()(T& a, const T& c) { a ^= c; }
};
template<class T> struct And_assign {
    void operator()(T& a, const T& c) { a &= c; }
};

template<class T> struct Not_assign {
    void operator()(T& a) { a = !a; }
};

template<class T> struct Not {
    T operator()(T& a) { return !a; }
};

template<class T> struct Unary_minus {
    T operator()(T& a) { return -a; }
};

template<class T> struct Complement {
    T operator()(T& a) { return ~a; }
};

//-----------------------------------------------------------------------------

// Matrix_base represents the common part of the Matrix classes:
template<class T> class Matrix_base {
    // matrixs store their memory (elements) in Matrix_base and have copy semantics
    // Matrix_base does element-wise operations
protected:
    T* elem;    // vector? no: we couldn't easily provide a vector for a slice
    const Index sz;    
    mutable bool owns;
    mutable bool xfer;
public:
    Matrix_base(Index n) :elem(new T[n]()), sz(n), owns(true), xfer(false)
        // matrix of n elements (default initialized)
    {
        // std::cerr << "new[" << n << "]->" << elem << "\n";
    }

    Matrix_base(Index n, T* p) :elem(p), sz(n), owns(false), xfer(false)
        // descriptor for matrix of n elements owned by someone else
    {
    }

    ~Matrix_base()
    {
        if (owns) {
            // std::cerr << "delete[" << sz << "] " << elem << "\n";
            delete[]elem;
        }
    }

    // if necessay, we can get to the raw matrix:
          T* data()       { return elem; }
    const T* data() const { return elem; }
    Index    size() const { return sz; }

    void copy_elements(const Matrix_base& a)
    {
        if (sz!=a.sz) error("copy_elements()");
        for (Index i=0; i<sz; ++i) elem[i] = a.elem[i];
    }

    void base_assign(const Matrix_base& a) { copy_elements(a); }

    void base_copy(const Matrix_base& a)
    {
        if (a.xfer) {          // a is just about to be deleted
                               // so we can transfer ownership rather than copy
            // std::cerr << "xfer @" << a.elem << " [" << a.sz << "]\n";
            elem = a.elem;
            a.xfer = false;    // note: modifies source
            a.owns = false;
        }
        else {
            elem = new T[a.sz];
            // std::cerr << "base copy @" << a.elem << " [" << a.sz << "]\n";
            copy_elements(a);
        }
        owns = true;
        xfer = false;
    }

    // to get the elements of a local matrix out of a function without copying:
    void base_xfer(Matrix_base& x)
    {
        if (owns==false) error("cannot xfer() non-owner");
        owns = false;     // now the elements are safe from deletion by original owner
        x.xfer = true;    // target asserts temporary ownership
        x.owns = true;
    }

    template<class F> void base_apply(F f) { for (Index i = 0; i<size(); ++i) f(elem[i]); }
    template<class F> void base_apply(F f, const T& c) { for (Index i = 0; i<size(); ++i) f(elem[i],c); }
private:
    void operator=(const Matrix_base&);    // no ordinary copy of bases
    Matrix_base(const Matrix_base&);
};

//-----------------------------------------------------------------------------

template<class T> class Matrix<T,1> : public Matrix_base<T> {
    const Index d1;

protected:
    // for use by Row:
    Matrix(Index n1, T* p) : Matrix_base<T>(n1,p), d1(n1)
    {
        // std::cerr << "construct 1D Matrix from data\n";
    }

public:

    Matrix(Index n1) : Matrix_base<T>(n1), d1(n1) { }

    Matrix(Row<T,1>& a) : Matrix_base<T>(a.dim1(),a.p), d1(a.dim1()) 
    { 
        // std::cerr << "construct 1D Matrix from Row\n";
    }

    // copy constructor: let the base do the copy:
    Matrix(const Matrix& a) : Matrix_base<T>(a.size(),0), d1(a.d1)
    {
        // std::cerr << "copy ctor\n";
        this->base_copy(a);
    }

    template<int n> 
    Matrix(const T (&a)[n]) : Matrix_base<T>(n), d1(n)
        // deduce "n" (and "T"), Matrix_base allocates T[n]
    {
        // std::cerr << "matrix ctor\n";
        for (Index i = 0; i<n; ++i) this->elem[i]=a[i];
    }

    Matrix(const T* p, Index n) : Matrix_base<T>(n), d1(n)
        // Matrix_base allocates T[n]
    {
        // std::cerr << "matrix ctor\n";
        for (Index i = 0; i<n; ++i) this->elem[i]=p[i];
    }

    template<class F> Matrix(const Matrix& a, F f) : Matrix_base<T>(a.size()), d1(a.d1)
        // construct a new Matrix with element's that are functions of a's elements:
        // does not modify a unless f has been specifically programmed to modify its argument
        // T f(const T&) would be a typical type for f
    {
        for (Index i = 0; i<this->sz; ++i) this->elem[i] = f(a.elem[i]); 
    }

    template<class F, class Arg> Matrix(const Matrix& a, F f, const Arg& t1) : Matrix_base<T>(a.size()), d1(a.d1)
        // construct a new Matrix with element's that are functions of a's elements:
        // does not modify a unless f has been specifically programmed to modify its argument
        // T f(const T&, const Arg&) would be a typical type for f
    {
        for (Index i = 0; i<this->sz; ++i) this->elem[i] = f(a.elem[i],t1); 
    }

    Matrix& operator=(const Matrix& a)
        // copy assignment: let the base do the copy
    {
        // std::cerr << "copy assignment (" << this->size() << ',' << a.size()<< ")\n";
        if (d1!=a.d1) error("length error in 1D=");
        this->base_assign(a);
        return *this;
    }

    ~Matrix() { }

    Index dim1() const { return d1; }    // number of elements in a row

    Matrix xfer()    // make an Matrix to move elements out of a scope
    {
        Matrix x(dim1(),this->data()); // make a descriptor
        this->base_xfer(x);                  // transfer (temporary) ownership to x
        return x;
    }

    void range_check(Index n1) const
    {
        // std::cerr << "range check: (" << d1 << "): " << n1 << "\n"; 
        if (n1<0 || d1<=n1) error("1D range error: dimension 1");
    }

    // subscripting:
          T& operator()(Index n1)       { range_check(n1); return this->elem[n1]; }
    const T& operator()(Index n1) const { range_check(n1); return this->elem[n1]; }

    // slicing (the same as subscripting for 1D matrixs):
          T& operator[](Index n)       { return row(n); }
    const T& operator[](Index n) const { return row(n); }

          T& row(Index n)       { range_check(n); return this->elem[n]; }
    const T& row(Index n) const { range_check(n); return this->elem[n]; }

    Row<T,1> slice(Index n)
        // the last elements from a[n] onwards
    {
        if (n<0) n=0;
        else if(d1<n) n=d1;// one beyond the end
        return Row<T,1>(d1-n,this->elem+n);
    }

    const Row<T,1> slice(Index n) const
        // the last elements from a[n] onwards
    {
        if (n<0) n=0;
        else if(d1<n) n=d1;// one beyond the end
        return Row<T,1>(d1-n,this->elem+n);
    }

    Row<T,1> slice(Index n, Index m)
        // m elements starting with a[n]
    {
        if (n<0) n=0;
        else if(d1<n) n=d1;    // one beyond the end
        if (m<0) m = 0;
        else if (d1<n+m) m=d1-n;
        return Row<T,1>(m,this->elem+n);
    }

    const Row<T,1> slice(Index n, Index m) const
        // m elements starting with a[n]
    {
        if (n<0) n=0;
        else if(d1<n) n=d1;    // one beyond the end
        if (m<0) m = 0;
        else if (d1<n+m) m=d1-n;
        return Row<T,1>(m,this->elem+n);
    }

    // element-wise operations:
    template<class F> Matrix& apply(F f) { this->base_apply(f); return *this; }
    template<class F> Matrix& apply(F f,const T& c) { this->base_apply(f,c); return *this; }

    Matrix& operator=(const T& c)  { this->base_apply(Assign<T>(),c);       return *this; }

    Matrix& operator*=(const T& c) { this->base_apply(Mul_assign<T>(),c);   return *this; }
    Matrix& operator/=(const T& c) { this->base_apply(Div_assign<T>(),c);   return *this; }
    Matrix& operator%=(const T& c) { this->base_apply(Mod_assign<T>(),c);   return *this; }
    Matrix& operator+=(const T& c) { this->base_apply(Add_assign<T>(),c);   return *this; }
    Matrix& operator-=(const T& c) { this->base_apply(Minus_assign<T>(),c); return *this; }

    Matrix& operator&=(const T& c) { this->base_apply(And_assign<T>(),c);   return *this; }
    Matrix& operator|=(const T& c) { this->base_apply(Or_assign<T>(),c);    return *this; }
    Matrix& operator^=(const T& c) { this->base_apply(Xor_assign<T>(),c);   return *this; }

    Matrix operator!() { return xfer(Matrix(*this,Not<T>())); }
    Matrix operator-() { return xfer(Matrix(*this,Unary_minus<T>())); }
    Matrix operator~() { return xfer(Matrix(*this,Complement<T>()));  }

    template<class F> Matrix apply_new(F f) { return xfer(Matrix(*this,f)); }
    
    void swap_rows(Index i, Index j)
        // swap_rows() uses a row's worth of memory for better run-time performance
        // if you want pairwise swap, just write it yourself
    {
        if (i == j) return;
    /*
        Matrix<T,1> temp = (*this)[i];
        (*this)[i] = (*this)[j];
        (*this)[j] = temp;
    */
        Index max = (*this)[i].size();
        for (Index ii=0; ii<max; ++ii) std::swap((*this)(i,ii),(*this)(j,ii));
    }

    
    void swap_column(Index i, Index j)
    {
        if (i == j) return;
        //一维矩阵（一维向量）只有一行
        std::swap((*this)(i), (*this)(j));
    }
};

//-----------------------------------------------------------------------------

template<class T> class Matrix<T,2> : public Matrix_base<T> {
    const Index d1;
    const Index d2;

protected:
    // for use by Row:
    Matrix(Index n1, Index n2, T* p) : Matrix_base<T>(n1*n2,p), d1(n1), d2(n2) 
    {
        // std::cerr << "construct 3D Matrix from data\n";
    }

public:

    Matrix(Index n1, Index n2) : Matrix_base<T>(n1*n2), d1(n1), d2(n2) { }

    Matrix(Row<T,2>& a) : Matrix_base<T>(a.dim1()*a.dim2(),a.p), d1(a.dim1()), d2(a.dim2())
    { 
        // std::cerr << "construct 2D Matrix from Row\n";
    }

    // copy constructor: let the base do the copy:
    Matrix(const Matrix& a) : Matrix_base<T>(a.size(),0), d1(a.d1), d2(a.d2)
    {
        // std::cerr << "copy ctor\n";
        this->base_copy(a);
    }

    template<int n1, int n2> 
    Matrix(const T (&a)[n1][n2]) : Matrix_base<T>(n1*n2), d1(n1), d2(n2)
        // deduce "n1", "n2" (and "T"), Matrix_base allocates T[n1*n2]
    {
        // std::cerr << "matrix ctor\n";
        for (Index i = 0; i<n1; ++i)
            for (Index j = 0; j<n2; ++j) this->elem[i*n2+j]=a[i][j];
    }

    template<class F> Matrix(const Matrix& a, F f) : Matrix_base<T>(a.size()), d1(a.d1), d2(a.d2)
        // construct a new Matrix with element's that are functions of a's elements:
        // does not modify a unless f has been specifically programmed to modify its argument
        // T f(const T&) would be a typical type for f
    {
        for (Index i = 0; i<this->sz; ++i) this->elem[i] = f(a.elem[i]); 
    }

    template<class F, class Arg> Matrix(const Matrix& a, F f, const Arg& t1) : Matrix_base<T>(a.size()), d1(a.d1), d2(a.d2)
        // construct a new Matrix with element's that are functions of a's elements:
        // does not modify a unless f has been specifically programmed to modify its argument
        // T f(const T&, const Arg&) would be a typical type for f
    {
        for (Index i = 0; i<this->sz; ++i) this->elem[i] = f(a.elem[i],t1); 
    }

    Matrix& operator=(const Matrix& a)
        // copy assignment: let the base do the copy
    {
        // std::cerr << "copy assignment (" << this->size() << ',' << a.size()<< ")\n";
        if (d1!=a.d1 || d2!=a.d2) error("length error in 2D =");
        this->base_assign(a);
        return *this;
    }

    ~Matrix() { }
    
    Index dim1() const { return d1; }    // number of elements in a row
    Index dim2() const { return d2; }    // number of elements in a column

    Matrix xfer()    // make an Matrix to move elements out of a scope
    {
        Matrix x(dim1(),dim2(),this->data()); // make a descriptor
        this->base_xfer(x);            // transfer (temporary) ownership to x
        return x;
    }

    void range_check(Index n1, Index n2) const
    {
        // std::cerr << "range check: (" << d1 << "," << d2 << "): " << n1 << " " << n2 << "\n";
        if (n1<0 || d1<=n1) error("2D range error: dimension 1");
        if (n2<0 || d2<=n2) error("2D range error: dimension 2");
    }

    // subscripting:
          T& operator()(Index n1, Index n2)       { range_check(n1,n2); return this->elem[n1*d2+n2]; }
    const T& operator()(Index n1, Index n2) const { range_check(n1,n2); return this->elem[n1*d2+n2]; }

    // slicing (return a row):
          Row<T,1> operator[](Index n)       { return row(n); }
    const Row<T,1> operator[](Index n) const { return row(n); }

          Row<T,1> row(Index n)       { range_check(n,0); return Row<T,1>(d2,&this->elem[n*d2]); }
    const Row<T,1> row(Index n) const { range_check(n,0); return Row<T,1>(d2,&this->elem[n*d2]); }

    Row<T,2> slice(Index n)
        // rows [n:d1)
    {
        if (n<0) n=0;
        else if(d1<n) n=d1;    // one beyond the end
        return Row<T,2>(d1-n,d2,this->elem+n*d2);
    }

    const Row<T,2> slice(Index n) const
        // rows [n:d1)
    {
        if (n<0) n=0;
        else if(d1<n) n=d1;    // one beyond the end
        return Row<T,2>(d1-n,d2,this->elem+n*d2);
    }

    Row<T,2> slice(Index n, Index m)
        // the rows [n:m)
    {
        if (n<0) n=0;
        if(d1<m) m=d1;    // one beyond the end
        return Row<T,2>(m-n,d2,this->elem+n*d2);

    }

    const Row<T,2> slice(Index n, Index m) const
        // the rows [n:sz)
    {
        if (n<0) n=0;
        if(d1<m) m=d1;    // one beyond the end
        return Row<T,2>(m-n,d2,this->elem+n*d2);
    }

    // Column<T,1> column(Index n); // not (yet) implemented: requies strides and operations on columns

    // element-wise operations:
    template<class F> Matrix& apply(F f)            { this->base_apply(f);   return *this; }
    template<class F> Matrix& apply(F f,const T& c) { this->base_apply(f,c); return *this; }

    Matrix& operator=(const T& c)  { this->base_apply(Assign<T>(),c);       return *this; }

    Matrix& operator*=(const T& c) { this->base_apply(Mul_assign<T>(),c);   return *this; }
    Matrix& operator/=(const T& c) { this->base_apply(Div_assign<T>(),c);   return *this; }
    Matrix& operator%=(const T& c) { this->base_apply(Mod_assign<T>(),c);   return *this; }
    Matrix& operator+=(const T& c) { this->base_apply(Add_assign<T>(),c);   return *this; }
    Matrix& operator-=(const T& c) { this->base_apply(Minus_assign<T>(),c); return *this; }

    Matrix& operator&=(const T& c) { this->base_apply(And_assign<T>(),c);   return *this; }
    Matrix& operator|=(const T& c) { this->base_apply(Or_assign<T>(),c);    return *this; }
    Matrix& operator^=(const T& c) { this->base_apply(Xor_assign<T>(),c);   return *this; }

    Matrix operator!() { return xfer(Matrix(*this,Not<T>())); }
    Matrix operator-() { return xfer(Matrix(*this,Unary_minus<T>())); }
    Matrix operator~() { return xfer(Matrix(*this,Complement<T>()));  }

    template<class F> Matrix apply_new(F f) { return xfer(Matrix(*this,f)); }
    
    void swap_rows(Index i, Index j)
        // swap_rows() uses a row's worth of memory for better run-time performance
        // if you want pairwise swap, just write it yourself
    {
        if (i == j) return;
    /*
        Matrix<T,1> temp = (*this)[i];
        (*this)[i] = (*this)[j];
        (*this)[j] = temp;
    */
        Index max = (*this)[i].size();
        for (Index ii=0; ii<max; ++ii) std::swap((*this)(i,ii),(*this)(j,ii));
    }

    void swap_column(Index i, Index j)
    {
        if (i == j) return;

        Index r_max = dim1();
        for (Index ri = 0; ri < r_max; ++ri)
            std::swap((*this)(ri, i), (*this)(ri, j));
    }
};

//-----------------------------------------------------------------------------

template<class T> class Matrix<T,3> : public Matrix_base<T> {
    const Index d1;
    const Index d2;
    const Index d3;

protected:
    // for use by Row:
    Matrix(Index n1, Index n2, Index n3, T* p) : Matrix_base<T>(n1*n2*n3,p), d1(n1), d2(n2), d3(n3) 
    {
        // std::cerr << "construct 3D Matrix from data\n";
    }

public:

    Matrix(Index n1, Index n2, Index n3) : Matrix_base<T>(n1*n2*n3), d1(n1), d2(n2), d3(n3) { }

    Matrix(Row<T,3>& a) : Matrix_base<T>(a.dim1()*a.dim2()*a.dim3(),a.p), d1(a.dim1()), d2(a.dim2()), d3(a.dim3())
    { 
        // std::cerr << "construct 3D Matrix from Row\n";
    }

    // copy constructor: let the base do the copy:
    Matrix(const Matrix& a) : Matrix_base<T>(a.size(),0), d1(a.d1), d2(a.d2), d3(a.d3)
    {
        // std::cerr << "copy ctor\n";
        this->base_copy(a);
    }

    template<int n1, int n2, int n3> 
    Matrix(const T (&a)[n1][n2][n3]) : Matrix_base<T>(n1*n2), d1(n1), d2(n2), d3(n3)
        // deduce "n1", "n2", "n3" (and "T"), Matrix_base allocates T[n1*n2*n3]
    {
        // std::cerr << "matrix ctor\n";
        for (Index i = 0; i<n1; ++i)
            for (Index j = 0; j<n2; ++j)
                for (Index k = 0; k<n3; ++k)
                    this->elem[i*n2*n3+j*n3+k]=a[i][j][k];
    }

    template<class F> Matrix(const Matrix& a, F f) : Matrix_base<T>(a.size()), d1(a.d1), d2(a.d2), d3(a.d3)
        // construct a new Matrix with element's that are functions of a's elements:
        // does not modify a unless f has been specifically programmed to modify its argument
        // T f(const T&) would be a typical type for f
    {
        for (Index i = 0; i<this->sz; ++i) this->elem[i] = f(a.elem[i]); 
    }

    template<class F, class Arg> Matrix(const Matrix& a, F f, const Arg& t1) : Matrix_base<T>(a.size()), d1(a.d1), d2(a.d2), d3(a.d3)
        // construct a new Matrix with element's that are functions of a's elements:
        // does not modify a unless f has been specifically programmed to modify its argument
        // T f(const T&, const Arg&) would be a typical type for f
    {
        for (Index i = 0; i<this->sz; ++i) this->elem[i] = f(a.elem[i],t1); 
    }

    Matrix& operator=(const Matrix& a)
        // copy assignment: let the base do the copy
    {
        // std::cerr << "copy assignment (" << this->size() << ',' << a.size()<< ")\n";
        if (d1!=a.d1 || d2!=a.d2 || d3!=a.d3) error("length error in 2D =");
        this->base_assign(a);
        return *this;
    }

    ~Matrix() { }

    Index dim1() const { return d1; }    // number of elements in a row
    Index dim2() const { return d2; }    // number of elements in a column
    Index dim3() const { return d3; }    // number of elements in a depth

    Matrix xfer()    // make an Matrix to move elements out of a scope
    {
        Matrix x(dim1(),dim2(),dim3(),this->data()); // make a descriptor
        this->base_xfer(x);            // transfer (temporary) ownership to x
        return x;
    }

    void range_check(Index n1, Index n2, Index n3) const
    {
        // std::cerr << "range check: (" << d1 << "," << d2 << "): " << n1 << " " << n2 << "\n";
        if (n1<0 || d1<=n1) error("3D range error: dimension 1");
        if (n2<0 || d2<=n2) error("3D range error: dimension 2");
        if (n3<0 || d3<=n3) error("3D range error: dimension 3");
    }

    // subscripting:
          T& operator()(Index n1, Index n2, Index n3)       { range_check(n1,n2,n3); return this->elem[d2*d3*n1+d3*n2+n3]; }; 
    const T& operator()(Index n1, Index n2, Index n3) const { range_check(n1,n2,n3); return this->elem[d2*d3*n1+d3*n2+n3]; };

    // slicing (return a row):
          Row<T,2> operator[](Index n)       { return row(n); }
    const Row<T,2> operator[](Index n) const { return row(n); }

          Row<T,2> row(Index n)       { range_check(n,0,0); return Row<T,2>(d2,d3,&this->elem[n*d2*d3]); }
    const Row<T,2> row(Index n) const { range_check(n,0,0); return Row<T,2>(d2,d3,&this->elem[n*d2*d3]); }

    Row<T,3> slice(Index n)
        // rows [n:d1)
    {
        if (n<0) n=0;
        else if(d1<n) n=d1;    // one beyond the end
        return Row<T,3>(d1-n,d2,d3,this->elem+n*d2*d3);
    }

    const Row<T,3> slice(Index n) const
        // rows [n:d1)
    {
        if (n<0) n=0;
        else if(d1<n) n=d1;    // one beyond the end
        return Row<T,3>(d1-n,d2,d3,this->elem+n*d2*d3);
    }

    Row<T,3> slice(Index n, Index m)
        // the rows [n:m)
    {
        if (n<0) n=0;
        if(d1<m) m=d1;    // one beyond the end
        return Row<T,3>(m-n,d2,d3,this->elem+n*d2*d3);

    }

    const Row<T,3> slice(Index n, Index m) const
        // the rows [n:sz)
    {
        if (n<0) n=0;
        if(d1<m) m=d1;    // one beyond the end
        return Row<T,3>(m-n,d2,d3,this->elem+n*d2*d3);
    }

    // Column<T,2> column(Index n); // not (yet) implemented: requies strides and operations on columns

    // element-wise operations:
    template<class F> Matrix& apply(F f)            { this->base_apply(f);   return *this; }
    template<class F> Matrix& apply(F f,const T& c) { this->base_apply(f,c); return *this; }

    Matrix& operator=(const T& c)  { this->base_apply(Assign<T>(),c);       return *this; }
                                                                            
    Matrix& operator*=(const T& c) { this->base_apply(Mul_assign<T>(),c);   return *this; }
    Matrix& operator/=(const T& c) { this->base_apply(Div_assign<T>(),c);   return *this; }
    Matrix& operator%=(const T& c) { this->base_apply(Mod_assign<T>(),c);   return *this; }
    Matrix& operator+=(const T& c) { this->base_apply(Add_assign<T>(),c);   return *this; }
    Matrix& operator-=(const T& c) { this->base_apply(Minus_assign<T>(),c); return *this; }

    Matrix& operator&=(const T& c) { this->base_apply(And_assign<T>(),c);   return *this; }
    Matrix& operator|=(const T& c) { this->base_apply(Or_assign<T>(),c);    return *this; }
    Matrix& operator^=(const T& c) { this->base_apply(Xor_assign<T>(),c);   return *this; }

    Matrix operator!() { return xfer(Matrix(*this,Not<T>())); }
    Matrix operator-() { return xfer(Matrix(*this,Unary_minus<T>())); }
    Matrix operator~() { return xfer(Matrix(*this,Complement<T>()));  }

    template<class F> Matrix apply_new(F f) { return xfer(Matrix(*this,f)); }
    
    void swap_rows(Index i, Index j)
        // swap_rows() uses a row's worth of memory for better run-time performance
        // if you want pairwise swap, just write it yourself
    {
        if (i == j) return;
        
        Matrix<T,2> temp = (*this)[i];
        (*this)[i] = (*this)[j];
        (*this)[j] = temp;
    }
};

//-----------------------------------------------------------------------------

template<class T> Matrix<T> scale_and_add(const Matrix<T>& a, T c, const Matrix<T>& b)
    //  Fortran "saxpy()" ("fma" for "fused multiply-add").
    // will the copy constructor be called twice and defeat the xfer optimization?
{
    if (a.size() != b.size()) error("sizes wrong for scale_and_add()");
    Matrix<T> res(a.size());
    for (Index i = 0; i<a.size(); ++i) res[i] += a[i]*c+b[i];
    return res.xfer();
}

//-----------------------------------------------------------------------------

template<class T> T dot_product(const Matrix<T>&a , const Matrix<T>& b)
{
    if (a.size() != b.size()) error("sizes wrong for dot product");
    T sum = 0;
    for (Index i = 0; i<a.size(); ++i) sum += a[i]*b[i];
    return sum;
}

//-----------------------------------------------------------------------------

template<class T, int N> Matrix<T,N> xfer(Matrix<T,N>& a)
{
    return a.xfer();
}

//-----------------------------------------------------------------------------

template<class F, class A>            A apply(F f, A x)        { A res(x,f);   return xfer(res); }
template<class F, class Arg, class A> A apply(F f, A x, Arg a) { A res(x,f,a); return xfer(res); }

//-----------------------------------------------------------------------------

// The default values for T and D have been declared before.
template<class T, int D> class Row {
    // general version exists only to allow specializations
private:
        Row();
};

//-----------------------------------------------------------------------------

template<class T> class Row<T,1> : public Matrix<T,1> {
public:
    Row(Index n, T* p) : Matrix<T,1>(n,p)
    {
    }

    Matrix<T,1>& operator=(const T& c) { this->base_apply(Assign<T>(),c); return *this; }

    Matrix<T,1>& operator=(const Matrix<T,1>& a)
    {
        return *static_cast<Matrix<T,1>*>(this)=a;
    }
};

//-----------------------------------------------------------------------------

template<class T> class Row<T,2> : public Matrix<T,2> {
public:
    Row(Index n1, Index n2, T* p) : Matrix<T,2>(n1,n2,p)
    {
    }
        
    Matrix<T,2>& operator=(const T& c) { this->base_apply(Assign<T>(),c); return *this; }

    Matrix<T,2>& operator=(const Matrix<T,2>& a)
    {
        return *static_cast<Matrix<T,2>*>(this)=a;
    }
};

//-----------------------------------------------------------------------------

template<class T> class Row<T,3> : public Matrix<T,3> {
public:
    Row(Index n1, Index n2, Index n3, T* p) : Matrix<T,3>(n1,n2,n3,p)
    {
    }

    Matrix<T,3>& operator=(const T& c) { this->base_apply(Assign<T>(),c); return *this; }

    Matrix<T,3>& operator=(const Matrix<T,3>& a)
    {
        return *static_cast<Matrix<T,3>*>(this)=a;
    }
};

//-----------------------------------------------------------------------------

template<class T, int N> Matrix<T,N-1> scale_and_add(const Matrix<T,N>& a, const Matrix<T,N-1> c, const Matrix<T,N-1>& b)
{
    Matrix<T> res(a.size());
    if (a.size() != b.size()) error("sizes wrong for scale_and_add");
    for (Index i = 0; i<a.size(); ++i) res[i] += a[i]*c+b[i];
    return res.xfer();
}

//-----------------------------------------------------------------------------

template<class T, int D> Matrix<T,D> operator*(const Matrix<T,D>& m, const T& c) { Matrix<T,D> r(m); return r*=c; }
template<class T, int D> Matrix<T,D> operator/(const Matrix<T,D>& m, const T& c) { Matrix<T,D> r(m); return r/=c; }
template<class T, int D> Matrix<T,D> operator%(const Matrix<T,D>& m, const T& c) { Matrix<T,D> r(m); return r%=c; }
template<class T, int D> Matrix<T,D> operator+(const Matrix<T,D>& m, const T& c) { Matrix<T,D> r(m); return r+=c; }
template<class T, int D> Matrix<T,D> operator-(const Matrix<T,D>& m, const T& c) { Matrix<T,D> r(m); return r-=c; }

template<class T, int D> Matrix<T,D> operator&(const Matrix<T,D>& m, const T& c) { Matrix<T,D> r(m); return r&=c; }
template<class T, int D> Matrix<T,D> operator|(const Matrix<T,D>& m, const T& c) { Matrix<T,D> r(m); return r|=c; }
template<class T, int D> Matrix<T,D> operator^(const Matrix<T,D>& m, const T& c) { Matrix<T,D> r(m); return r^=c; }

//-----------------------------------------------------------------------------

}
#endif

MatrixIO.h

#include <iostream>
#include "Matrix.h"

namespace Numeric_lib {

//-----------------------------------------------------------------------------

template<class T> std::ostream& operator<<(std::ostream& os, const Matrix<T>& v)
{
    os << '{';

    for (Index i = 0; i<v.dim1(); ++i) {
        os << "  ";
        os << v(i);
    }

    os << '}';
    return os;
}

//-----------------------------------------------------------------------------

template<class T> std::ostream& operator<<(std::ostream& os, const Matrix<T,2>& m)
{
    os << "{\n";

    for (Index i = 0; i<m.dim1(); ++i)
        os << m[i] << '\n';

    os << '}';
    return os;
}

//-----------------------------------------------------------------------------

template<class T> std::istream& operator>>(std::istream& is, Matrix<T>& v)
{
    char ch;
    is >> ch;

    if (ch!='{') error("'{' missing in Matrix<T,1> input");

    for (Index i = 0; i<v.dim1(); ++i)
        is >> v(i);

    is >> ch;

    if (ch!='}') error("'}' missing in Matrix<T,1> input");
  
    return is;
}

//-----------------------------------------------------------------------------

template<class T> std::istream& operator>>(std::istream& is, Matrix<T,2>& m)
{
    char ch;
    is >> ch;

    if (ch!='{') error("'{' missing in Matrix<T,2> input");

    for (Index i = 0; i<m.dim1(); ++i)
    {
        Matrix<T,1> tmp(m.dim2());
        is >> tmp;
        m[i] = tmp;
    }

    is >> ch;

    if (ch!='}') error("'}' missing in Matrix<T,2> input");

    return is;
}

//-----------------------------------------------------------------------------
}

24.1

#include <iostream>
#include "../../Matrix.h"
#include "../../MatrixIO.h"

using namespace Numeric_lib;
using namespace std;

//用于 a.apply(f) 的 triple()函数
void triple1(int& elem)
{
	elem *= 3;
}


//用于 apply(f, a) 的 triple()函数
int triple2(const int& elem)
{
	return elem * 3;
}


//同时用于 a.apply(f) 和 apply(f, a) 的 triple()函数
int triple3(int& elem)
{
	elem *= 3;
	return elem;
}

int main()
{
	Matrix<int, 1> mi1({ 1,2,3,4,5 });
	Matrix<int, 1> mi2 = apply(triple2, mi1);

	mi1.apply(triple1);
	cout << mi1 << endl;
	cout << mi2 << endl;

	mi1.apply(triple3);
	cout << mi1 << endl;

	Matrix<int, 1> mi4 = apply(triple3, mi2);
	cout << mi4 << endl;

	return 0;
}

24.2

#include <iostream>
#include "../../Matrix.h"
#include "../../MatrixIO.h"

using namespace Numeric_lib;
using namespace std;

//用于 a.apply(f) 的 triple()函数对象
struct triple1 {
	void operator()(int& elem)
	{
		elem *= 3;
	}
};


//用于 apply(f, a) 的 triple()函数
class triple2 {
public:
	int operator()(const int& elem)
	{	return elem * 3;	}
};


//同时用于 a.apply(f) 和 apply(f, a) 的 triple()函数
class triple3 {
public:
	int operator()(int& elem)
	{
		elem *= 3;
		return elem;
	}
};


int main()
try{
	Matrix<int, 1> mi1({ 1,2,3,4,5 });
	Matrix<int, 1> mi2 = apply(triple2(), mi1);

	mi1.apply(triple1());
	cout << mi1 << endl;
	cout << mi2 << endl;

	mi1.apply(triple3());
	cout << mi1 << endl;

	Matrix<int, 1> mi4 = apply(triple3(), mi2);
	cout << mi4 << endl;

	return 0;
}
catch (std::exception& e) {
	cerr << "Exception: " << e.what() << '\n';
}
catch (...) {
	cerr << "Exception\n";
}

24.4, 24.5, 24.6

//#include "../../std_lib_facilities.h"
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <random>
#include "../../Matrix.h"
#include "../../MatrixIO.h"

using namespace std;
using Numeric_lib::Index;

//------------------------------------------------------------------------------
//声明error相关函数
inline void error(const string& errormessage);
inline void error(string s1, string s2);
inline void error(string s1, int n);
//------------------------------------------------------------------------------


//------------------------------------------------------------------------------

typedef Numeric_lib::Matrix<double, 2> Matrix;
typedef Numeric_lib::Matrix<double, 1> Vector;

//------------------------------------------------------------------------------

void classical_elimination(Matrix& A, Vector& b);

//------------------------------------------------------------------------------

void elim_with_partial_pivot(Matrix& A, Vector& b);

//------------------------------------------------------------------------------

Vector back_substitution(const Matrix& A, const Vector& b);

//------------------------------------------------------------------------------

template<class T>
string to_string(const T& t)
{
    ostringstream os;
    os << t;
    return os.str();
}

//------------------------------------------------------------------------------

// An exception of this type is thrown when elimination fails.
struct Elim_failure : std::domain_error {
    Elim_failure(std::string s) : std::domain_error(s) { }
};

//------------------------------------------------------------------------------

// An exception of this type is thrown when back substitution fails.
struct Back_subst_failure : std::domain_error {
    Back_subst_failure(std::string s) : std::domain_error(s) { }
};

//------------------------------------------------------------------------------

Vector classical_gaussian_elimination(Matrix A, Vector b)   //第一版 经典高斯消去法
{
    classical_elimination(A, b);
    return back_substitution(A, b);
}

//------------------------------------------------------------------------------

Vector pivotal_elimination(Matrix A, Vector b)      //第二版 选取主元的消去法
{
    elim_with_partial_pivot(A, b);
    return back_substitution(A, b);
}

//------------------------------------------------------------------------------

void classical_elimination(Matrix& A, Vector& b)
{
    const Index n = A.dim1();

    // traverse from 1st column to the next-to-last
    // filling zeros into all elements under the diagonal:
    for (Index j = 0; j < n - 1; ++j) {
        const double pivot = A(j, j);
        if (pivot == 0) throw Elim_failure("Elimination failure in row " + to_string(j));

        // fill zeros into each element under the diagonal of the ith row:
        for (Index i = j + 1; i < n; ++i) {
            const double mult = A(i, j) / pivot;
            //A[i].slice(j) = scale_and_add(A[j].slice(j), -mult, A[i].slice(j));
            //将 slice_and_add() 替换为循环，并测试
            for (Index k = j; k < n; ++k)       //k 是列索引，索引A(i,k)
            {
                A(i, k) += -mult * A(j, k);
            }
            b(i) -= mult * b(j); // make the corresponding change to b
        }
    }
}

//------------------------------------------------------------------------------

void elim_with_partial_pivot(Matrix& A, Vector& b)
{
    const Index n = A.dim1();

    for (Index j = 0; j < n; ++j) {
        Index pivot_row = j;

        // look for a suitable pivot:
        for (Index k = j + 1; k < n; ++k)
            if (abs(A(k, j)) > abs(A(pivot_row, j))) pivot_row = k;

        // swap the rows if we found a better pivot:
        if (pivot_row != j) {
            A.swap_rows(j, pivot_row);
            std::swap(b(j), b(pivot_row));
        }

        // elimination:
        for (Index i = j + 1; i < n; ++i) {
            const double pivot = A(j, j);
            if (pivot == 0) error("can't solve: pivot==0");
            const double mult = A(i, j) / pivot;
            //A[i].slice(j) = scale_and_add(A[j].slice(j), -mult, A[i].slice(j));
            //将 slice_and_add() 替换为循环，并测试
            for (Index k = j; k < n; ++k)       //k 是列索引，索引A(i,k)
            {
                A(i, k) += -mult * A(j, k);
            }
            b(i) -= mult * b(j);
        }
    }
}

//------------------------------------------------------------------------------

Vector back_substitution(const Matrix& A, const Vector& b)
{
    const Index n = A.dim1();
    Vector x(n);

    for (Index i = n - 1; i >= 0; --i) {
        //double s = b(i) - dot_product(A[i].slice(i + 1), x.slice(i + 1));
        //将 dot_product() 替换为循环，并测试
        double s = 0.0;
        for (Index j = i + 1; j < A.dim2(); ++j)
            s += A(i, j) * x(j);
        s = b(i) - s;
        if (double m = A(i, i))
            x(i) = s / m;
        else
            //throw Back_subst_failure("Back substitution failure in row " + to_string(i));   //第一版 classical_elimination 下的异常
            error("can't solve: diagonal element is zero during back substitution");        //第二版 elim_with_partial_pivot 下的异常
    }

    return x;
}

//------------------------------------------------------------------------------

Vector random_vector(Index n)
{
    Vector v(n);
    const int max = 50;
    default_random_engine ran{ };                    //生成整数
    uniform_real_distribution<> ureal{ 0,max };     //将 int 映射为 [0:max) 中的 double
    for (Index i = 0; i < n; ++i)
        v(i) = ureal(ran);

    return v;
}

//------------------------------------------------------------------------------

Vector random_vector_with_seed(Index n, unsigned seed)
{
    Vector v(n);
    const int max = 50;
    default_random_engine ran{ seed };                    //生成整数
    uniform_real_distribution<> ureal{ 0,max };     //将 int 映射为 [0:max) 中的 double
    for (Index i = 0; i < n; ++i)
        v(i) = ureal(ran);

    return v;
}

//------------------------------------------------------------------------------

Matrix random_matrix(Index n)
{
    Matrix m(n, n);

    for (Index i = 0; i < n; ++i)
        m[i] = random_vector_with_seed(n, i);

    return m;
}

//------------------------------------------------------------------------------

Vector operator*(const Matrix& m, const Vector& u)
{
    const Index n = m.dim1();
    Vector v(n);
    for (Index i = 0; i < n; ++i) 
        v(i) = dot_product(m[i], u);
    return v;
}

//------------------------------------------------------------------------------

void solve_random_system(Index n)
{
    Matrix A = random_matrix(n);    // see ?4.6
    Vector b = random_vector(n);

    cout << "A = " << A << endl;
    cout << "b = " << b << endl;

    try {
        Vector x = classical_gaussian_elimination(A, b);
        cout << "classical elim solution is x = " << x << endl;
        Vector y = pivotal_elimination(A, b);
        cout << "pivotal elim solution is y = " << y << endl;
        Vector v = A * x;
        cout << " A * x = " << v << endl;
        Vector w = A * y;
        cout << " A * y = " << w << endl;
    }
    catch (const exception& e) {
        cerr << e.what() << std::endl;
    }
}

//------------------------------------------------------------------------------

int main()
try
{
    //double a[3][3] =
    //{
    //   {1,2,3},
    //   {2,3,4},
    //   {3,4,1}
    //};
    //double b[3] = { 14,20,14 };

    //double a[2][2] =
    //{
    //   {0,1},
    //   {1,0}
    //};
    //double b[2] = { 5,6 };

    //Matrix A(a);
    //Vector B(b);

    //cout << "Solving A*x=B" << endl;
    //cout << "A =\n" << A << endl;
    //cout << "B = " << B << endl;

    //Vector x = classical_gaussian_elimination(A, B);
    //Vector y = pivotal_elimination(A, B);
    //cout << "x = " << x << endl;
    //cout << "y = " << y << endl;


    solve_random_system(3);
    solve_random_system(4);
    solve_random_system(5);

    return 0;
}
catch (exception& e) {
    cerr << "error: " << e.what() << '\n';
    return 1;
}
catch (...) {
    cerr << "Oops: unknown exception!\n";
    return 2;
}




inline void error(const string& errormessage)
{
    cerr << errormessage << endl;
    throw runtime_error(errormessage);
}

//------------------------------------------------------------------------------

inline void error(string s1, string s2)
{
    error(s1 + s2);
}

//------------------------------------------------------------------------------

inline void error(string s1, int n)
{
    ostringstream os;
    os << n;
    error(s1, os.str());
}

24.7

#include <iostream>
#include <sstream>
#include <stdexcept>
#include <random>
#include <vector>


using namespace std;


//------------------------------------------------------------------------------
//声明error相关函数
inline void error(const string& errormessage);
inline void error(string s1, string s2);
inline void error(string s1, int n);
//------------------------------------------------------------------------------

using Index = long;

//------------------------------------------------------------------------------

typedef vector<vector<double>> Matrix;
typedef vector<double> Vector;

//------------------------------------------------------------------------------

void classical_elimination(Matrix& A, Vector& b);

//------------------------------------------------------------------------------

void elim_with_partial_pivot(Matrix& A, Vector& b);

//------------------------------------------------------------------------------

Vector back_substitution(const Matrix& A, const Vector& b);

//------------------------------------------------------------------------------

template<class T>
string to_string(const T& t)
{
    ostringstream os;
    os << t;
    return os.str();
}

//------------------------------------------------------------------------------

// An exception of this type is thrown when elimination fails.
struct Elim_failure : std::domain_error {
    Elim_failure(std::string s) : std::domain_error(s) { }
};

//------------------------------------------------------------------------------

// An exception of this type is thrown when back substitution fails.
struct Back_subst_failure : std::domain_error {
    Back_subst_failure(std::string s) : std::domain_error(s) { }
};

//------------------------------------------------------------------------------

Vector classical_gaussian_elimination(Matrix A, Vector b)   //第一版 经典高斯消去法
{
    classical_elimination(A, b);
    return back_substitution(A, b);
}

//------------------------------------------------------------------------------

Vector pivotal_elimination(Matrix A, Vector b)      //第二版 选取主元的消去法
{
    elim_with_partial_pivot(A, b);
    return back_substitution(A, b);
}

//------------------------------------------------------------------------------

void classical_elimination(Matrix& A, Vector& b)
{
    //const Index n = A.dim1();
    const Index n = A.size();

    // traverse from 1st column to the next-to-last
    // filling zeros into all elements under the diagonal:
    for (Index j = 0; j < n - 1; ++j) {
        const double pivot = A[j][j];
        if (pivot == 0) throw Elim_failure("Elimination failure in row " + to_string(j));

        // fill zeros into each element under the diagonal of the ith row:
        for (Index i = j + 1; i < n; ++i) {
            const double mult = A[i][j] / pivot;
            for (Index k = j; k < n; ++k)       //k 是列索引，索引A(i,k)
            {
                A[i][k] += -mult * A[j][k];
            }
            b[i] -= mult * b[j]; // make the corresponding change to b
        }
    }
}

//------------------------------------------------------------------------------

void elim_with_partial_pivot(Matrix& A, Vector& b)
{
    const Index n = A.size();

    for (Index j = 0; j < n; ++j) {
        Index pivot_row = j;

        // look for a suitable pivot:
        for (Index k = j + 1; k < n; ++k)
            if (abs(A[k][j]) > abs(A[pivot_row][j])) pivot_row = k;

        // swap the rows if we found a better pivot:
        if (pivot_row != j)
        {
            std::swap(A[j], A[pivot_row]);
            std::swap(b[j], b[pivot_row]);
        }

        // elimination:
        for (Index i = j + 1; i < n; ++i) {
            const double pivot = A[j][j];
            if (pivot == 0) error("can't solve: pivot==0");
            const double mult = A[i][j] / pivot;
            for (Index k = j; k < n; ++k)       //k 是列索引，索引A(i,k)
            {
                A[i][k] += -mult * A[j][k];
            }
            b[i] -= mult * b[j];
        }
    }
}

//------------------------------------------------------------------------------

Vector back_substitution(const Matrix& A, const Vector& b)
{
    const Index n = A.size();
    Vector x(n);

    for (Index i = n - 1; i >= 0; --i) {
        double s = 0.0;
        for (Index j = i + 1; j < A[0].size(); ++j)
            s += A[i][j] * x[j];
        s = b[i] - s;
        if (double m = A[i][i])
            x[i] = s / m;
        else
            //throw Back_subst_failure("Back substitution failure in row " + to_string(i));   //第一版 classical_elimination 下的异常
            error("can't solve: diagonal element is zero during back substitution");        //第二版 elim_with_partial_pivot 下的异常
    }

    return x;
}

//------------------------------------------------------------------------------

Vector random_vector(Index n)
{
    Vector v(n);
    const int max = 50;
    default_random_engine ran{ };                    //生成整数
    uniform_real_distribution<> ureal{ 0,max };     //将 int 映射为 [0:max) 中的 double
    for (Index i = 0; i < n; ++i)
        v[i] = ureal(ran);

    return v;
}

//------------------------------------------------------------------------------

Vector random_vector_with_seed(Index n, unsigned seed)
{
    Vector v(n);
    const int max = 50;
    default_random_engine ran{ seed };                    //生成整数
    uniform_real_distribution<> ureal{ 0,max };     //将 int 映射为 [0:max) 中的 double
    for (Index i = 0; i < n; ++i)
        v[i] = ureal(ran);

    return v;
}

//------------------------------------------------------------------------------

Matrix random_matrix(Index n)
{
    Matrix m(n);

    for (Index i = 0; i < n; ++i)
        m[i] = random_vector_with_seed(n, i);

    return m;
}

//------------------------------------------------------------------------------

Vector operator*(const Matrix& m, const Vector& u)
{
    const Index n = m.size();
    Vector v(n);    //double 默认初始化为0
    for (Index i = 0; i < n; ++i)
        for (Index j = 0; j < m[i].size(); ++j)
            v[i] += m[i][j] * u[j];
    return v;
}

//------------------------------------------------------------------------------

ostream& operator<<(ostream& os, const Vector& v)
{
    os << "{ ";
    for (const auto& x : v)
        os << x << ' ';

    os << '}';

    return os;
}

//------------------------------------------------------------------------------

ostream& operator<<(ostream& os, const Matrix& m)
{
    os << "{\n";
    for (const auto& x : m)
    {
        os << x << '\n';
    }
    os << '\n';
    return os;
}

//------------------------------------------------------------------------------

void solve_random_system(Index n)
{
    Matrix A = random_matrix(n);    // see ?4.6
    Vector b = random_vector(n);

    cout << "A = " << A << endl;
    cout << "b = " << b << endl;

    try {
        Vector x = classical_gaussian_elimination(A, b);
        cout << "classical elim solution is x = " << x << endl;
        Vector y = pivotal_elimination(A, b);
        cout << "pivotal elim solution is y = " << y << endl;
        Vector v = A * x;
        cout << " A * x = " << v << endl;
        Vector w = A * y;
        cout << " A * y = " << w << endl;
    }
    catch (const exception& e) {
        cerr << e.what() << std::endl;
    }
}

//------------------------------------------------------------------------------

int main()
try
{
    //double a[3][3] =
    //{
    //   {1,2,3},
    //   {2,3,4},
    //   {3,4,1}
    //};
    //double b[3] = { 14,20,14 };

    //double a[2][2] =
    //{
    //   {0,1},
    //   {1,0}
    //};
    //double b[2] = { 5,6 };

    //Matrix A(a);
    //Vector B(b);

    //cout << "Solving A*x=B" << endl;
    //cout << "A =\n" << A << endl;
    //cout << "B = " << B << endl;

    //Vector x = classical_gaussian_elimination(A, B);
    //Vector y = pivotal_elimination(A, B);
    //cout << "x = " << x << endl;
    //cout << "y = " << y << endl;


    solve_random_system(3);
    solve_random_system(4);
    solve_random_system(5);

    return 0;
}
catch (exception& e) {
    cerr << "error: " << e.what() << '\n';
    return 1;
}
catch (...) {
    cerr << "Oops: unknown exception!\n";
    return 2;
}



inline void error(const string& errormessage)
{
    cerr << errormessage << endl;
    throw runtime_error(errormessage);
}

//------------------------------------------------------------------------------

inline void error(string s1, string s2)
{
    error(s1 + s2);
}

//------------------------------------------------------------------------------

inline void error(string s1, int n)
{
    ostringstream os;
    os << n;
    error(s1, os.str());
}

24.9 利用decltype推导模板实参

#include"../../Matrix.h"

using namespace Numeric_lib;

template<typename F, typename A>
auto my_apply(F f, A x)
{
    Matrix<decltype(f(x[0]))> res(x.size());

    for (auto i = 0; i < x.size(); ++i) 
        res.data()[i] = f(x.data()[i]);     //在不知道维数的情况下，利用 矩阵在内存中的 行主序 排列。

    return xfer(res);
}

24.10

#include"../../std_lib_facilities.h"


using namespace std;

//------------------------------------------------------------------------------


//------------------------------------------------------------------------------

int main()
try
{
    int n{ 0 }, d{ 0 };

    cout << "Enter two integers n and d, ctrl + z to quit\n";
    while (cin >> n >> d)
    {
        if (n <= 0 || d <= 0)
        {
            cout << "n or d should greater than zero\n";
        }
        else
        {
            vector<int> vi(n);  //默认初始化为0
            while (d-- > 0)
            {
                ++vi[randint(n-1)];     //randint的范围是[0,n-1]，最大随机数是其输入参数
            }
            for (int i = 0; i < n; ++i)
            {
                cout << i << ":  ";
                for (int j = 0; j < vi[i]; ++j)
                    cout << '*';
                cout << "  " << vi[i] << '\n';
            }
        }

        cout << "try again\n";
    }
    return 0;
}
catch (exception& e) {
    cerr << "error: " << e.what() << '\n';
    return 1;
}
catch (...) {
    cerr << "Oops: unknown exception!\n";
    return 2;
}

24.11

#include <iostream>
#include<algorithm>
#include "../../Matrix.h"



using namespace Numeric_lib;

//成员函数 swap_column 已写于 Matrix.h 中

//1维
void swap_column(Index i, Index j)
{
    if (i == j) return;
    //一维矩阵（一维向量）只有一行
    std::swap((*this)(i), (*this)(j));
}

//2维
void swap_column(Index i, Index j)
{
    if (i == j) return;

    Index r_max = dim1();
    for (Index ri = 0; ri < r_max; ++ri)
        std::swap((*this)(ri, i), (*this)(ri, j));
}


//非成员函数 swap_column
template <class T, int N = 1>
void swap_columns(Matrix<T, N>& m, Index i, Index j)
{
    if (i == j) return;

    std::swap(m(i), m(j));
}

template <class T, int N = 2>
void swap_columns(Matrix<T, N>& m, Index i, Index j)
{
    if (i == j) return;

    Index r_max = m.dim1();
    for (Index ri = 0; ri < r_max; ++ri)
        std::swap(m(ri, i), m(ri, j));
}

24.12 二维矩阵和一维矩阵的乘法，两个相同维度矩阵的加法

#include <iostream>
#include <sstream>
#include <stdexcept>

#include "../../Matrix.h"
#include "../../MatrixIO.h"

using namespace Numeric_lib;

Matrix<double> operator*(Matrix<double, 2>& m2d, Matrix<double>& m1d)
{
	if (m2d.dim2() != m1d.dim1())
		error("二维矩阵的列数和一维矩阵的行数不匹配");

	const Index n = m2d.dim1();		//二维矩阵行数
	Matrix<double> res(n);		//存放结果

	for (Index i = 0; i < n; ++i)
		res(i) = dot_product(m2d[i], m1d);

	return res;
}

template<int N>
Matrix<double, N> operator+(Matrix<double, N>& m1, Matrix<double, N>& m2)
{
	//两个矩阵的维数应该同为N，但其余N-1维度可能不同
    // 以下尝试不可行，无法通过编译
	//switch (N) {
	//case 1:
	//	if (m1.dim1() != m2.dim1())
	//		error("两个一维矩阵不匹配");
	//	break;
	//case 2:
	//	if (m1.dim2() != m1.dim2())
	//		error("两个二维矩阵不匹配");
	//	break;
	//case 3:
	//	if (m1.dim2() != m1.dim2() || m1.dim3() != m1.dim3())
	//		error("两个三维矩阵不匹配");
	//	break;
	//default:
	//	error("至多三维矩阵");
	//	break;
	//}
    //换一种方法，只判断元素数量是否一致
    if (m1.size() != m2.size())
        error("两个矩阵不匹配");

	Matrix<double, N> res(m1);
	const Index sz = m1.size();

	//利用矩阵在内存中的“行主序”排列
	double* pres = res.data();
	double* p2 = m2.data();
	for (Index i = 0; i < sz; ++i)
		pres[i] += p2[i];

	return res;
}


int main()
try {

    // single dimension swap
    Matrix<double> m1(5);
    for (auto i = 0; i < m1.dim1(); ++i)
        m1(i) = i;

    std::cout << "Initial state of m1: " << m1 << '\n';

    // two dimension swap
    Matrix<double, 2> m2(4, 5);
    for (auto i = 0; i < m2.dim1(); ++i)
        for (auto j = 0; j < m2.dim2(); ++j)
            m2(i, j) = 10 + i * 10 + j;

    std::cout << "Initial state of m2: " << m2 << '\n';


    // multi-dimension addition
    Matrix<double, 2> m3(4, 5);         // same size as m2
    m3 = 3;                             // set all values to 3
    Matrix<double, 2> m4 = m2 + m3;
    std::cout << "Results of addition: " << m4 << '\n';

    // single dimension addition
    Matrix<double> m5(5);
    m5 = 5;
    Matrix<double> m6 = m1 + m5;
    std::cout << "Results of addition: " << m6 << '\n';

    // Matrix multiplication
    auto m7 = m2 * m1;
    std::cout << "Results of multiplication: " << m7 << '\n';

    // more fun and tests
    double ad1[][2]{ { 2, -5 }, { -2, 4 } };
    double ad2[2]{ 7, -6 };
    Matrix<double, 2> m8{ ad1 };
    Matrix<double> m9{ ad2 };
    auto m10 = m8 * m9;
    std::cout << "More multiplication: " << m10 << '\n';
}
catch (std::exception& e) {
    std::cerr << "Exception: " << e.what() << '\n';
    return 1;
}
catch (...) {
    std::cerr << "Unknown exception\n";
    return 2;
}