@<font color="red”> 本文为"Algorithms"一书的笔记 @</font>

fibonacci

书上给出了另外一种我想不到的解法:

#include <stdio.h>

#define MAXN 100

int fib[MAXN];

int get_fib(int num);

int main(int argc, char *argv[])
{
    printf("%d\n", get_fib(10));

    return 0;
}

int get_fib(int num)
{
    int i;

    fib[0] = 0;
    fib[1] = 1;

    for (i = 2; i <= num; i++) {
        fib[i] = fib[i - 1] + fib[i - 2];
    }

    return fib[num];
}

我以前写过递归的,记忆性递归的,迭代的,但是这种以循环代替递归的还没写过, 也没见过,虽然理解起来非常简单,但是为什么我就想不到呢? 不过这个程序,如果要求多个数的fibonacci值,还要在循环里面判断一下:

if (fib[i] != 0) {
    fib[i] = fib[i - 1] + fib[i - 2];
}

这样可以避免每次计算,以空间换取时间。

Big-oh notation

每一本算法书都会讲这个东西,离散数学也讲这东西,基本讲一次忘一次,这次算是 回忆,却让我理解更透彻了。 Ο(n),英文读作:big oh of n 定义如下:


Let f(n) and g(n) be functions from positive integers to positive reals.

We say f = Ο(g) if there is a constant c > 0 such that f(n) <= c * g(n).


f = Ο(g)表明f增长得没有g快,是f <= g的一个简单模拟。 比如:


f1(n) = n ^ 2, f2(n) = 2 * n + 20,

f2 = Ο(f1), 因为:

f2(n) / f1(n) = (2 * n + 20) / n ^ 2 <= 22

可以取1到22之间的任意一个常数。

而f1 != Ο(f2),因为:

f1(n) / f2(n) = n ^ 2/ (2 * n + 20) 趋于无穷


有关Ο 的一些常识:

  1. 常数相乘可以忽略,如:14 * n ^ 2 –> n ^ 2
  2. 如果有n ^ a和n ^ b,且a > b,忽略n ^ b,如:n ^ 2 + n –> n ^ 2
  3. 如果有指数和多项式,则忽略多项式,如:n ^ 3 + n ^ 5 –> n ^ 3
  4. 如果有多项式和对数,则忽略对数,如:n + (logn) ^ 3 –> n,再如:

n ^ 2 + nlogn –> n ^ 2 作为扩展,有如下定义:


如果f = Ο(g),则g = Ω(f)

如果f = Ο(g), 且f = Ω(g),则f = Θ(g)


more fibonacci algorithms

第一种方法: 用矩阵来求

我们知道:

0 * bullet F0 + F1 = F1
F0 + F1 = F2

根据线性方程组转矩阵的方法,可以变成:

\begin{equation} \left( \begin{array}{cc} 0 & 1 \
1 & 1 \end{array} \right) \bullet \left( \begin{array}{c} F0 \
F1 \end{array} \right)

\left( \begin{array}{c} F1 \
F2 \end{array} \right) \end{equation}

通过数学归纳法,可以推出:

\begin{equation} \left( \begin{array}{cc} 0 & 1 \
1 & 1 \end{array} \right)^{n} \bullet \left( \begin{array}{c} F0 \
F1 \end{array} \right)

\left( \begin{array}{c} F_{n} \
F_{n+1} \end{array} \right) \end{equation}

于是,可以写代码来实现这个算法了。

main.c

#include "matrix.h"

/* return the fibonacci number of NUM */
float fibonacci(int num);

int main(int argc, char *argv[])
{
    printf("%g\n", fibonacci(10));

    return 0;
}

float fibonacci(int num)
{
    Matrix *m1, *m2, *exp;
    float **arr, **vector, result;
    int i, j;

    arr = (float **)malloc(2 * sizeof(float *));
    for (i = 0; i < 2; i++) {
        arr[i] = (float *)malloc(2 * sizeof(float));
        for (j = 0; j < 2; j++) {
            arr[i][j] = 1;
        }
    }
    arr[0][0] = 0;

    vector = (float **)malloc(2 * sizeof(float *));
    for (i = 0; i < 2; i++) {
        vector[i] = (float *)malloc(sizeof(float));
    }
    vector[0][0] = 0;
    vector[1][0] = 1;

    m1 = matrix_from_array(arr, 2, 2);
    m2 = matrix_from_array(vector, 2, 1);
    exp = matrix_from_file(NULL, 2, 2);

    matrix_power(m1, num, exp);
    matrix_multiple(exp, m2, m2);
    result = *(m2->elem);

    matrix_delete(m1);
    matrix_delete(m2);
    matrix_delete(exp);
    for (i = 0; i < 2; i++) {
        free(arr[i]);
    }
    free(arr);
    for (i = 0; i < 2; i++) {
        free(vector[i]);
    }
    free(vector);

    return result;
}

