Segmentation Fault in MAtrix Multiplication

Question

I am a pretty much beginner in C/C++. I was trying to implement a program to return the value of the n-th Fibonacci number using Matrix Exponentiation.I wrote the following code :

#include <cmath>
#include <iostream>

using namespace std;

typedef struct Matrix {
    unsigned long long a_00;
    unsigned long long a_01;
    unsigned long long a_10;
    unsigned long long a_11;
}Matrix;

Matrix* init_matrix(unsigned long long n);
Matrix* init_matrix(unsigned long long a_00, unsigned long long a_01,
                    unsigned long long a_10, unsigned long long a_11);
Matrix* mult_matrix(Matrix* A, Matrix* B);      
Matrix* expo_matrix(Matrix* A,unsigned int n);
unsigned long long fib_matrix_expo(unsigned int n); 

int main() {

    unsigned int n;
    cin >> n;

    unsigned long long result = fib_matrix_expo(n);
    cout << result << endl;
    return 0;
}


Matrix* init_matrix(unsigned long long n) {
    Matrix* ret_matrix = (Matrix*)::operator new (sizeof(Matrix));

    (ret_matrix->a_00) = n;
    (ret_matrix->a_01) = n;
    (ret_matrix->a_10) = n;
    (ret_matrix->a_11) = n;

    return ret_matrix;
}

Matrix* init_matrix(unsigned long long a_00, unsigned long long a_01,
                    unsigned long long a_10, unsigned long long a_11) {
/*==================================================================================*/
    Matrix* ret_matrix = (Matrix*)::operator new (sizeof(Matrix));

    (ret_matrix->a_00) = a_00;
    (ret_matrix->a_01) = a_01;
    (ret_matrix->a_10) = a_10;
    (ret_matrix->a_11) = a_11;

    return ret_matrix;  

}


Matrix* matrix_mult(Matrix* A, Matrix* B) {
    Matrix* result = (Matrix*)::operator new (sizeof(Matrix));

    (result->a_00) = ((A->a_00) * (B->a_00) + (A->a_01)*(B->a_10));
    (result->a_01) = ((A->a_00) * (B->a_01) + (A->a_01)*(B->a_11));
    (result->a_10) = ((A->a_10) * (B->a_00) + (A->a_11)*(B->a_10));
    (result->a_11) = ((A->a_10) * (B->a_01) + (A->a_11)*(B->a_11));

    return result;
}

Matrix* matrix_expo(Matrix* base_matrix, unsigned int index) {
    Matrix* result;

    if (index == 0) {
        result = init_matrix(1);
    }
    else if (index == 1) {
        result = base_matrix;
    }
    else {
        Matrix* temp = matrix_expo(base_matrix,static_cast<unsigned int>(index/2));
        result = matrix_mult( matrix_expo(temp,2), matrix_expo(base_matrix,(index % 2)) );
    }
    return result;
}

unsigned long long fib_matrix_expo(unsigned int n) {

    unsigned long long result = 0;
    Matrix* base = init_matrix(1,1,1,0);
    if (n == 0) {
        result = 0;
    }
    else {
        base = matrix_expo(base,(n-1));
        result = (base->a_00);
    }
    return result;      
}

But this code results in segmentation fault for values >= 3. I am feeling that the problem is with the mult_matrix() function but I am unable to find it. Please help.

