24. Numerics

24.1. Introduction

수치 계산은 중요하다.

24.2. Size, precision, and overflow

정수형과 부동 소수점 자료형은 정수와 실수에 대한 근사만을 할 수 있기 때문에 다음과 같은 부정확한 계산이 나온다.

float x = 1.0/333;
float sum = 0;
for (int i=0; i<333; ++i) sum += x;
std::cout << std::setprecision(15) << sum << "\n";

이는 1.0이 정확히 나오지 않는다.

정수형에 대해서는 다음과 같이 오버플로우가 나는 문제가 있다.

short int y = 40000;
int i = 1000000;
std::cout << y << "   " << i*i << "\n";

C++은 여러 산술 자료형들을 제공하며 그 크기는 일정하지 않다. 보통은 그냥 char, int, double을 중심으로 쓰는 것이 좋다.

정수형을 부동 소수점 자료형에 할당하면 경우에 따라 정밀도를 잃는다. 부동 소수점 자료형을 정수형에 할당하면 소수점 부분을 잃는다. 이는 버림으로 적용된다.

void f(int i, double fpd)
{
          char c = i;            // yes: chars really are very small integers
          short s = i;          // beware: an int may not fit in a short int
          i = i+1;                // what if i was the largest int?
          long lg = i*i;      // beware: a long may not be any larger than an int
          float fps = fpd;  // beware: a large double may not fit in a float
          i = fpd;               // truncates: e.g., 5.7 –> 5
          fps = i;               // you can lose precision (for very large int values)
}

24.2.1. Numeric limits

<limits>에는 C++에서 기본 자료형에 대한 구현체 의존적 특성들이 저장되어 있다. 보통 다음과 같다.

std::cout << "number of bytes in an int: " << sizeof(int) << '\n';
std::cout << "largest int: " << INT_MAX << '\n';
std::cout << "smallest int value: " << std::numeric_limits<int>::min() << '\n';

if (std::numeric_limits<char>::is_signed)
          cout << "char is signed\n";
else
          cout << "char is unsigned\n";

char ch = std::numeric_limits<char>::min() ;         // smallest positive value
cout << "the char with the smallest positive value: " << ch << '\n';
cout << "the int value of the char with the smallest positive value: "
          << int(ch) << '\n';

24.3. Arrays

배열은 인덱스로 항목을 접근할 수 있는 항목의 나열이다. N차원 배열의 구현에 대해 생각해 보자.

24.4. C-style multidimensional arrays

C 스타일의 다차원 배열이 C++에도 똑같이 존재한다.

int ai[4];                     // 1-dimensional array
double ad[3][4];      // 2-dimensional array
char ac[3][4][5];       // 3-dimensional array
ai[1] = 7;
ad[2][3] = 7.2;
ac[2][3][4] = 'c';

이는 버그에 취약하기 때문에 쓰지 않는 것이 좋다. 포인터로의 묵시적 형변환 등의 문제점이 있기 때문이다.

24.5. The Matrix library

행렬 라이브러리를 생각해 보자.

#ifndef MATRIX_LIB
#define MATRIX_LIB

#include<string>
#include<algorithm>

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));
    }
};

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

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));
    }
};

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

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

이는 다음과 같이 사용된다.

#include "Matrix.h"
using namespace Numeric_lib;

void f(int n1, int n2, int n3)
{
          Matrix<double,1> ad1(n1);        // elements are doubles; one dimension
          Matrix<int,1> ai1(n1);                  // elements are ints; one dimension
          ad1(7) = 0;                                      // subscript using ( ) — Fortran style
          ad1[7] = 8;                                      // [ ] also works — C style

          Matrix<double,2> ad2(n1,n2);          // 2-dimensional
          Matrix<double,3> ad3(n1,n2,n3);    // 3-dimensional
          ad2(3,4) = 7.5;                                       // true multidimensional subscripting
          ad3(3,4,5) = 9.2;
}

void f(int n1, int n2, int n3)
{
          Matrix<int,0> ai0;           // error: no 0D matrices

          Matrix<double,1> ad1(5);
          Matrix<int,1> ai(5);
          Matrix<double,1> ad11(7);

          ad1(7) = 0;                       // Matrix_error exception (7 is out of range)
          ad1 = ai;                           // error: different element types
          ad1 = ad11;                     // Matrix_error exception (different dimensions)

          Matrix<double,2> ad2(n1);         // error: length of 2nd dimension missing
          ad2(3) = 7.5;                                   // error: wrong number of subscripts
          ad2(1,2,3) = 7.5;                            // error: wrong number of subscripts

          Matrix<double,3> ad3(n1,n2,n3);
          Matrix<double,3> ad33(n1,n2,n3);
          ad3 = ad33;                     // OK: same element type, same dimensions
}

