## digitalmars.D.learn - matrix and fibonacci

- "newcomer[bob]" <bob dontknow.com> Mar 08 2012
- "newcomer[bob]" <bob dontknow.com> Mar 09 2012
- sclytrack <sclytrack hotmail.com> Mar 09 2012
- "newcomer[bob]" <bob dontknow.com> Mar 09 2012
- Matthias Walter <xammy xammy.homelinux.net> Mar 09 2012
- Timon Gehr <timon.gehr gmx.ch> Mar 09 2012
- Timon Gehr <timon.gehr gmx.ch> Mar 10 2012
- "newcomer[bob]" <bob dontknow.com> Mar 10 2012

The following is a matrix implementation of the fibonacci algorithm: int fib(int n) { int M[2][2] = {{1,0},{0,1}} for (int i = 1; i < n; i++) M = M * {{1,1},{1,0}} return M[0][0]; } problem is I don't really understand how matrix multiplication works so i cannot translate it to an equivalent solution in D. Thank for you assistance in advance. newcomer[bob]

Mar 08 2012

On Friday, 9 March 2012 at 05:50:03 UTC, newcomer[bob] wrote:The following is a matrix implementation of the fibonacci algorithm: int fib(int n) { int M[2][2] = {{1,0},{0,1}} for (int i = 1; i < n; i++) M = M * {{1,1},{1,0}} return M[0][0]; } problem is I don't really understand how matrix multiplication works so i cannot translate it to an equivalent solution in D. Thank for you assistance in advance. newcomer[bob]

I attempted the following but it does not work: int fib(int n) ulong[2][2] M = [[1,0], [0,1]]; ulong[2][2] C = [[1,1], [1,0]]; foreach (i; 0 .. n) { M[0][0] = M[0][0] * C[0][0] + M[0][0] * C[0][1]; M[0][1] = M[0][1] * C[0][1] + M[0][1] * C[1][1]; M[1][0] = M[0][1] * C[0][0] + M[1][1] * C[0][1]; M[1][1] = M[1][1] * C[0][1] + M[1][1] * C[1][1]; } return M[0][0]; } any ideas what I'm doing wrong? Thanks, Bob

Mar 09 2012

On 03/09/2012 10:51 AM, newcomer[bob] wrote:On Friday, 9 March 2012 at 09:22:47 UTC, newcomer[bob] wrote:On Friday, 9 March 2012 at 05:50:03 UTC, newcomer[bob] wrote:The following is a matrix implementation of the fibonacci algorithm: int fib(int n) { int M[2][2] = {{1,0},{0,1}} for (int i = 1; i < n; i++) M = M * {{1,1},{1,0}} return M[0][0]; } problem is I don't really understand how matrix multiplication works so i cannot translate it to an equivalent solution in D. Thank for you assistance in advance. newcomer[bob]

I attempted the following but it does not work: int fib(int n) long[2][2] M = [[1,0], [0,1]]; ulong[2][2] C = [[1,1], [1,0]]; foreach (i; 0 .. n) { M[0][0] = M[0][0] * C[0][0] + M[0][0] * C[0][1]; M[0][1] = M[0][1] * C[0][1] + M[0][1] * C[1][1]; M[1][0] = M[0][1] * C[0][0] + M[1][1] * C[0][1]; M[1][1] = M[1][1] * C[0][1] + M[1][1] * C[1][1]; } return M[0][0]; } any ideas what I'm doing wrong? Thanks, Bob

Turns out that this this is an algorithm for calculating the powers of two up to 2^63. I still cannot find how to modify it to produce the fibonacci sequence though. Any advice would be appreciated. Thanks, Bob

http://en.wikipedia.org/wiki/Matrix_multiplication#Non-technical_example http://en.wikipedia.org/wiki/Fibonacci_number#Matrix_form I don't understand what you are trying to do above.

Mar 09 2012

On Friday, 9 March 2012 at 09:22:47 UTC, newcomer[bob] wrote:On Friday, 9 March 2012 at 05:50:03 UTC, newcomer[bob] wrote:

