www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - matrix and fibonacci

reply "newcomer[bob]" <bob dontknow.com> writes:
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
next sibling parent reply "newcomer[bob]" <bob dontknow.com> writes:
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
parent reply "newcomer[bob]" <bob dontknow.com> writes:
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
Mar 09 2012
next sibling parent sclytrack <sclytrack hotmail.com> writes:
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
prev sibling parent Matthias Walter <xammy xammy.homelinux.net> writes:
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
 


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
prev sibling parent reply Timon Gehr <timon.gehr gmx.ch> writes:
On 03/09/2012 06:50 AM, 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]


// 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
parent reply "newcomer[bob]" <bob dontknow.com> writes:
On Friday, 9 March 2012 at 14:07:10 UTC, Timon Gehr wrote:

 On 03/09/2012 06:50 AM, 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]


 // 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
parent Timon Gehr <timon.gehr gmx.ch> writes:
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:

 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]




 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