Fast Fibonacci algorithm

| Forums | Site Guide |

Although the text wraps on this page, if you copy and paste it into an editor it should come out with correct formatting

/*
 * Algorithm to compute the nth fibonacci number in a logarithmic number of steps
 * By Matt Squire - www.InsideReality.net
 *
 * Inspired by exercise 1.19, from the Structure and Interpretation of Computer Programs (2nd ed)
 *
 * This implementation can only calculate for up to n = 93 before an integer overflow occurs.
 * Moggers87 has produced another version using the GNU Multiple Precision library
 * Source code available from: http://www.m87.co.uk/fib.c
 *
 */

#include <stdio.h>

typedef long long int int_64;

int main(int argc, char* argv[])
{
	//Check input
	if(argc < 2)
	{
		printf("Usage: fibonacci n (computes nth fibonacci number)\n");
		return 1;
	}

	//Set up initial state
	int_64 a = 1;
	int_64 b = 0;
	int_64 p = 0;
	int_64 q = 1;
	int n = atoi(argv[1]);

	//Check input size
	if(n > 93 || n < 0)
	{
		printf("Input must be at least zero and no more than 93\n");
		return 1;
	}

	//Main loop - n >> 1 is a binary shift right operation, i.e. division by 2
	while((n >>= 1) > 0)
	{
		//Compute next p and q values using the mapping
		//	q <- q^2 + 2pq
		//	p <- p^2 + q^2
		int_64 next_q = q*q + 2*p*q;
		int_64 next_p = q*q + p*p;
		p = next_p;
		q = next_q;
		
		//If n is odd...
		if(n % 2 == 1)
		{
			//Compute next a and b values using the mapping
			//	a <- bq + aq + ap
			//	b <- bp + aq
			int_64 next_a = b*q + a*q + a*p;
			int_64 next_b = b*p + a*q;
			a = next_a;
			b = next_b;
		}

		//Debugging line
		//printf("(%d, %d, %d, %d, %d),  ", n, p, q, a, b);
	}

	//Output result
	n = atoi(argv[1]); //Get back original input
	//The loop will produce a pair of integers (a, b)
	//If n is even, the nth fibonacci number is b
	//Otherwise, it is a
	if(n % 2 == 0)
		printf("%llu\n", b);
	else
		printf("%llu\n", a);

	return 0;
}