matrix.h

#ifndef _MATRIX_H_
#define _MATRIX_H_

#include <stdio.h>
#include <stdlib.h>

//matrix element's type
typedef double matrix_t;
typedef struct matrix Matrix;

struct matrix {
    int row_no;       //row number
    int col_no;       //column number
    matrix_t *elem;   //point to the elements
};

/* create a matrix from a file pointed by FP, with
 * ROW_NO rows and COL_NO columns, return the pointer
 * to the matrix.
 * [note]:if FP is NULL, return an empty matrix */
Matrix *
matrix_from_file(FILE *fp, int row_no, int col_no);

/* create a matrix from an two dimension pointer ARR with
 * ROW_NO rows and COL_NO columns, return the pointer
 * to the matrix.
 * [note]:if ARR is NULL, return an empty matrix */
Matrix *
matrix_from_array(float **arr, int row_no, int col_no);

/* free all memory */
void matrix_delete(Matrix *m);

/* initial the matrix M as identity matrix */
void matrix_identity(Matrix *m);

/* multiple matrix LEFT and RIGHT, the result is matrix RESULT.
 * with a stupid algorithm. */
void matrix_multiple(Matrix *left, Matrix *right, Matrix *result);

/* calculate the matrix's power of exponent EXP, the result is
 * the matrix RESULT. */   
void matrix_power(Matrix *m, int exp, Matrix *result);

#endif /* _MATRIX_H_ */

matrix.c

#include "matrix.h"

static void assert_matrix(const Matrix *);
static int can_mul(const Matrix *left, const Matrix *right,
                   const Matrix *result);
static int is_square(const Matrix *m);

Matrix *
matrix_from_file(FILE *fp, int row_no, int col_no)
{
    int i;
    Matrix *m = (Matrix *)malloc(sizeof(Matrix));

    if (m == NULL) {
        fprintf(stderr, "memory no enough!");
        exit(1);
    }
    m->row_no = row_no;
    m->col_no = col_no;
    m->elem = (matrix_t *)malloc(row_no * col_no * sizeof(matrix_t));
    if (m->elem == NULL) {
        fprintf(stderr, "memory no enough!");
        exit(1);
    }
    if (fp != NULL) {
        for (i = 0; i < row_no * col_no; i++) {
            fscanf(fp, "%lf", (m->elem + i));
        }
    }

    return m;
}

Matrix *
matrix_from_array(float **arr, int row_no, int col_no)
{
    int i, j;
    Matrix *m = (Matrix *)malloc(sizeof(Matrix));

    if (m == NULL) {
        fprintf(stderr, "memory no enough!");
        exit(1);
    }
    m->row_no = row_no;
    m->col_no = col_no;
    m->elem = (matrix_t *)malloc(row_no * col_no * sizeof(matrix_t));
    if (m->elem == NULL) {
        fprintf(stderr, "memory no enough!");
        exit(1);
    }
    if (arr != NULL) {
        for (i = 0; i < row_no; i++) {
            for (j = 0; j < col_no; j++) {
                *(m->elem + i * col_no + j) = arr[i][j];
            }
        }
    }

    return m;
}

void matrix_delete(Matrix *m)
{
    assert_matrix(m);

    if (m->elem) {
        free(m->elem);
    }
    free(m);
}

void matrix_identity(Matrix *m)
{
    int i;

    assert_matrix(m);

    if (!is_square(m)) {
        fprintf(stderr, "identity: the matrix is not square!");
        exit(1);
    }
    for (i = 0; i < m->row_no * m->col_no; i++) {
        if (((i / m->col_no) * m->row_no) + (i / m->col_no) == i) {
            *(m->elem + i) = 1;
        } else {
            *(m->elem + i) = 0;
        }
    }
}

