HDU 3117 由n^k的前4位，推出（斐波那契）矩阵^k的前四位

http://acm.hdu.edu.cn/showproblem.php?pid=3117

Fibonacci Numbers

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)

Total Submission(s): 2570    Accepted Submission(s): 1073

Problem Description

The Fibonacci sequence is the sequence of numbers such that every element is equal to the sum of the two previous elements, except for the first two elements f0 and f1 which are respectively zero and one.

What is the numerical value of the nth Fibonacci number?

 

Input

For each test case, a line will contain an integer i between 0 and 10⁸ inclusively, for which you must compute the ith Fibonacci number fi. Fibonacci numbers get large pretty quickly, so whenever the answer has more than 8 digits, output only the first and last 4 digits of the answer, separating the two parts with an ellipsis (“...”).

There is no special way to denote the end of the of the input, simply stop when the standard input terminates (after the EOF).

 

Sample Input

  
  

   
   0
1
2
3
4
5
35
36
37
38
39
40
64
65
  
  

 

Sample Output

  
  

   
   0
1
1
2
3
5
9227465
14930352
24157817
39088169
63245986
1023...4155
1061...7723
1716...7565
  
  

题意：求出第n项斐波那契数，如果位数小于8直接输出，大于8的只输出前4位和后4位。

思路：后4位，很好做，矩阵快速幂就可以了，（至于矩阵快速幂的原理就不介绍了）

