趣味程序 - 分数类与矩阵类的实现

来源:互联网 发布:nat端口复用 编辑:程序博客网 时间:2024/06/14 19:12

分数和矩阵其实都是Matlab玩剩下的,程序本身纯粹算是C++的一个练习,实现一些基本的功能接口。大体来说就是利用类和重载,最主要的还是体会C++代码复用、封装和标准化接口的思想。C++单论代码效率上当然无法和C相提并论,但是从工程的角度上来说,C++堪称优秀。

唯一的亮点(至少个人看来)在于,行列式的值算法其实很值得研究,如果按照行列式按行/列展开计算则需要递归,复杂度高达 O(n!)。这毫无疑问不是我们可以忍受的效率。我在这里利用了行列式行间相消的特点将行列式先简化为一个上三角行列式,然后求对角元素之积即可。时间开销集中在简化行列式的过程中的三重循环,复杂度是 O(n3)。思想很朴素,也没有多做优化,如果对读者能够有所启发就再好不过。

不过现在行列式计算仍然没有非常高效的算法,效率最高的 Fast matrix multiplication[1]需要的时间复杂度为O(n2.373)。此外比较有名的还有Strassen算法[2]。诸如Matlab这类专业数学工具则对这类计算在利用有名的算法的基础上进行了非常多的优化,甚至优化到编译器和操作系统的级别,基本上是 O(n2+x)(0<x<1)的级别。事实上因为求行列式的值必须是要遍历每一个元素的,因此这个问题本身的复杂度下界就是 O(n2)。有兴趣的朋友可以深入研究。下附代码。