void matrix_multiple(Matrix *left, Matrix *right, Matrix *result)
{
    int i, j, k, left_elem, right_elem;
    matrix_t sum;
    Matrix *temp;

    assert_matrix(left);
    assert_matrix(right);
    assert_matrix(result);

    if (!can_mul(left, right, result)) {
        fprintf(stderr, "multiple: row number or column number not match!");
        exit(1);
    }
    if (result == left || result == right) {
        temp = matrix_from_file(NULL, result->row_no, result->col_no);
        for (i = 0; i < result->row_no * result->col_no; i++) {
            *(temp->elem + i) = *(result->elem + i);
        }
    } else {
        temp = result;
    }
    for (i = 0; i < left->row_no; i++) {
        for (j = 0; j < right->col_no; j++) {
            sum = 0.0;
            for (k = 0; k < left->col_no; k++) {
                left_elem = *(left->elem + i * left->col_no + k);
                right_elem = *(right->elem + j * right->col_no + k);
                if ((left_elem != 0) && (right_elem != 0)) {
                    sum += left_elem * right_elem;
                }
            }
            *(temp->elem + i * temp->col_no + j) = sum;
        }
    }
    if (result == left || result == right) {
        for (i = 0; i < result->row_no * result->col_no; i++) {
            *(result->elem + i) = *(temp->elem + i);
        }
    }
}

void matrix_power(Matrix *m, int exp, Matrix *result)
{
    matrix_identity(result);

    while (exp > 0) {
        if ((exp & 1) == 1) {
            matrix_multiple(m, result, result);
        }
        matrix_multiple(m, m, m);
        exp >>= 1;
    }
}

static void assert_matrix(const Matrix *m)
{
    if (m == NULL) {
        fprintf(stderr, "the matrix is NULL!");
    }
}

static int can_mul(const Matrix *left, const Matrix *right,
                   const Matrix *result)
{
    if (left->col_no == right->row_no &&
        right->row_no == result->row_no &&
        right->col_no == result->col_no) {
        return 1;
    } else {
        return 0;
    }
}

static int is_square(const Matrix *m)
{
    return m->row_no == m->col_no;
}

完整的矩阵类,请参见这里

第二种方法:称为fast doubling,引用自这里

从上面的矩阵法中,可以推出下面这个等式:

\begin{equation} \left( \begin{array}{cc} 0 & 1 \
1 & 1 \end{array} \right)^n

\left( \begin{array}{cc} F(n-1) & F(n) \
F(n) & F(n+1) \end{array} \right)^n \end{equation}

这个等式是可以通过数学归纳法证明的,也比较容易记忆(虽然我很不喜欢记公式, 但是有的东西,能记住是非常好的)。

于是,我们可以根据这个等式推出fast doubling的算法:

\begin{equation} \begin{split} \left( \begin{array}{cc} F(2n-1) & F(2n) \
F(2n) & F(2n+1) \end{array} \right)

\left( \begin{array}{cc} 0 & 1 \
1 & 1 \end{array} \right)^{2n}

\left( \left( \begin{array}{cc} 0 & 1 \
1 & 1 \end{array} \right)^n \right)^2

\left( \begin{array}{cc} F(n-1) & F(n) \
F(n) & F(n+1) \end{array} \right)^2 \

\left( \begin{array}{cc} F(n-1)^2 + F(n)^2 & F(n-1)F(n) + F(n)F(n+1) \
F(n)F(n-1) + F(n+1)F(n) & F(n)^2 + F(n+1)^2 \end{array} \right) \end{split} \end{equation}

于是有:

F(2n+1) = F(n)^2^ + F(n+1)^2^

F(2n) = F(n-1)F(n) + F(n)F(n+1) = F(n)( F(n-1) + F(n+1) ) = F(n)( F(n+1)

  • F(n) + F(n+1) ) = F(n)( 2F(n+1) - F(n) )

所以,如果知道F(k+1)和F(k),可以知道:

F(2k+1) = F(k)^2^ + F(k+1)^2^

F(2k) = F(k)( 2F(k+1) - F(k))

这样就得到了通用的方法。

有了算法,实现不了也是没用,马上敲代码吧!

#include <stdio.h>

/* calculate the fibonacci number of n with fast doubling algorithm. */
int fibonacci(int n);

int main(int argc, char *argv[])
{
    int i;

    for (i = 0; i < 20; i++) {
        printf("%d\n", fibonacci(i));
    }

    return 0;
}

int fibonacci(int n)
{
    int fn, fn1;

    if (n == 0) {
        return 0;
    }
    if (n == 1 || n == 2) {
        return 1;
    }
    if (n & 1) {
        fn = fibonacci((n - 1) / 2);
        fn1 = fibonacci((n - 1) / 2 + 1);
        return fn * fn + fn1 * fn1;
    } else {
        fn = fibonacci(n / 2);
        fn1 = fibonacci(n / 2 + 1);
        return fn * (2 * fn1 - fn);
    }
}