oid init(Matrix<int,2>& a)       // initialize each element to a characteristic value
{
          for (int i=0; i<a.dim1(); ++i)
                    for (int j = 0; j<a.dim2(); ++j)
                              a(i,j) = 10*i+j;
}

void print(const Matrix<int,2>& a)    // print the elements row by row
{
          for (int i=0; i<a.dim1(); ++i) {
                    for (int j = 0; j<a.dim2(); ++j)
                              std::cout << a(i,j) <<'\t';
                    std::cout << '\n';
          }
}

24.5.2. 1D Matrix

1차원 행렬은 다음과 같이 쓸 수 있다.

Matrix<int,1> a1(8);     // a1 is a 1D Matrix of ints
Matrix<int> a(8);          // means Matrix<int,1> a(8);
a.size();                 // number of elements in Matrix
a.dim1();               // number of elements in 1st dimension
int* p = a.data();        // extract data as a pointer to an array
a(i);            // ith element (Fortran style), but range checked
a[i];            // ith element (C style), range checked
a(1,2);       // error: a is a 1D Matrix
a.slice(i);         // the elements from a[i] to the last
a.slice(i,n);     // the n elements from a[i] to a[i+n–1]
a.slice(4,4) = a.slice(0,4);       // assign first half of a to second half
a.slice(4) = a.slice(0,4);       // assign first half of a to second half
Matrix<int> a2 = a;         // copy initialization
a = a2;                               // copy assignment
a *= 7;                         // scaling: a[i]*=7 for each i (also +=, –=, /=, etc.)
a = 7;                          // a[i]=7 for each i
a.apply(f);               // a[i]=f(a[i]) for each element a[i]
a.apply(f,7);            // a[i]=f(a[i],7) for each element a[i]
b = a*7;                     // b[i] = a[i]*7 for each i
a *= 7;                       // a[i] = a[i]*7 for each i
y = apply(f,x);         // y[i] = f(x[i]) for each i
x.apply(f);               // x[i] = f(x[i]) for each i
b = apply(abs,a);    // make a new Matrix with b(i)==abs(a(i))
b = apply(f,a,x);                   // b[i]=f(a[i],x) for each i
double scale(double d, double s) { return d*s; }
b = apply(scale,a,7);                // b[i] = a[i]*7 for each i
void scale_in_place(double& d, double s) { d *= s; }
b.apply(scale_in_place,7);     // b[i] *= 7 for each i
Matrix<int> a3 = scale_and_add(a,8,a2);         // fused multiply and add
int r = dot_product(a3,a);                                   // dot product

void some_function(double* p, int n)
{
          double val[] = { 1.2, 2.3, 3.4, 4.5 };
          Matrix<double> data(p,n);
          Matrix<double> constants(val);
          // . . .
}

24.5.3. 2D Matrix

2차원 행렬은 다음과 같이 쓸 수 있다.

Matrix<int,2> a(3,4);

int s = a.size();                // number of elements
int d1 = a.dim1();           // number of elements in a row
int d2 = a.dim2();           // number of elements in a column
int* p = a.data();            // extract data as a pointer to a C-style array
a(i,j);                            // (i,j)th element (Fortran style), but range checked
a[i];                              // ith row (C style), range checked
a[i][j];                          // (i,j)th element (C style)
a.slice(i);                      // the rows from the a[i] to the last
a.slice(i,n);                  // the rows from the a[i] to the a[i+n–1]
Matrix<int,2> a2 = a;         // copy initialization
a = a2;                                  // copy assignment
a *= 7;                                  // scaling (and +=, –=, /=, etc.)
a.apply(f);                           // a(i,j)=f(a(i,j)) for each element a(i,j)
a.apply(f,7);                        // a(i,j)=f(a(i,j),7) for each element a(i,j)
b=apply(f,a);                      // make a new Matrix with b(i,j)==f(a(i,j))
b=apply(f,a,7);                  // make a new Matrix with b(i,j)==f(a(i,j),7)
a.swap_rows(1,2);            // swap rows a[1] <–> a[2]
enum Piece { none, pawn, knight, queen, king, bishop, rook };
Matrix<Piece,2> board(8,8);          // a chessboard

const int white_start_row = 0;
const int black_start_row = 7;

Matrix<Piece> start_row
          = {rook, knight, bishop, queen, king, bishop, knight, rook};

