home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / magazine / drdobbs / 1991 / 06 / matrices.asc < prev    next >
Text File  |  1991-05-02  |  4KB  |  109 lines
  1. _EFFICIENTLY RAISING MATRICES TO AN INTEGER POWER_
  2. by Victor Duvanenko
  3.  
  4. [LISTING ONE]
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8. #include <string.h>
  9.  
  10. #define N        2
  11.  
  12. /* Procedure to multiply two square matrices */
  13. void matrix_mult( mat1, mat2, result, n )
  14. double *mat1, *mat2;   /* pointers to the matrices to be multiplied */
  15. double *result;        /* pointer to the result matrix */
  16. int  n;                /* n x n matrices */
  17. {
  18.    register  i, j, k;
  19.    for( i = 0; i < n; i++ )
  20.       for( j = 0; j < n; j++, result++ ) {
  21.          *result = 0.0;
  22.          for( k = 0; k < n; k++ )
  23.             *result += *( mat1 + i * n + k ) * *( mat2 + k * n + j );
  24.           /* result[i][j] += mat1[i][k] * mat2[k][j]; */
  25.       }
  26. }
  27. /* Procedure to copy square matrices quickly.  Assumes the elements are
  28.    double precision floating-point type. */
  29. void matrix_copy( src, dest, n )
  30. double *src, *dest;
  31. int  n;
  32. {
  33.    memcpy( (void *)dest, (void *)src, (unsigned)( n * n * sizeof( double )));
  34. }
  35. /* Procedure to compute Fibonacci numbers in linear time */
  36. double fibonacci( n )
  37. int n;
  38. {
  39.    register i;
  40.    double fib_n, fib_n1, fib_n2;
  41.  
  42.    if ( n <= 0 )   return( 0.0 );
  43.    if ( n == 1 )   return( 1.0 );
  44.    fib_n2 = 0.0;                   /* initial value of Fibonacci( n - 2 ) */
  45.    fib_n1 = 1.0;                   /* initial value of Fibonacci( n - 1 ) */
  46.    for( i = 2; i <= n; i++ ) {
  47.       fib_n = fib_n1 + fib_n2;
  48.       fib_n2 = fib_n1;
  49.       fib_n1 = fib_n;
  50.    }
  51.    return( fib_n );
  52. }
  53. /* Procedure to compute Fibonacci numbers recursively (inefficiently) */
  54. double fibonacci_recursive( n )
  55. int n;
  56. {
  57.    if ( n <= 0 )   return( 0.0 );
  58.    if ( n == 1 )   return( 1.0 );
  59.    return( (double)( fibonacci_recursive(n - 1) + fibonacci_recursive(n - 2)));
  60. }
  61. /* Procedure to compute Fibonacci numbers in logarithmic time (fast!) */
  62. double fibonacci_log( n )
  63. int n;
  64. {
  65.    register k, bit, num_bits;
  66.    double a[2][2], b[2][2], c[2][2], d[2][2];
  67.    if ( n <= 0 )   return( 0.0 );
  68.    if ( n == 1 )   return( 1.0 );
  69.    if ( n == 2 )   return( 1.0 );
  70.  
  71.    n--;                             /* need only a^(n-1) */
  72.    a[0][0] = 1.0;  a[0][1] = 1.0;   /* initialize the Fibonacci matrix */
  73.    a[1][0] = 1.0;  a[1][1] = 0.0;
  74.    c[0][0] = 1.0;  c[0][1] = 0.0;   /* initialize the result to identity */
  75.    c[1][0] = 0.0;  c[1][1] = 1.0;
  76.  
  77.    /* need to convert log bases as only log base e (or ln) is available */
  78.    num_bits = ceil( log((double)( n + 1 )) / log( 2.0 ));
  79.  
  80.    /* Result will be in matrix 'c'. Result (c) == a if bit0 is 1. */
  81.    bit = 1;
  82.    if ( n & bit )   matrix_copy( a, c, N );
  83.    for( bit <<= 1, k = 1; k < num_bits; k++ )   /* Do bit1 through num_bits. */
  84.    {
  85.       matrix_mult( a, a, b, N );    /* square matrix a; result in matrix b */
  86.       if ( n & bit ) {              /* adjust the result */
  87.          matrix_mult( b, c, d, N );
  88.          matrix_copy( d, c, N );
  89.       }
  90.       matrix_copy( b, a, N );
  91.       bit <<= 1;                    /* next bit */
  92.    }
  93.    return( c[0][0] );
  94. }
  95. main()
  96. {
  97.    int  n;
  98.    for(;;) {
  99.       printf( "Enter the Fibonacci number to compute ( 0 to exit ): " );
  100.       scanf( "%d", &n );
  101.       if ( n == 0 )   break;
  102.       printf("\nMatrix    method: Fibonacci( %d ) = %le\n",n,fibonacci_log(n));
  103.       printf("\nLinear    method: Fibonacci( %d ) = %le\n",n,fibonacci(n));
  104.       printf("\nRecursive method: Fibonacci( %d ) = %le\n",n,fibonacci_recursive(n));
  105.    }
  106.    return(0);
  107. }
  108.  
  109.