此题对数学知识的要求比较高，9楼版主所说的循环，程序上的验证是正确的，从数学原理上的解释，我还没有弄明白，望能解惑。
下面是我从另一个角度的探索，从斐波那契的递推公式Fn=Fn-1+Fn-2，其中F1=F2=1可以得出以下两个关系，
F(n+m-1)=F(n-1)*F(m-1)+F(n)*F(m)
F(n+m)=F(n-1)*F(m)+F(n)*(F(m-1)+F(m))
这两个公式可以做到，从n-1项，n项，m-1项,m项这4项的值得到m+n-1项,m+n项的值，大概就是可以做到两项之和，这样在计算某一项时，就可以不用一项一项的去求，从某一项开始，直接加上一个合适的项就可以了。以下是代码实现：

#include<stdio.h>
#include<stdlib.h>
#define modN 10007
int main()
{
   unsigned int f0,f1 ,rs0,rs1;
   unsigned int i,n,m,end,temp0,temp1;
       
       while(1)
       {
          if(scanf("%d",&n)==EOF )break;    
          m=n;
          f0=0;     
           f1=1;
           rs0=1;
          rs1=0;
          i=1; 
          end=(m&0xffff0000)==0?16:32;//m<=0xffff0000时只计算至二进制数低16位  
          while(i<=end)
         {
            
                  if (m&1)
               {
                    //    F(n+m-1)=F(n-1)*F(m-1)+F(n)*F(m)
                       // F(n+m)=F(n-1)*F(m)+F(n)*(F(m-1)+F(m))
                       // 从rs0,rs1,f0,f1得到rs1+f1-1,rs1+f1 各项值的累加 
                    temp0=(rs0*f0+rs1*f1)%modN;
                    temp1=(rs0*f1+rs1*(f0+f1))%modN;
                    rs0=temp0;
                    rs1=temp1;
             
                }
                //    F(n+m-1)=F(n-1)*F(m-1)+F(n)*F(m)
                    //     F(n+m)=F(n-1)*F(m)+F(n)*(F(m-1)+F(m))
                    //    从f0,f1,f0,f1得到f0+f1-1,f0+f1 项数增加一倍，也就是二进制位向上一位 f1 始终记录2^i项的值 
                temp0=(f0*f0+f1*f1)%modN;
                    temp1=(f0*f1+f1*(f0+f1))%modN;
                    f0=temp0;
                    f1=temp1;    
                m/=2;
                i++; 
            }
            printf("%u\n",rs1);
       }
   
       return 1;
}

以前看到的算法

我写一个可以递归的

设 矩阵 M22 = [1, 1][1, 0]

[F(n), F(n-1)] = [F(n-1), F(n-2)] * M = [F(n-2), F(n-3)] * M * M ...... = [1, 0] * M ^ (n - 1)

分治求后半部分

#include <stdio.h>

#define N 2
#define MOD 10007

typedef struct _tag_Node
{
    int ans[N][N];
} Matrix;

Matrix sa = { { { 1, 0 }, { 0, 1 } } };          // F[1]
Matrix sb = { { { 1, 1 }, { 1, 0 } } };          // 元组

Matrix multiply(Matrix a, Matrix b)
{
    int i, j, k;

    Matrix s = { 0 };                      

    for (i = 0;i < N;i++)
    {
        for (j = 0;j < N;j++)
        {
            for (k = 0;k < N;k++)
            {
                s.ans[i][j] += a.ans[i][k] * b.ans[k][j];
                
                s.ans[i][j] %= MOD;
            }
        }
    }

    return s;
}

Matrix fun(int n)
{
    Matrix div, ans;
    if (n <= 1)
    {
        return n ? sb : sa;
    }

    div = fun(n / 2);
    ans = multiply(div, div);
    if (n % 2)
    {
        return multiply(ans, sb);
    }
    else
    {
        return ans;
    }
}

Matrix fibonacci(int n)
{                                           
    return multiply(sa, fun(n - 1));
}

void print(Matrix s)
{
    int i, j;

    for (i = 0;i < N;i++)
    {
        for (j = 0;j < N;j++)
        {
            printf(" %d", s.ans[i][j]);
        }
        puts("");
    }
    printf("---------------\n");
}

int main()
{
    print(fibonacci(22));
    return 0;
}

[此贴子已经被作者于2016-3-12 01:09编辑过]

用矩阵方法求斐波那契n项的值，网上有很多介绍，求取 矩阵{{1，1}，{1,0}}n次方所得矩阵中1行2列就是n项的值，矩阵计算有分治递归的，有快速幂方法，矩阵相乘大都使用循环的常用方法，实际上还是可以再优化的，根据斐波那契数列的特点，所得矩阵中的值是有如下关系的f(1,2) = f(2,1)  ,f(1,1) = f(2,1) + f(2,2)  这样，根据矩阵乘法规则可以写出f(1,2)和f(2,2）的计算式，两数之和就是f(1,1),从中也可以看出矩阵中的数据只有两个是有效的，如果使用快速幂法，保留2^n项的数据时，有两个就够了。按这个方法写了代码，这个程序与我以前写的很相似。

#include<stdio.h>
#include<stdlib.h>
#define modN 10007
struct fb
{
  unsigned int f1;//F(n)项的值
  unsigned int f0;//F(n-1)项的值            
};
int main()
{
   struct fb f[32],rs;
   unsigned int i,n,m,end,temp0,temp1;
   
   //生成斐波那契数位序为2^i 的数据  
   f[0].f0=0;     
   f[0].f1=1; 
   for(i=1;i<32;i++)
       { 
      f[i].f0=(f[i-1].f0*f[i-1].f0+f[i-1].f1*f[i-1].f1)%modN;
      f[i].f1=(f[i-1].f0*f[i-1].f1+f[i-1].f1*(f[i-1].f0+f[i-1].f1))%modN;
       }    
       while(1)
       {
          if(scanf("%d",&n)==EOF )break;    
          m=n;
       i=0;
       rs.f0=1;
       rs.f1=0;
           while(i<32)
          {
                if (m&1)
          {
             temp0=(rs.f0*f[i].f0+rs.f1*f[i].f1)%modN;
             temp1=(rs.f0*f[i].f1+rs.f1*(f[i].f0+f[i].f1))%modN;
             rs.f0=temp0;
             rs.f1=temp1;
          }
          m/=2;
          i++; 
           }
           printf("%u\n",rs.f1);
        }
        return 1;
}