Matrix<Piece> clear_row(8) ;   // 8 elements of the default value
board[white_start_row] = start_row;                    // reset white pieces
for (int i = 1; i<7; ++i) board[i] = clear_row;      // clear middle of the board
board[black_start_row] = start_row;                   // reset black pieces

24.5.4. Matrix I/O

행렬 입출력은 다음과 같이 수행한다.

Matrix<double> a(4);
std::cin >> a;
std::cout << a;

해당하는 입출력 연산자가 정의되어 있다.

24.5.5. 3D Matrix

3차원 행렬은 다음과 같이 사용한다.

Matrix<int,3> a(10,20,30);

a.size();                            // number of elements
a.dim1();                          // number of elements in dimension 1
a.dim2();                          // number of elements in dimension 2
a.dim3();                          // number of elements in dimension 3
int* p = a.data();             // extract data as a pointer to a C-style array
a(i,j,k);                              // (i,j,k)th element (Fortran style), but range checked
a[i];                                    // ith row (C style), range checked
a[i][j][k];                            // (i,j,k)th element (C style)
a.slice(i);                           // the rows from the ith to the last
a.slice(i,j);                         // the rows from the ith to the jth
Matrix<int,3> a2 = a;      // copy initialization
a = a2;                                // copy assignment
a *= 7;                                // scaling (and +=, –=, /=, etc.)
a.apply(f);                          // a(i,j,k)=f(a(i,j,k)) for each element a(i,j,k)
a.apply(f,7);                       // a(i,j,k)=f(a(i,j,k),7) for each element a(i,j,k)
b=apply(f,a);                     // make a new Matrix with b(i,j,k)==f(a(i,j,k))
b=apply(f,a,7);                  // make a new Matrix with b(i,j,k)==f(a(i,j,k),7)
a.swap_rows(7,9);            // swap rows a[7] <–> a[9]

24.6. An example: solving linear equations

선형 방정식을 푸는 데 쓸 수 있다.

24.6.1. Classical Gaussian elimination

가우스 소거법을 다음과 같이 구현해 보자.

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

Vector classical_gaussian_elimination(Matrix A, Vector b)
{
          classical_elimination(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(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));
                              b(i) –= mult*b(j);       // make the corresponding change to b
                    }
          }
}

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));

                    if (double m = A(i,i))
                              x(i) = s/m;
                    else
                              throw Back_subst_failure(i);
          }

          return x;
}

24.6.2. Pivoting

피보팅이 서브루틴으로 쓰인다.

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 (std::abs(A(k,j)) > std::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));
                              b(i) –= mult*b(j);
                    }
          }
}

24.6.3. Testing

테스팅은 다음과 같이 수행한다.

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

          std::cout << "A = " << A << '\n';
          std::cout << "b = " << b << '\n';

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

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;
}

24.7. Random numbers

C++에서는 좋은 난수 생성기를 지원한다.

Vector random_vector(Index n)
{
          Vector v(n);
          std::default_random_engine ran{};                        // generates integers
          std::uniform_real_distribution<> ureal{0,max};   // maps ints into doubles
                                                                                                 // in [0:max)

          for (Index i = 0; i < n; ++i)
                    v(i) = ureal(ran);

          return v;
}

int randint(int min, int max)
{
          static std::default_random_engine ran;
          return std::uniform_int_distribution<>{min,max}(ran);
}

int randint(int max)
{
          return randint(0,max);
}

auto gen = std::bind(std::normal_distribution<double>{15,4.0},
                              std::default_random_engine{});

24.8. The standard mathematical functions

<cmath>에는 많은 수학 함수가 있다.

errno = 0;
double s2 = sqrt(–1);
if (errno) std::cerr << "something went wrong with something somewhere";
if (errno == EDOM)            // domain error
          std::cerr << "sqrt() not defined for negative argument";
std::pow(very_large,2);            // not a good idea
if (errno==ERANGE)         // range error
          std::cerr << "pow(" << very_large << ",2) too large for a double";

이는 errno로 에러를 리포트한다.

24.9. Complex numbers

C++에서는 복소수 std::complex를 지원한다.

24.10. References

더 많은 정보를 찾아보고 싶다면:

https://en.cppreference.com/w/cpp/numeric

답글 남기기

아래 항목을 채우거나 오른쪽 아이콘 중 하나를 클릭하여 로그 인 하세요:

WordPress.com 로고

WordPress.com의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

Google photo

Google의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

Twitter 사진

Twitter의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

Facebook 사진

Facebook의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

%s에 연결하는 중