EDIT : The bug causing the overflow has been found but there seems to be other logical errors here as well. In particular the edited code outputs wrong values and the outputs to consecutive inputs like 1 and 2 , 3 and 4, 5 and 6 are same. :(

The edited code :

#include <cmath>
#include <iostream>

using namespace std;

typedef struct Matrix {
    unsigned long long a_00;
    unsigned long long a_01;
    unsigned long long a_10;
    unsigned long long a_11;
}Matrix;

Matrix* init_matrix(unsigned long long n);
Matrix* init_matrix(unsigned long long a_00, unsigned long long a_01,
                    unsigned long long a_10, unsigned long long a_11);
Matrix* mult_matrix(Matrix* A, Matrix* B);      
Matrix* expo_matrix(Matrix* A,unsigned int n);
unsigned long long fib_matrix_expo(unsigned int n); 

int main() {

    unsigned int n;
    cin >> n;

    unsigned long long result = fib_matrix_expo(n);
    cout << result << endl;
    return 0;
}


Matrix* init_matrix(unsigned long long n) {
    Matrix* ret_matrix = (Matrix*)::operator new (sizeof(Matrix));

    (ret_matrix->a_00) = n;
    (ret_matrix->a_01) = n;
    (ret_matrix->a_10) = n;
    (ret_matrix->a_11) = n;

    return ret_matrix;
}

Matrix* init_matrix(unsigned long long a_00, unsigned long long a_01,
                    unsigned long long a_10, unsigned long long a_11) {
/*==================================================================================*/
    Matrix* ret_matrix = (Matrix*)::operator new (sizeof(Matrix));

    (ret_matrix->a_00) = a_00;
    (ret_matrix->a_01) = a_01;
    (ret_matrix->a_10) = a_10;
    (ret_matrix->a_11) = a_11;

    return ret_matrix;  

}


Matrix* matrix_mult(Matrix* A, Matrix* B) {
    Matrix* result = (Matrix*)::operator new (sizeof(Matrix));

    (result->a_00) = ((A->a_00) * (B->a_00) + (A->a_01)*(B->a_10));
    (result->a_01) = ((A->a_00) * (B->a_01) + (A->a_01)*(B->a_11));
    (result->a_10) = ((A->a_10) * (B->a_00) + (A->a_11)*(B->a_10));
    (result->a_11) = ((A->a_10) * (B->a_01) + (A->a_11)*(B->a_11));

    return result;
}

Matrix* matrix_expo(Matrix* base_matrix, unsigned int index) {
    Matrix* result;

    if (index == 0) {
        result = init_matrix(1);
    }
    else if (index == 1) {
        result = base_matrix;
    }
    else {
        Matrix* temp = matrix_expo(base_matrix,static_cast<unsigned int>(index/2));
        result = matrix_mult(temp,temp);
        result = matrix_mult(result,matrix_expo(base_matrix,(index % 2)));
    }
    return result;
}

unsigned long long fib_matrix_expo(unsigned int n) {

    unsigned long long result = 0;
    Matrix* base = init_matrix(1,1,1,0);
    if (n == 0) {
        result = 0;
    }
    else {
        base = matrix_expo(base,(n-1));
        result = (base->a_00);
    }
    return result;      
}

EDIT : This bug has also been resolved.But as pointed out by Paul this has lots of memory leaks. I shall try to implement this without using manual dynamic memory allocation and repost it. Then as advised by n.m. I shall update the code after using constructors and removing other unnecessary complex features.The current working code:

EDIT:The code free of memory leaks:

#include <cmath>
#include <iostream>

using namespace std;

struct Matrix {
    unsigned long long a_00;
    unsigned long long a_01;
    unsigned long long a_10;
    unsigned long long a_11;
};

Matrix init_matrix(unsigned long long n);
Matrix init_matrix(unsigned long long a_00, unsigned long long a_01,
                   unsigned long long a_10, unsigned long long a_11);
Matrix mult_matrix(const Matrix& A, const Matrix& B);      
Matrix expo_matrix(const Matrix& A,unsigned int n);
unsigned long long fib_matrix_expo(unsigned int n); 

int main() {

    unsigned int n;
    cin >> n;

    unsigned long long result = fib_matrix_expo(n);
    cout << result << endl;
    return 0;
}


Matrix init_matrix(unsigned long long n) {
    Matrix ret_matrix;

    ret_matrix.a_00 = n;
    ret_matrix.a_01 = n;
    ret_matrix.a_10 = n;
    ret_matrix.a_11 = n;

    return ret_matrix;
}

Matrix init_matrix(unsigned long long a_00, unsigned long long a_01,
                    unsigned long long a_10, unsigned long long a_11) {
/*==================================================================================*/
    Matrix ret_matrix;

    ret_matrix.a_00 = a_00;
    ret_matrix.a_01 = a_01;
    ret_matrix.a_10 = a_10;
    ret_matrix.a_11 = a_11;
    return ret_matrix;  
}


Matrix matrix_mult(const Matrix& A, const Matrix& B) {
    Matrix result;

    result.a_00 = ((A.a_00) * (B.a_00) + (A.a_01)*(B.a_10));
    result.a_01 = ((A.a_00) * (B.a_01) + (A.a_01)*(B.a_11));
    result.a_10 = ((A.a_10) * (B.a_00) + (A.a_11)*(B.a_10));
    result.a_11 = ((A.a_10) * (B.a_01) + (A.a_11)*(B.a_11));

    return result;
}

Matrix matrix_expo(const Matrix& base_matrix, unsigned int index) {
    Matrix result;

    if (index == 0) {
        result = init_matrix(1);
    }
    else if (index == 1) {
        result = base_matrix;
    }
    else {
        Matrix temp = matrix_expo(base_matrix,static_cast<unsigned int>(index/2));
        result = matrix_mult(temp,temp);
        if(index % 2 == 1) {
            result = matrix_mult(result,base_matrix);
        }
    }
    return result;
}

unsigned long long fib_matrix_expo(unsigned int n) {

    unsigned long long result = 0;
    Matrix base = init_matrix(1,1,1,0);
    if (n == 0) {
        result = 0;
    }
    else {
        base = matrix_expo(base,(n-1));
        result = base.a_00;
    }
    return result;      
}

Answer 1

Here is a C++ implementation using a class for Matrix. Unfortunately, I don't know enough of matrices to point out "the" problem with your original code; simply rewriting it with a class seems to have gotten rid of the problem. I couldn't interpret your expo function inner workings so I wrote it as a simple loop instead.

Using a class means that all operations are "on" a single Matrix only; no need to make multiple copies anymore! (Except for the 1,1,1,0 helper in the expo function; C++ will automatically delete the temp1 Matrix at the end of that function.)

A new Matrix gets initialized by default with 1,1,1,0, so the init in expo and fib_matrix_expo can be removed.

This also means you can easily expand the Matrix class to perform additional functions.

#include <cmath>
#include <iostream>

using namespace std;

class Matrix {
    unsigned long long a_00;
    unsigned long long a_01;
    unsigned long long a_10;
    unsigned long long a_11;

public:
    Matrix();
    void init(unsigned long long n);
    void mult(Matrix* B);
    void expo (unsigned int n);
    unsigned long long value ()
    {
        return a_00;
    }
};

unsigned long long fib_matrix_expo(unsigned int n); 

int main() {

    unsigned int n;
    cin >> n;

    unsigned long long result = fib_matrix_expo(n);
    cout << result << endl;
    return 0;
}

/* In this, and all other class functions, the specifier 'this->' is
   not required. Per definition, all variables from the class are
   immediately accessible. However, including 'this->..' makes it
   clear you *are* modifying the class variables, and not local ones
   with the same name */
Matrix::Matrix () {
    this->a_00 = 1;
    this->a_01 = 1;
    this->a_10 = 1;
    this->a_11 = 0;
}


void Matrix::init (unsigned long long n) {
    this->a_00 = n;
    this->a_01 = n;
    this->a_10 = n;
    this->a_11 = n;
}

void Matrix::mult(Matrix* B)
{
    unsigned long long t00,t01,t10,t11;

    t00 = this->a_00 * B->a_00 + this->a_01 * B->a_10;
    t01 = this->a_00 * B->a_01 + this->a_01 * B->a_11;
    t10 = this->a_10 * B->a_00 + this->a_11 * B->a_10;
    t11 = this->a_10 * B->a_01 + this->a_11 * B->a_11;

    this->a_00 = t00;
    this->a_01 = t01;
    this->a_10 = t10;
    this->a_11 = t11;

    // printf ("%llu %llu %llu %llu\n", a_00, a_01, a_10, a_11);
}

void Matrix::expo(unsigned int index)
{
    if (index == 0)
    {
        this->init(1);
        return;
    }
    if (index == 1)
    {
        return;
    }

    Matrix temp1;
    // temp1 will initialize itself with (1,1,1,0),
    // nothing more needed
    while (--index > 0)
    {
        this->mult (&temp1);
    }
}

unsigned long long fib_matrix_expo (unsigned int n)
{
    Matrix base;
    if (n == 0)
        return 0;

    // No need to initialize to the default values, as this
    // is done in the default constructor.
    base.expo(n-1);
    return base.value();
    // Since 'base' is a locally created class object, you
    // do not need to delete it.
}

^{Constructive criticism is welcomed, as I'm not a daily C++ programmer ;)}

Answer 2

You have a stack overflow in the matrix_expo function because of recursion here

Matrix* temp = matrix_expo(base_matrix,static_cast<unsigned int>(index/2));
result = matrix_mult( matrix_expo(temp,2), matrix_expo(base_matrix,(index % 2)) );

lets see why

Matrix matrix_expo(Matrix base_matrix, unsigned int index) {
    Matrix result;

    if (index == 0) {
        result = init_matrix(1);
    }
    else if (index == 1) {
        result = base_matrix;
    }
    else {
        Matrix temp = matrix_expo(base_matrix,static_cast<unsigned int>(index/2));
        result = matrix_mult( matrix_expo(temp,2), matrix_expo(base_matrix,(index % 2)) );
    }
    return result;
}

if index == 2 then matrix_expo(temp, 2) is always reached, and that is causing infinite recursion.

A few other problems are with new and never using delete and here actually loosing the chance to delete

base = matrix_expo(base,(n-1));

since base could be allocated using new and you reassing the base pointer which was pointing to new allocated memory.

I say could be allocated using new because of this line

result = base_matrix;

here you are returning base it self in the case if the previously mentioned case, but in any other case, the other assignments to result in this function will allocate using new this one will return a pointer to the originally allocated memory, hence if you use the delete operator to free this memory, you have a chance to double free.