matrix and fibonacci

[Timon Gehr]

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);
}