_EFFICIENTLY RAISING MATRICES TO AN INTEGER POWER_
by Victor Duvanenko

[LISTING ONE]

#include <stdio.h>
#include <math.h>
#include <string.h>

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