F2=2 , F1=1 , F0=1 , F(-1)=0 , F(-2)=1;按照斐波那契规律再往前推两项  (这样就避免了特判，或者是避免了n--了）

初始矩阵                    系数矩阵     第n个矩阵的左上角的数字代表第n项斐波那契数，第n个矩阵=初始矩阵*系数矩阵的^n;（都%10000）

F(-1)0   F(-2)=1                1  1 

0           0                          1  0 

前4位的方法，跟网上其他代码不一样，我按照自己的思路打的代码。

首先介绍一下，求n^k的前4位怎么求，n是数字，k是次方。首先我令10^p=n^k;则有p=logn^k=k*logn;（log都是以10为底的）

p求着很简单，有了这个p我就能算出n^k，现在double  x=(int)p;  y=p-x;  y是p的小数部分，x是p的整数部分。

p=x+y,   so:  10^x  *  10^y  =10^p  =n^k;(这个式子没毛病吧）  那么10^x代表的是什么，也就是1后面x个0而已。

现在我只需要知道n^k的前4位数字，因为10^y=  (n^k) / (10^x); 因为10^x是10的整数倍，(n^k) / (10^x)就只是n^k把小数点往左移x位，数字是不变的，所以只需要知道y即可知道n^k的前几位数字。

(int)(10^y *1000)  就表示n^k的前4位数字。

计算方法：

 double  p=k*log(n);

p=p-(int)p;

int ans=(int)(pow(10,p)*1000);

懂了求n^k的前4位，再说这道题。这道题的n不是数字了，n表示一个矩阵。

只要求出p即可;n现在是矩阵了，如果我令10^p=n,然后再快速幂也可以(可是这个式子明显是不对的，n是矩阵，p还是没法求），由这个式子再思考一下，我们能不能反过来，让n=10^p;（这个式子就有可能对了），n表示一个矩阵，矩阵里面有数字，我都表示成10的p次方不就行了（矩阵里面的数字都是次方数p），这时候用快速幂求出n^k，不就知道了p。

所以现在构造的矩阵n里的数字已经不能是表示数字的了，应该是该数字转换成10^p的（double）p了，所以矩阵对应的数字改一下即可。

原来的  初始矩阵      系数矩阵     

              0  1                 1  1 

              0  0                 1  0 

(因为10的任何次方也不可能为0，又因为斐波那契额数不会小于1大于0的，所以把原来的0统统写成负数)    (是负数则特判就行了)

变过来就成了

初始矩阵(double)      系数矩阵(double)

-1  0                             0  0     (0表示1及10^0)  （若是2的话就表示100及10^2）

-1  -1                            0  -1    （-1表示原来的数字是0，遇到-1则特判即可）

原来的矩阵乘法也会改变，数字乘法  会变成  次方的乘法。数字加法 会变成 次方的加法

100*100=10000,及2+2=4；（乘法）

100+100=200，及2+log（pow（10，2-2）+1）=(double)log2+2;（加法）

再举个加法例子

数字加法 10^5.5+10^3.6=10^3.6  *  (10^1.9+1)  

次方加法  5.5+3.6=3.6+log ( pow( 10 , 1.9 ) + 1 );跟上面的一一对应。

奉上我的AC代码：

#include <iostream>
#include <queue>
#include<cstring>
#include<cstdio>
#include <cmath>
#include<algorithm>
#define LL long long
#define mod 2008
using namespace std;
int f[100];
struct node2//计算前4位的矩阵
{
    double b[5][5];
};
node2 operator *(node2 a,node2 b)//跟矩阵乘法类似
{
    node2 c;
    for(int i=0; i<2; i++)
        for(int j=0; j<2; j++)
        {
            double ma,mb;
            int flag1=0,flag2=0;//标记，记录是否有-1出现。
            if(a.b[i][0]<0||b.b[0][j]<0)  flag1=1;
            else  ma=a.b[i][0]+b.b[0][j];

             if(a.b[i][1]<0||b.b[1][j]<0)  flag2=1;
            else   mb=a.b[i][1]+b.b[1][j];

             if(flag1&&flag2)   c.b[i][j]=-1.0;
             else if(flag1&&!flag2)  c.b[i][j]=mb;
             else if(!flag1&&flag2)  c.b[i][j]=ma;
             else
             {
                 double mc;
                 if(mb>=ma) swap(mb,ma);
                    mc=mb; ma-=mb;  mb=0.0;
                 c.b[i][j]=mc+log10(1+pow(10,ma));
             }
        }
        return c;
}

node2 ans2,a2;
void init2()//初始化
{
    memset(ans2.b,0,sizeof(ans2.b));
    memset(a2.b,0,sizeof(a2.b));
    ans2.b[0][0]=-1.0;
    ans2.b[1][0]=-1.0;ans2.b[1][1]=-1.0;
    a2.b[1][1]=-1.0;
}
void pow2(int s)//快速幂
{
    init2();
    while(s)
    {
        if(s%2) ans2=ans2 * a2;
        a2= a2 * a2;
        s=s>>1;
    }
}

struct node//计算后四位的矩阵
{
    int a[5][5];
};
node operator *(node a,node b)
{
    node c;
    memset(c.a,0,sizeof(c.a));
    for(int i=0; i<2; i++)
        for(int j=0; j<2; j++)
        {
            for(int k=0; k<2; k++)
                c.a[i][j]+=a.a[i][k]*b.a[k][j];
            c.a[i][j]%=10000;
        }
        return c;
}
node ans1,a1;
void init1()
{
    memset(ans1.a,0,sizeof(ans1.a));
    memset(a1.a,0,sizeof(a1.a));
    ans1.a[0][1]=1;
    a1.a[0][0]=1;a1.a[0][1]=1;
    a1.a[1][0]=1;
}
void pow1(int s)
{
    init1();
    while(s)
    {
        if(s%2) ans1=ans1*a1;
        a1=a1*a1;
        s=s>>1;
    }
}
int main()
{
    f[0]=0;
    f[1]=1;
    for(int i=2;i<=39;i++)
        f[i]=f[i-1]+f[i-2];
    int k;
    while(~scanf("%d",&k))
    {
        if(k<=39)  printf("%d\n",f[k]);
        else
        {
           pow2(k);
           double p=ans2.b[0][0];
           p=p-(int)p;
           p=pow(10,p);
           printf("%d...",(int)(p*1000));
           pow1(k);
           printf("%04d\n",ans1.a[0][0]);
        }
    }
}