/** MIT license* Copyright © 2016 Maristie* * Permission is hereby granted, free of charge, to any person obtaining a copy of this software* and associated documentation files (the "Software"), to deal in the Software without restriction,* including without limitation the rights to use, copy, modify, merge, publish, distribute, * sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, * subject to the following conditions:* The above copyright notice and this permission notice shall be included in all copies or * substantial portions of the Software.* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED* TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE* OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.*/#include<iostream>#include<math.h>#include<time.h>using namespace std;class Fraction;class Matrix;long long findGCD(long long, long long); //find greatest common divisorostream& operator<<(ostream&, Fraction&);istream& operator >> (istream&, Fraction&);ostream& operator<<(ostream&, Matrix&);istream& operator >> (istream&, Matrix&);Matrix operator*(int, Matrix);class Fraction//create Fraction class{private:    long long Up, Down; //numerator and denominator    int whatSign();//judge sign of fraction    bool isInt();//simplize fraction and judge if it's a intpublic:    Fraction();    Fraction(int, int);    Fraction(const Fraction&);//3 constructors    void Simplize();//simplize fraction by divide GCD    Fraction operator+(Fraction&);    Fraction operator-(Fraction&);    Fraction operator*(Fraction&);    Fraction operator/(Fraction&);//4 overloaded functions    Fraction operator*(int);    bool operator==(Fraction);    bool operator==(int);    bool operator==(long long);    bool operator>(Fraction);    bool operator<(Fraction);    bool operator>=(Fraction);    bool operator<=(Fraction);    friend ostream& operator<<(ostream&, Fraction&);    friend istream& operator >> (istream&, Fraction&);    friend class Matrix;};class Matrix//create Matrix class{private:    int Row, Col;    Fraction** Data;    Matrix AlgCof(int, int);//get algebraic cofactor    Fraction detAlgCof(int, int);//get the determinant value of algebraic cofactor    int maxIndexRow();//get row of the max element    int maxIndexCol();//get col of the max elementpublic:    Matrix();    Matrix(int, int);    Matrix(const Matrix&);    ~Matrix();    void setFree();//set Data pointer to free status    void Reset(int, int);    void Reset(const Matrix&);//reset Data pointer, or rather reallocate space for Data    Fraction getElem(int, int); //get element    void putElem(int, int, Fraction); //put element in    void operator=(Matrix&);//matrix assignment    bool operator==(Matrix&);    bool operator!=(Matrix&);//compare functions    Matrix operator+(Matrix&);    Matrix operator-(Matrix&);    Matrix operator*(Matrix&);    Matrix operator*(int);    friend Matrix operator*(int, Matrix); //5 overloaded calculating functions    Matrix Power(int);    void resize(int, int);    void swaprows(int, int);    void swapcols(int, int);    int size();//return size    friend ostream& operator<<(ostream&, Matrix&);    friend istream& operator >> (istream&, Matrix&);    void transpose();//the result will be output    void diag();//here I'll output all elements on the diagonal in order if Row equals Col    Fraction det();//a fraction will be returned, and the algorithm complexity is O(n^3)    Matrix inv();//a matrix will be returned as the inverse matrix    void norm();//the output type is double    Fraction max();    void max_index();//output the index in the format (x,y), x=row+1, y=col+1};Fraction::Fraction(){    Up = 0;    Down = 1;}Fraction::Fraction(int inputUp, int inputDown){    if (inputDown != 0)    {        Up = inputUp;        Down = inputDown;        Simplize();    }    else        cout << "Denominator can't be zero." << endl;}Fraction::Fraction(const Fraction& origFrac){    Up = origFrac.Up;    Down = origFrac.Down;}Fraction Fraction::operator+(Fraction& origFrac){    Fraction Res;    Res.Up = Up * origFrac.Down + Down * origFrac.Up;    Res.Down = Down * origFrac.Down;    Res.Simplize();    return Res;}Fraction Fraction::operator-(Fraction& origFrac){    Fraction Res;    Res.Up = Up * origFrac.Down - Down * origFrac.Up;    Res.Down = Down * origFrac.Down;    Res.Simplize();    return Res;}Fraction Fraction::operator*(Fraction& origFrac){    Fraction Res;    Res.Up = Up * origFrac.Up;    Res.Down = Down * origFrac.Down;    Res.Simplize();    return Res;}Fraction Fraction::operator/(Fraction& origFrac){    if (origFrac == 0)    {        cout << "0 is divided." << endl;        return origFrac;    }    else    {        Fraction Res;        Res.Up = Up * origFrac.Down;        Res.Down = Down * origFrac.Up;        Res.Simplize();        return Res;    }}Fraction Fraction::operator*(int Num){    Fraction Res(*this);    Res.Up *= Num;    Res.Simplize();    return Res;}int Fraction::whatSign(){    if (Up * Down == 0)        return 0;    else if (Up * Down > 0)        return 1;    else return -1;}void Fraction::Simplize(){    if (whatSign() == 0)    {        Down = 1;        return;    }    else    {        int GCD = findGCD(abs(Up), abs(Down));        Up = whatSign() * abs(Up) / GCD;        Down = abs(Down) / GCD;        return;    }}bool Fraction::operator==(Fraction anothFrac){    Fraction myFrac(*this);    if (myFrac.Up == anothFrac.Up && myFrac.Down == anothFrac.Down)        return true;    else return false;}bool Fraction::operator==(int Num){    if (Up == Down * Num)        return true;    else return false;}bool Fraction::operator==(long long Num){    if (Up == Down * Num)        return true;    else return false;}bool Fraction::operator>(Fraction origFrac){    Fraction myFrac(*this);    if (myFrac.Up * origFrac.Down > origFrac.Up * myFrac.Down)        return true;    else return false;}bool Fraction::operator>=(Fraction origFrac){    Fraction myFrac(*this);    if (myFrac.Up * origFrac.Down >= origFrac.Up * myFrac.Down)        return true;    else return false;}bool Fraction::operator<(Fraction origFrac){    Fraction myFrac(*this);    if (myFrac.Up * origFrac.Down < origFrac.Up * myFrac.Down)        return true;    else return false;}bool Fraction::operator<=(Fraction origFrac){    Fraction myFrac(*this);    if (myFrac.Up * origFrac.Down <= origFrac.Up * myFrac.Down)        return true;    else return false;}bool Fraction::isInt(){    if (Up == 0 || Down == 1)        return true;    else return false;}ostream& operator<<(ostream& cout, Fraction& myFrac){    if (myFrac.isInt() == true)        cout << myFrac.Up;    else        cout << myFrac.Up << '/' << myFrac.Down;    return cout;}istream& operator >> (istream& cin, Fraction& myFrac){    char Input[1000] = { 0 };    int isFrac = 0;    int Pos = 0;    cin >> Input;    int InputLen = strlen(Input);    for (Pos; Pos < InputLen; ++Pos)        if (Input[Pos] == '/')        {            isFrac = 1;            break;        }    if (isFrac == 0)    {        myFrac.Up = atoi(Input);        myFrac.Down = 1;    }    else    {        Input[Pos] = 0;        char InputUp[1000] = { 0 };        strcpy(InputUp, Input);        char InputDown[1000] = { 0 };        strcpy(InputDown, Input + Pos + 1);        myFrac.Up = atoi(InputUp);        myFrac.Down = atoi(InputDown);    }    return cin;}long long findGCD(long long Num1, long long Num2) //find GCD by Euclidean Algorithm{    if (Num1 < Num2)    {        long long SwapTemp = Num1;        Num1 = Num2;        Num2 = SwapTemp;    }    long long Remain = Num1 % Num2;    while (Remain)    {        Num1 = Num2;        Num2 = Remain;        Remain = Num1 % Num2;    }    return Num2;}Matrix::Matrix(){    Row = Col = 0;    Data = NULL;}Matrix::Matrix(int InRow, int InCol){    Reset(InRow, InCol);}Matrix::Matrix(const Matrix& origMat){    Reset(origMat);}void Matrix::Reset(int InRow, int InCol){    Row = InRow;    Col = InCol;    Data = new Fraction*[Row];    for (int i = 0; i < Row; ++i)        Data[i] = new Fraction[Col];    for (int i = 0; i < Row; ++i)        for (int j = 0; j < Col; ++j)            Data[i][j] = Fraction(0, 1);}void Matrix::Reset(const Matrix& origMat){    Row = origMat.Row;    Col = origMat.Col;    Data = new Fraction*[Row];    for (int i = 0; i < Row; ++i)    {        Data[i] = new Fraction[Col];        for (int j = 0; j < Col; ++j)            Data[i][j] = origMat.Data[i][j];    }}Matrix::~Matrix(){    setFree();}void Matrix::setFree(){    for (int i = 0; i < Row; ++i)        delete[] Data[i];    delete[] Data;    Data = NULL;}Fraction Matrix::getElem(int i, int j){    return Data[i][j];}void Matrix::putElem(int i, int j, Fraction origFrac){    Data[i][j] = origFrac;    return;}void Matrix::operator=(Matrix& origMat){    setFree();    Reset(origMat.Row, origMat.Col);    for (int i = 0; i < Row; ++i)        for (int j = 0; j < Col; ++j)            Data[i][j] = origMat.Data[i][j];    return;}bool Matrix::operator==(Matrix& origMat){    if (Row == origMat.Row && Col == origMat.Col)    {        for (int i = 0; i < Row; ++i)            for (int j = 0; j < Col; ++j)                if (!(Data[i][j] == origMat.Data[i][j]))                    return false;        return false;    }    else cout << "Illegal input." << endl;    return true;}bool Matrix::operator!=(Matrix& origMat){    if (Row == origMat.Row && Col == origMat.Col)    {        for (int i = 0; i < Row; ++i)            for (int j = 0; j < Col; ++j)                if (!(Data[i][j] == origMat.Data[i][j]))                    return true;        return false;    }    else cout << "Illegal input." << endl;    return false;}Matrix Matrix::operator+(Matrix& origMat){    if (Row == origMat.Row && Col == origMat.Col)    {        Matrix resMat(Row, Col);        for (int i = 0; i < Row; ++i)            for (int j = 0; j < Col; ++j)                resMat.Data[i][j] = Data[i][j] + origMat.Data[i][j];        return resMat;    }    else cout << "Illegal input." << endl;    return origMat;}Matrix Matrix::operator-(Matrix& origMat){    if (Row == origMat.Row && Col == origMat.Col)    {        Matrix resMat(Row, Col);        for (int i = 0; i < Row; ++i)            for (int j = 0; j < Col; ++j)                resMat.Data[i][j] = Data[i][j] - origMat.Data[i][j];        return resMat;    }    else cout << "Illegal input." << endl;    return origMat;}Matrix Matrix::operator*(Matrix& origMat){    if (Col == origMat.Row)    {        Matrix resMat(Row, origMat.Col);        for (int i = 0; i < resMat.Row; ++i)            for (int j = 0; j < resMat.Col; ++j)                for (int k = 0; k < Col; ++k)                    resMat.Data[i][j] = resMat.Data[i][j] + Data[i][k] * origMat.Data[k][j];        return resMat;    }    else cout << "Illegal input." << endl;    return origMat;}Matrix Matrix::operator*(int Num){    Matrix resMat(Row, Col);    for (int i = 0; i < Row; ++i)        for (int j = 0; j < Col; ++j)            resMat.Data[i][j] = Data[i][j] * Num;    return resMat;}Matrix operator*(int Num, Matrix origMat){    Matrix resMat(origMat.Row, origMat.Col);    for (int i = 0; i < origMat.Row; ++i)        for (int j = 0; j < origMat.Col; ++j)            resMat.Data[i][j] = origMat.Data[i][j] * Num;    return resMat;}Matrix Matrix::Power(int n){    if (Row != Col || n <= 0)    {        cout << "Illegal input." << endl;        return *this;    }    else    {        Matrix resMat(*this);        for (int i = 1; i < n; ++i)            resMat = resMat * (*this);        return resMat;    }}void Matrix::resize(int sizeRow, int sizeCol){    if (Row * Col != sizeRow * sizeCol)        cout << "Illegal input." << endl;    else    {        Fraction* MemStor = new Fraction[Row * Col];        for (int i = 0; i < Row; ++i)            for (int j = 0; j < Col; ++j)                MemStor[i * Col + j] = Data[i][j];        setFree();        Reset(sizeRow, sizeCol);        for (int i = 0; i < Row; ++i)            for (int j = 0; j < Col; ++j)                Data[i][j] = MemStor[i * Col + j];        delete[] MemStor;        MemStor = NULL;    }    return;}void Matrix::swaprows(int Row1, int Row2){    Fraction swapTemp;    for (int i = 0; i < Col; ++i)    {        swapTemp = Data[Row1][i];        Data[Row1][i] = Data[Row2][i];        Data[Row2][i] = swapTemp;    }    return;}void Matrix::swapcols(int Col1, int Col2){    Fraction swapTemp;    for (int j = 0; j < Row; ++j)    {        swapTemp = Data[j][Col1];        Data[j][Col1] = Data[j][Col2];        Data[j][Col2] = swapTemp;    }    return;}int Matrix::size(){    return Row * Col;}istream& operator >> (istream& cin, Matrix& myMat){    for (int i = 0; i < myMat.Row; ++i)        for (int j = 0; j < myMat.Col; ++j)            cin >> myMat.Data[i][j];    return cin;}ostream& operator<<(ostream& cout, Matrix& myMat){    for (int i = 0; i < myMat.Row; ++i)    {        for (int j = 0; j < myMat.Col; ++j)            cout << myMat.Data[i][j] << '\t';        cout << endl;    }    return cout;}void Matrix::transpose(){    Fraction* memStor = new Fraction[Row * Col];    for (int i = 0; i < Row; ++i)        for (int j = 0; j < Col; ++j)            memStor[i * Col + j] = Data[i][j];    setFree();    Reset(Col, Row);    for (int i = 0; i < Col; ++i)        for (int j = 0; j < Row; ++j)            Data[j][i] = memStor[i * Row + j];    delete[] memStor;    memStor = NULL;    cout << "After transpose, the matrix turns to:" << endl;    cout << *this << endl;    return;}void Matrix::diag(){    if (Row != Col)    {        cout << "No diag element existing." << endl;        return;    }    int i;    cout << "Output the diag elements in order:" << endl;    for (i = 0; i<Row - 1; ++i)        cout << getElem(i, i) << ' ';    cout << getElem(i, i) << endl;    return;}Fraction Matrix::det(){    Matrix myMat(*this);    if (myMat.Row != myMat.Col)    {        cout << "Row not equal to Col." << endl;        return myMat.Data[0][0];    }    else    {        int i, j, k;        int SwapTime = 0;        for (i = 0; i < myMat.Row; ++i)        {            if (myMat.Data[i][i] == 0)            {                int isDetZero = 1;                for (j = i + 1; j < myMat.Row; ++j)                    if (!(myMat.Data[j][i] == 0))                    {                        myMat.swaprows(i, j);                        ++SwapTime;                        isDetZero = 0;                        break;                    }                if (isDetZero == 1)                {                    Fraction resFrac(0, 1);                    return resFrac;                }            }            for (j = i + 1; j < myMat.Row; ++j)                for (k = myMat.Col - 1; k >= i; --k)                    myMat.Data[j][k] = myMat.Data[j][k] - myMat.Data[i][k] * myMat.Data[j][i] / myMat.Data[i][i];        }        Fraction resFrac(1, 1);        for (i = 0; i < Row; ++i)            resFrac = resFrac * myMat.Data[i][i];        if (SwapTime % 2 == 1)            resFrac = resFrac * (-1);        return resFrac;    }}Matrix Matrix::AlgCof(int InRow, int InCol){    Matrix resMat(Row - 1, Col - 1);    int resRow = 0, resCol = 0;    for (int i = 0; i < Row; ++i)    {        if (i == InRow)            continue;        for (int j = 0; j < Col; ++j)        {            if (j == InCol)                continue;            resMat.Data[resRow][resCol] = Data[i][j];            ++resCol;        }        ++resRow;        resCol = 0;    }    return resMat;}Fraction Matrix::detAlgCof(int InRow, int InCol){    Matrix AlgMat = AlgCof(InRow, InCol);    Fraction detValue = AlgMat.det();    if ((InRow + InCol) % 2 == 1)        detValue = detValue * -1;    return detValue;}Matrix Matrix::inv(){    if (Row != Col)    {        cout << "Row unequal to Col." << endl;        return *this;    }    else    {        Matrix resMat(Row, Col);        Fraction detVal = det();        if (detVal == 0)        {            cout << "No inverse matrix." << endl;            return *this;        }        for (int i = 0; i < Row; ++i)            for (int j = 0; j < Col; ++j)                resMat.Data[i][j] = detAlgCof(j, i) / detVal;        return resMat;    }}void Matrix::norm(){    Fraction Sum;    for (int i = 0; i < Row; ++i)        for (int j = 0; j < Col; ++j)            Sum = Sum + Data[i][j] * Data[i][j];    cout << "The norm value is:" << endl;    cout << sqrt(1.0*Sum.Up / Sum.Down) << endl << endl;    return;}Fraction Matrix::max(){    Fraction Res = Data[0][0];    for (int i = 0; i < Row; ++i)        for (int j = 0; j < Col; ++j)            if (Data[i][j] > Res)                Res = Data[i][j];    return Res;}int Matrix::maxIndexRow(){    int Val = 0;    Fraction Res = Data[0][0];    for (int i = 0; i < Row; ++i)        for (int j = 0; j < Col; ++j)            if (Data[i][j] > Res)            {                Res = Data[i][j];                Val = i;            }    return Val;}int Matrix::maxIndexCol(){    int Val = 0;    Fraction Res = Data[0][0];    for (int i = 0; i < Row; ++i)        for (int j = 0; j < Col; ++j)            if (Data[i][j] > Res)            {                Res = Data[i][j];                Val = j;            }    return Val;}void Matrix::max_index(){    int row = maxIndexRow();    int col = maxIndexCol();    cout << "The max index is: (" << row + 1 << ',' << col + 1 << ")." << endl;    return;}int main(){    Matrix a(4, 4);    srand((unsigned int)time(0));    for (int i = 0; i < 4; ++i)        for (int j = 0; j < 4; ++j)            a.putElem(i, j, Fraction(rand() % 100 + 10 * i + j + 1, rand() % 100 + 10 * j + i + 1));    cout << a << endl << a.inv() << endl << a.det() << endl << endl;    a.transpose();    a.norm();    cout << a.inv().inv() << endl;    a.transpose();    return 0;}

参考资料:
[1] - Computational complexity of mathematical operations, Wikipedia, https://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations
[2] - Strassen algorithm, Wikipedia, https://en.wikipedia.org/wiki/Strassen_algorithm

0 0