I attempted the following but it does not work: int fib(int n) long[2][2] M = [[1,0], [0,1]]; ulong[2][2] C = [[1,1], [1,0]]; foreach (i; 0 .. n) { M[0][0] = M[0][0] * C[0][0] + M[0][0] * C[0][1]; M[0][1] = M[0][1] * C[0][1] + M[0][1] * C[1][1]; M[1][0] = M[0][1] * C[0][0] + M[1][1] * C[0][1]; M[1][1] = M[1][1] * C[0][1] + M[1][1] * C[1][1]; } return M[0][0]; } any ideas what I'm doing wrong? Thanks, Bob

Turns out that this this is an algorithm for calculating the powers of two up to 2^63. I still cannot find how to modify it to produce the fibonacci sequence though. Any advice would be appreciated. Thanks, Bob

Mar 09 2012

On 03/09/2012 10:51 AM, newcomer[bob] wrote:On Friday, 9 March 2012 at 09:22:47 UTC, newcomer[bob] wrote:On Friday, 9 March 2012 at 05:50:03 UTC, newcomer[bob] wrote:

I attempted the following but it does not work: int fib(int n) long[2][2] M = [[1,0], [0,1]]; ulong[2][2] C = [[1,1], [1,0]]; foreach (i; 0 .. n) { M[0][0] = M[0][0] * C[0][0] + M[0][0] * C[0][1]; M[0][1] = M[0][1] * C[0][1] + M[0][1] * C[1][1]; M[1][0] = M[0][1] * C[0][0] + M[1][1] * C[0][1]; M[1][1] = M[1][1] * C[0][1] + M[1][1] * C[1][1]; } return M[0][0]; } any ideas what I'm doing wrong? Thanks, Bob

Turns out that this this is an algorithm for calculating the powers of two up to 2^63. I still cannot find how to modify it to produce the fibonacci sequence though. Any advice would be appreciated. Thanks, Bob

The idea is not that bad but it contains a bug: Once you modify M[0][1] in the second step of your loop, its changed content is plugged into the third step. You might simply save the contents of M in 4 additional variables and then use these to fill M again (with the result of M*C). In fact it suffices to store three values as all matrices M (and C) are symmetric, that is, M[0][1] == M[1][0]. But I recommend to do this *after* you made the current version work. Matthias

Mar 09 2012

On 03/09/2012 06:50 AM, newcomer[bob] wrote:

// simple 2x2 matrix type alias int[2][2] Mat; // 2x2 matrix multiplication Mat matmul(Mat a, Mat b){ return [[a[0][0]*b[0][0]+a[0][1]*b[1][0], a[0][0]*b[0][1]+a[0][1]*b[1][1]], [a[1][0]*b[0][0]+a[1][1]*b[1][0], a[0][1]*b[0][1]+a[1][1]*b[1][1]]]; } // implementation of your algorithm int fib(int n)in{assert(n>=0 && n<46);}body{ int M[2][2] = [[1,0],[0,1]]; int F[2][2] = [[1,1],[1,0]]; foreach (i; 0..n) M = matmul(F,M); return M[0][0]; } // faster way of computing matrix power int fi(int n)in{assert(n>=0 && n<46);}body{ Mat M = [[1,0],[0,1]]; Mat F = [[1,1],[1,0]]; for(;n;n>>=1){ if(n&1) M = matmul(F,M); F = matmul(F,F); } return M[0][0]; } // closed form derived from basis transform to eigenbasis int f(int n)in{assert(n>=0 && n<46);}body{ enum sqrt5=sqrt(5.0); return cast(int)((((1+sqrt5)/2)^^(n+1)-((1-sqrt5)/2)^^(n+1))/sqrt5); }

Mar 09 2012

On 03/10/2012 11:00 PM, newcomer[bob] wrote:On Friday, 9 March 2012 at 14:07:10 UTC, Timon Gehr wrote:On 03/09/2012 06:50 AM, newcomer[bob] wrote:

