_EFFICIENTLY RAISING MATRICES TO AN INTEGER POWER_
by Victor Duvanenko

[LISTING ONE]

#include
#include
#include

#define N 2

/* Procedure to multiply two square matrices */
void matrix_mult( mat1, mat2, result, n )
double *mat1, *mat2; /* pointers to the matrices to be multiplied */
double *result; /* pointer to the result matrix */
int n; /* n x n matrices */
{
register i, j, k;
for( i = 0; i < n; i++ )
for( j = 0; j < n; j++, result++ ) {
*result = 0.0;
for( k = 0; k < n; k++ )
*result += *( mat1 + i * n + k ) * *( mat2 + k * n + j );
/* result[i][j] += mat1[i][k] * mat2[k][j]; */
}
}
/* Procedure to copy square matrices quickly. Assumes the elements are
double precision floating-point type. */
void matrix_copy( src, dest, n )
double *src, *dest;
int n;
{
memcpy( (void *)dest, (void *)src, (unsigned)( n * n * sizeof( double )));
}
/* Procedure to compute Fibonacci numbers in linear time */
double fibonacci( n )
int n;
{
register i;
double fib_n, fib_n1, fib_n2;

if ( n <= 0 ) return( 0.0 );
if ( n == 1 ) return( 1.0 );
fib_n2 = 0.0; /* initial value of Fibonacci( n - 2 ) */
fib_n1 = 1.0; /* initial value of Fibonacci( n - 1 ) */
for( i = 2; i <= n; i++ ) {
fib_n = fib_n1 + fib_n2;
fib_n2 = fib_n1;
fib_n1 = fib_n;
}
return( fib_n );
}
/* Procedure to compute Fibonacci numbers recursively (inefficiently) */
double fibonacci_recursive( n )
int n;
{
if ( n <= 0 ) return( 0.0 );
if ( n == 1 ) return( 1.0 );
return( (double)( fibonacci_recursive(n - 1) + fibonacci_recursive(n - 2)));
}
/* Procedure to compute Fibonacci numbers in logarithmic time (fast!) */
double fibonacci_log( n )
int n;
{
register k, bit, num_bits;
double a[2][2], b[2][2], c[2][2], d[2][2];
if ( n <= 0 ) return( 0.0 );
if ( n == 1 ) return( 1.0 );
if ( n == 2 ) return( 1.0 );

n--; /* need only a^(n-1) */
a[0][0] = 1.0; a[0][1] = 1.0; /* initialize the Fibonacci matrix */
a[1][0] = 1.0; a[1][1] = 0.0;
c[0][0] = 1.0; c[0][1] = 0.0; /* initialize the result to identity */
c[1][0] = 0.0; c[1][1] = 1.0;

/* need to convert log bases as only log base e (or ln) is available */
num_bits = ceil( log((double)( n + 1 )) / log( 2.0 ));

/* Result will be in matrix 'c'. Result (c) == a if bit0 is 1. */
bit = 1;
if ( n & bit ) matrix_copy( a, c, N );
for( bit <<= 1, k = 1; k < num_bits; k++ ) /* Do bit1 through num_bits. */
{
matrix_mult( a, a, b, N ); /* square matrix a; result in matrix b */
if ( n & bit ) { /* adjust the result */
matrix_mult( b, c, d, N );
matrix_copy( d, c, N );
}
matrix_copy( b, a, N );
bit <<= 1; /* next bit */
}
return( c[0][0] );
}
main()
{
int n;
for(;;) {
printf( "Enter the Fibonacci number to compute ( 0 to exit ): " );
scanf( "%d", &n );
if ( n == 0 ) break;
printf("\nMatrix method: Fibonacci( %d ) = %le\n",n,fibonacci_log(n));
printf("\nLinear method: Fibonacci( %d ) = %le\n",n,fibonacci(n));
printf("\nRecursive method: Fibonacci( %d ) = %le\n",n,fibonacci_recursive(n));
}
return(0);
}

