By duckladydinh, history, 9 months ago, ,

Hi friends,

During the implementation of matrix fast exponentiation, I have encountered a segmentation fault in the recursive approach. Even if I have managed to fix it by changing it to a loop, but I really want to know what is wrong with my code. I wonder if you can help me.

This is my recursive implementation. I cannot understand what is wrong with it. With n = 1, it is fine, but if n is greater than 1, then segmentation fault.

const int N = 100;

struct Matrix {
llnum c[N+N][N+N];
int n = 0, m = 0;
Matrix(int _n, int _m) {
n = _n, m = _m;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < m; ++j) {
c[i][j] = 0;
}
}
}
void out() {
REP(i, n) {
REP(j, m) {
cout << c[i][j] << " \n"[j+1 == m];
}
}
cout << "\n";
}
};

Matrix fast_pow(Matrix A, int n) {
if (n == 0) return ident(A.m);
if (n == 1) return A;
Matrix E = fast_pow(A, (n/2));
if (n % 2 == 1) {
return mul(E, A);
}
return E;
}



Thank you.

•
• -5
•

 » 9 months ago, # |   0 Auto comment: topic has been updated by DuckLadyDinh (previous revision, new revision, compare).
 » 9 months ago, # |   0 Auto comment: topic has been updated by DuckLadyDinh (previous revision, new revision, compare).
 » 9 months ago, # |   0 Auto comment: topic has been updated by DuckLadyDinh (previous revision, new revision, compare).
 » 9 months ago, # |   +16 I think the third line should be Matrix E = fast_pow(mul(A, A), (n/2));
•  » » 9 months ago, # ^ |   0 Thank you, but that should just return wrong result, shouldn't it?Still, can you explain why it returns a segmentation fault with n = 2? Actually I discover that by passing A as a reference &A will solve it, but I don't understand why.
•  » » » 9 months ago, # ^ |   0 You should probably show your Matrix class definition (I'm particularly interested in its variables, copy constructor, and destructor if there's any explicit).
•  » » » » 9 months ago, # ^ | ← Rev. 2 →   0 I have updated (llnum == long double) and... there is no destructor
•  » » » » 9 months ago, # ^ |   0 And also, the following non-recursive code works with the structure: Matrix fast_pow(Matrix a, llint n) { Matrix ans = ident(a.m); while (n > 0) { if (n & 1) ans = mul(ans, a); a = mul(a, a), n >>= 1; } return ans; } 
•  » » » 9 months ago, # ^ | ← Rev. 2 →   +4 Just a guess: All the memory is being allocated from stack memory.Every Matrix object takes 320kb of memory. (Assuming 200 * 200 long longs)At every level of recursion, you create 2 or 5 new Matrix objects (remember, pass by value create a new copy. In mul you pass 2 objects and receive one back), so lets assume 1.5MB per recursion depth in the worst case.So even if your depth of recursion is 20, you will have allocated over 30MB memory, which may be higher than the default stack limit of whatever system you ran this on . (Windows 1MB, Linux 8MB according to this)You can overcome this by using dynamic memory allocation i.e malloc / new (that will fetch memory from the heap)PS: You have a bug in your code. You forgot to square the result. Add another 1MB per recursion level (2 objects are passed by value to mul(), one is returned).
•  » » » » 9 months ago, # ^ |   0 Thank you. To be honest I am still a bit unclear about it, since it gives exception even when n = 2.But I have just analyzed another problem and recognize that too deep recursion can really cause segmentation fault (by manually increasing the depth). So I think that is the problem. Thank you so much.
 » 9 months ago, # |   0 Auto comment: topic has been updated by DuckLadyDinh (previous revision, new revision, compare).