Thanks very much for the assist. All three of these methods though, seem to have a bug.

I just start the sequence from 1, because a quick glance at your implementation showed that it starts from 1. But it seems like actually it would produce the sequence 1,1,1,2,3,...fib(), fi() and f() all produce the same incorrect result which for lack of better word I'll call an "off by one" bug. A call to any of these functions with any integer value between 2 and 44 yields the return value of the next higher integer, while 45 overflows.

I don't observe any overflowing. But the assertions are a little too tight, 46 would be the first one that does not work.For example, 2 => 2, 10 => 89, and 15 => 987 which are the actual values for 3, 11 and 16. Thanks, Bob

Probably this is what you want then: // simple 2x2 matrix type alias int[2][2] Mat; // 2x2 matrix multiplication Mat matmul(Mat a, Mat b){ return [[a[0][0]*b[0][0]+a[0][1]*b[1][0], a[0][0]*b[0][1]+a[0][1]*b[1][1]], [a[1][0]*b[0][0]+a[1][1]*b[1][0], a[0][1]*b[0][1]+a[1][1]*b[1][1]]]; } // implementation of your algorithm int fib(int n)in{assert(n>=0 && n<47);}body{ if(!n) return 0; int M[2][2] = [[1,0],[0,1]]; int F[2][2] = [[1,1],[1,0]]; foreach (i; 1..n) M = matmul(F,M); return M[0][0]; } // faster way of computing matrix power int fi(int n)in{assert(n>=0 && n<47);}body{ if(!n) return 0; Mat M = [[1,0],[0,1]]; Mat F = [[1,1],[1,0]]; for(n--;n;n>>=1){ if(n&1) M = matmul(F,M); F = matmul(F,F); } return M[0][0]; } // closed form derived from basis transform to eigenbasis int f(int n)in{assert(n>=0 && n<47);}body{ enum sqrt5=sqrt(5.0); return cast(int)((((1+sqrt5)/2)^^n-((1-sqrt5)/2)^^n)/sqrt5); }

Mar 10 2012

On Friday, 9 March 2012 at 14:07:10 UTC, Timon Gehr wrote:On 03/09/2012 06:50 AM, newcomer[bob] wrote:

// simple 2x2 matrix type alias int[2][2] Mat; // 2x2 matrix multiplication Mat matmul(Mat a, Mat b){ return [[a[0][0]*b[0][0]+a[0][1]*b[1][0], a[0][0]*b[0][1]+a[0][1]*b[1][1]], [a[1][0]*b[0][0]+a[1][1]*b[1][0], a[0][1]*b[0][1]+a[1][1]*b[1][1]]]; } // implementation of your algorithm int fib(int n)in{assert(n>=0 && n<46);}body{ int M[2][2] = [[1,0],[0,1]]; int F[2][2] = [[1,1],[1,0]]; foreach (i; 0..n) M = matmul(F,M); return M[0][0]; } // faster way of computing matrix power int fi(int n)in{assert(n>=0 && n<46);}body{ Mat M = [[1,0],[0,1]]; Mat F = [[1,1],[1,0]]; for(;n;n>>=1){ if(n&1) M = matmul(F,M); F = matmul(F,F); } return M[0][0]; } // closed form derived from basis transform to eigenbasis int f(int n)in{assert(n>=0 && n<46);}body{ enum sqrt5=sqrt(5.0); return cast(int)((((1+sqrt5)/2)^^(n+1)-((1-sqrt5)/2)^^(n+1))/sqrt5); }

Thanks very much for the assist. All three of these methods though, seem to have a bug. fib(), fi() and f() all produce the same incorrect result which for lack of better word I'll call an "off by one" bug. A call to any of these functions with any integer value between 2 and 44 yields the return value of the next higher integer, while 45 overflows. For example, 2 => 2, 10 => 89, and 15 => 987 which are the actual values for 3, 11 and 16. Thanks, Bob

Mar 10 2012