斐波那契序列 集锦 (转)

斐波那契序列 集锦 (转)
[定理1] 标准Fibonacci序列（即第0项为0，第1项为1的序列）当N大于1时，一定有f(N)和f(N-1)互质

其实，结合“互质”的定义，和一个很经典的算法就可以轻松证明
对，就是辗转相除法
互质的定义就是最大公约数为1

数学归纳法是很有用的证明方法，我们接下来这个定理用数学归纳法就很好证明：
[定理2]若i为奇数， f(i)*f(i)=f(i-1)*f(i+1)+1，否则f(i)*f(i)=f(i-1)*f(i+1)-1
对，这个定理用数学归纳法可以轻松证明，大家有兴趣可以自己尝试

[定理3] f(n)=f(i)*f(n-i-1)+f(i+1)*f(n-i)

f(n)=f(1)*f(n-2)+ f(2)*f(n-1)
    =f(2)*f(n-3)+ f(3)*f(n-2)
    =f(3)*f(n-4)+ f(4)*f(n-3)
看来没有错，证明方法就是这样

这个公式也可以用来计算较大的fibonacci数除以某个数的余数

设i=n/2 不过，为了保证计算能延续下去 需要每次保留三个值
这样，下一次计算就可以利用这三个值来求出两个值，再相加就可以得到第三个值
譬如，计算出f(5),f(6),f(7)，可以计算出f(11)、f(12)，然后推出f(13)
就是刚才洛奇的悲鸣(364738334)所提到的矩阵方法
我们知道我们若要简单计算f(n)，有一种方法就是先保存
a=f(0),b=f(1)，然后每次设：
a'=b b'=a+b

并用新的a'和b'来继续这一运算

如果大家熟悉利用“矩阵”这一工具的话，就知道，如果把a、b写成一个向量[a,b]，完成上述操作相当于乘以矩阵
0 1
1 1
也就是说，如果我们要求第100个fibonacci数，只需要将矩阵
[0,1]乘上
0 1
1 1
的一百次方，再取出第一项

因为我们知道，矩阵运算满足结合律，一次次右乘那个矩阵完全可以用乘上那个矩阵的N次方代替，更进一步，那个矩阵的N次方就是这样的形式：
f(n-1) f(n)
f(n) f(n+1)

而求矩阵的N次方，由于矩阵乘法满足结合律，所以我们可以用log(N)的算法求出——这个算法大家都会么？
一个是二分，一个是基于二进制的求幂

二分的原理：要求矩阵的N次方A(N)，设i=N/2若N%2==1， 则 A(N)=A(i)*A(i)*A(1)若N%2==0， 则 A(N)=A(i)*A(i)

基于二进制的原理：将N拆为二进制数，譬如13=1101那么 A^13= A^8 * A^4 * A^1 （这里^表示幂运算）

也就是说，由A^1开始，自乘得到A^2，然后自乘得到A^4，如果N对应位为1，则将这个结果乘到目标上去

这样的话，将所有乘法改为模乘，就可以得到一个较大Fibonacci数除以M的余数

若不用递归，其实类似

http://acm.pku.edu.cn/JudgeOnline/problem?id=3070
这里用的fib矩阵略有不同，是
f(n+1) f(n)
f(n) f(n-1)
但实际上可以验证效果是一样的

这题是要求求F(n)的最后四位数，所有乘法过程增加一个模10000的步骤即可，大家可以收藏稍候AC

关于矩阵我们告一段落，等下会回来继续探讨利用矩阵来解决复杂些的Fibonacci问题

http://acm.hdu.edu.cn/showproblem.php?pid=1568
我们来看这题，这题要求求出Fibonacci某项的前四位

当然，用矩阵也可以解决这道题——只要将乘法改为乘并保留前四位

我们采用double 保留整数部分四位 这题最好还是double吧

不过显然有更好的解法——如果我们知道Fibonacci序列的通项公式

F(n) = (((1+Sqrt(5))/2)^n - ((1-Sqrt(5))/2)^n)*1/Sqrt(5)

不过组合数学里也有这一公式的推导方法 叫做“线性齐次递推式”

这个解法的核心是，通解是某个数的幂 将f(n)=x^n代入递推方程，可以解出三个通解 0和 (1+sqrt(5))/2

通常把“0”称作平凡解，那么特解就是通解的某个线性组合

再代入f(0)=0 f(1)=1，就可以得出我们刚才的公式

不过通常情况下，我们只需要记住那个公式就可以了

提醒大家，记忆公式的时候千万别忘记了系数1/sqrt(5)

因为(1-sqrt(5))/2的绝对值小于1

所以当i较大的时候，往往可以忽略掉这一项
f(i)≈((1+Sqrt(5))/2)^n/sqrt(5);

所以，刚才列举出的HDOJ的1568，可以很简单的30以内直接求解，30以上采用这个公式，还是用log(N)求幂的算法求解
恩，就是公式的前半部分

http://acm.hdu.edu.cn/showproblem.php?pid=1021
或 http://acm.zju.edu.cn/show_problem.php?pid=2060
Fibonacci某项是否被3整除

[定理5] 标准Fibonacci序列对任意大于2的正整数的余数序列，必然是以“0 1”为循环节开头的序列

显然0、1是序列开头，也就是说序列开头就是循环节开头

循环长度的计算貌似是个比较难的问题，我一时还没有想到有效解法，不过，要说明的是，计算复杂度时，这个循环节长度应该按复杂度O(N^2)计算

恩，证明方法是利用同余定理、反证法，还有我们之前证明过的相邻项一定互质的定理，留给大家家庭作业

http://acm.hdu.edu.cn/showproblem.php?pid=1588
这是前天比赛关于Fibonacci的一道题，大家先看看题。
Description看后半部分就行了

现在告诉大家一种正确解法，然后大家就可以去搞定这道题向别人炫耀了

首先，我们将问题整理一下，就是对等差数列 ai=k*i+b，求所有的f(ai)之和除以M的余数

当0<=i<N

大家有没有想到，因为ai是等差数列，倘若f(ai)也是个等什么序列，那说不定就有公式求了

f(ai)显然不是等差数列，直接看上去也不是等比数列

但是如果把f(ai)换成我们刚才所说的Fibonacci矩阵呢？

是的，可是我们对矩阵是直接求幂即可，由于矩阵加法的性质，我们要求A^ai的右上角元素之和，只要求A^ai之和的右上角元素

就矩阵这个东西来说，完全可以看作一个等比数列，
首项是：A^b
公比是：A^k
项数是：N

呵呵，我们可以把问题进一步简化

因为矩阵的加法对乘法也符合分配律，我们提出一个A^b来，形成这样的式子：
A^b*( I + A^k + (A^k)^2 + .... + (A^k)^(N-1) )

A^b 和 A^k 显然都可以用我们之前说过的方法计算出来，这剩下一部分累加怎么解决呢

简单起见，设A^k=B
要求 G(N)=I + ... + B^(N-1)，设i=N/2
若N为偶数，G(N)=G(i)+G(i)*B^i
若N为奇数，G(N)=I+ G(i)*B + G(i) * (B^(i+1))

呵呵，这个方法就是比赛当时ACRush用的
而农夫用的则是更精妙的方法，大家可想知道

我们来设置这样一个矩阵
B I
O I
其中O是零矩阵，I是单位矩阵

将它乘方，得到
B^2 I+B
O   I
乘三方，得到
B^3 I+B+B^2
O   I
乘四方，得到
B^4 I+B+B^2+B^3
O   I

既然已经转换成矩阵的幂了，继续用我们的二分或者二进制法，直接求出幂就可以了

http://online-judge.uva.es/p/v110/11089.html
大家来读读这一题

Fibinary数是指没有相邻的两个1的二进制数。给N，求出第N大的Fibinary数

相对于二进制中每一位的值是2的幂，十进制中每一位的值是十的幂，
Fibonacci进制是每一位的值是对应Fibonacci数的一种计数系统。
     8 5 3 2 1
1     1
2     1 0
3     1 0 0
4     1 0 1
5     1 0 0 0
6     1 0 0 1
7     1 0 1 0
8     1 0 0 0 0
9     1 0 0 0 1
10   1 0 0 1 0
11   1 0 1 0 0
12   1 0 1 0 1
以上是前12个数字对应的十进制到Fibonacci进制的表格

Fibonacci的运算方法很奇怪。首先，它每一位上非0即1，而且不同于二进制的逢二进一或者十进制的逢十进一，它的进位方法是逢连续两个1，则进1

譬如
1010110==> 1011000 ==> 1100000==>10000000

在最低位有个特殊情况，最低位既可以逢2进1，也可以和次低位一起逢相邻进1
这种奇怪的进位方法，换句话描述就是，不存在两个连续的1
因为Fibonacci数其实也增长很快，int范围内好像只有46个，本题只需要用最简单的办法转换成Fibonacii进制即可
其中一题是
http://www.mydrs.org/program/down/ ahoi2004da y1.pdf
中的第二题，叫做数字迷阵
还有一题是PKU上的很出名 的取石子问题
http://acm.pku.edu.cn/JudgeOnline/problem?id=1067

这题相当复杂，大家可以自己思考，往Fibonacci上想，也可以阅读这里的论文：
http://episte.math.ntu.edu.tw/articles/mm/mm_03_2_02/index.html
http://acm.pku.edu.cn/JudgeOnline/problem?id=2967

另外这题 可以利用Fibonacci判断数据范围进行优化设计。不过貌似可以水过去，仅仅给大家提供个思路罢

关于Fibonacci和黄金分割，还有很多更高明的结论和定理，希望大家也继续讨论，将自己的知识和他人共享
http://episte.math.ntu.edu.tw/articles/mm/mm_02_4_10/index.html
中例3详细讲述了用生成函数来计算Fibonacci数公式的运算过程。 http://acm.hdu.edu.cn/showproblem.php?pid=1568
Fibonacci 求fibonacci前4位

http://acm.hdu.edu.cn/showproblem.php?pid=1588
Gauss Fibonacci
http://acm.pku.edu.cn/JudgeOnline/problem?id=1067
取石子问题
http://episte.math.ntu.edu.tw/articles/mm/mm_03_2_02/index.html 是报告
http://acm.pku.edu.cn/JudgeOnline/problem?id=3070
Fibonacci矩阵
http://acm.hdu.edu.cn/showproblem.php?pid=1021
或
http://acm.zju.edu.cn/show_problem.php?pid=2060
Fibonacci某项是否被3整除
http://acm.pku.edu.cn/JudgeOnline/problem?id=2116
Fibonacci进制计算
http://acm.pku.edu.cn/JudgeOnline/problem?id=2967
利用Fibonacci判断数据范围进行优化设计。
http://online-judge.uva.es/p/v110/11089.html
Fi-binary numbers，Fibonacci进制。
http://www.mydrs.org/program/down/ahoi2004day1.pdf
第二题 数字迷阵   这些，是今天涉及到的资料和网页

posted @ 2009-09-04 00:09 Knuth_档案 阅读(50) | 评论 (0) | 编辑

费马小定理 素数判定 蒙哥马利算法
约定：
x%y为x取模y，即x除以y所得的余数，当x<y时，x%y=x，所有取模的运算对象都为整数。
x^y表示x的y次方。
乘方运算的优先级高于乘除和取模，加减的优先级最低。
见到x^y/z这样，就先算乘方，再算除法。
A/B，称为A除以B，也称为B除A。
若A%B=0，即称为A可以被B整除，也称B可以整除A。
A*B表示A乘以B或称A乘B，B乘A，B乘以A……都TMD的一样，靠！

复习一下小学数学
公因数：两个不同的自然数A和B，若有自然数C可以整除A也可以整除B，那么C就是A和B的公因数。
公倍数：两个不同的自然数A和B，若有自然数C可以被A整除也可以被B整除，那么C就是A和B的公倍数。
互质数：两个不同的自然数，它们只有一个公因数1，则称它们互质。

费马是法国数学家，又译“费尔马”，此人巨牛，他的简介请看下面。不看不知道，一看吓一跳。

费马小定理：
有N为任意正整数，P为素数，且N不能被P整除（显然N和P互质），则有：
N^P%P=N(即：N的P次方除以P的余数是N)

但是我查了很多资料见到的公式都是这个样子：
(N^(P-1))%P=1

后来分析了一下，两个式子其实是一样的，可以互相变形得到，原式可化为：
(N^P-N)%P=0(即：N的P次方减N可以被P整除，因为由费马小定理知道N的P次方除以P的余数是N)

把N提出来一个，N^P就成了你N*(N^(P-1))，那么(N^P-N)%P=0可化为：(N*(N^(P-1)-1))%P=0
请注意上式，含义是：N*(N^(P-1)-1)可以被P整除

又因为N*(N^(P-1)-1)必能整除N（这不费话么!）
所以，N*(N^(P-1)-1)是N和P的公倍数，小学知识了^_^

又因为前提是N与P互质，而互质数的最小公倍数为它们的乘积，所以一定存在正整数M使得等式成立：
N*(N^(P-1)-1)=M*N*P
两边约去N，化简之：
N^(P-1)-1=M*P
因为M是整数，显然：
(N^(P-1)-1)%P=0
即：
N^(P-1)%P=1
============================================
积模分解公式

先有一个引理，如果有：X%Z=0，即X能被Z整除，则有：
(X+Y)%Z=Y%Z
这个不用证了吧...

设有X、Y和Z三个正整数，则必有：(X*Y)%Z=((X%Z)*(Y%Z))%Z

想了很长时间才证出来，要分情况讨论才行：

1.当X和Y都比Z大时，必有整数A和B使下面的等式成立：
X=Z*I+A（1）
Y=Z*J+B（2）
不用多说了吧，这是除模运算的性质！
将（1）和（2）代入(X*Y)modZ得：((Z*I+A)(Z*J+B))%Z
乘开，再把前三项的Z提一个出来，变形为：(Z*(Z*I*J+I*A+I*B)+A*B)%Z（3）
因为Z*(Z*I*J+I*A+I*B)是Z的整数倍……晕，又来了。
概据引理，（3）式可化简为：(A*B)%Z
又因为：A=X%Z，B=Y%Z，代入上面的式子，就成了原式了。

2.当X比Z大而Y比Z小时，一样的转化：
X=Z*I+A
代入(X*Y)%Z得：
(Z*I*Y+A*Y)%Z
根据引理，转化得：(A*Y)%Z
因为A=X%Z，又因为Y=Y%Z，代入上式，即得到原式。
同理，当X比Z小而Y比Z大时，原式也成立。

3.当X比Z小，且Y也比Z小时，X=X%Z，Y=Y%Z，所以原式成立。
=====================================================
快速计算乘方的算法

如计算2^13，则传统做法需要进行12次乘法。

/*计算n^p*/
unsigned power(unsigned n,unsigned p)
{
    for(int i=0;i<p;i++) n*=n;
    return n;
}

该死的乘法，是时候优化一下了！把2*2的结果保存起来看看，是不是成了：4*4*4*4*4*4*2
再把4*4的结果保存起来：16*16*16*2
一共5次运算，分别是2*2、4*4和16*16*16*2

这样分析，我们算法因该是只需要计算一半都不到的乘法了。
为了讲清这个算法，再举一个例子2^7：2*2*2*2*2*2*2
两两分开：(2*2)*(2*2)*(2*2)*2
如果用2*2来计算，那么指数就可以除以2了，不过剩了一个，稍后再单独乘上它。
再次两两分开，指数除以2： ((2*2)*(2*2))*(2*2)*2
实际上最后一个括号里的2 * 2是这回又剩下的，那么，稍后再单独乘上它
现在指数已经为1了，可以计算最终结果了：16*4*2=128

优化后的算法如下：
unsigned Power(unsigned n,unsigned p)
{
   unsigned main=n; //用main保存结果
   unsigned odd=1; //odd用来计算“剩下的”乘积
   while (p>1)
   {//一直计算，直到指数小于或等于1
        if((p%2)!=0)
        {// 如果指数p是奇数，则说明计算后会剩一个多余的数，那么在这里把它乘到结果中
            odd*=main; //把“剩下的”乘起来
        }
        main*=main; //主体乘方
        p/=2; //指数除以2
   }
   return main*odd; //最后把主体和“剩下的”乘起来作为结果
}

够完美了吗？不，还不够！看出来了吗？main是没有必要的，并且我们可以有更快的代码来判断奇数。要知道除法或取模运算的效率很低，所以我们可以利用偶数的一个性质来优化代码，那就是偶数的二进制表示法中的最低位一定为0!

完美版：
unsigned Power(unsigned n, unsigned p)
{ // 计算n的p次方
    unsigned odd = 1; //oddk用来计算“剩下的”乘积
    while (p > 1)
    { // 一直计算到指数小于或等于1
       if (( p & 1 )!=0)
      { // 判断p是否奇数，偶数的最低位必为0
             odd *= n; // 若odd为奇数，则把“剩下的”乘起来
      }
      n *= n; // 主体乘方
      p /= 2; // 指数除以2
     }
    return n * odd; // 最后把主体和“剩下的”乘起来作为结果
}
========================================================
蒙格马利”快速幂模算法

后面我们会用到这样一种运算：(X^Y)%Z

问题是当X和Y很大时，只有32位的整型变量如何能够有效的计算出结果？
考虑上面那份最终的优化代码和再上面提到过的积模分解公式，我想你也许会猛拍一下脑门，吸口气说：“哦，我懂了！”。

下面的讲解是给尚没有做出这样动作的同学们准备的。X^Y可以看作Y个X相乘，即然有积模分解公式，那么我们就可以把Y个X相乘再取模的过程分解开来，比 如：(17^25)%29则可分解为：( ( 17 * 17 ) % 29 * ( 17 * 17 ) % 29 * ……
如果用上面的代码将这个过程优化，那么我们就得到了著名的“蒙格马利”快速幂模算法：
unsigned Montgomery(unsigned n, unsigned p, unsigned m)
{ // 快速计算 (n ^ e) % m 的值，与power算法极类似
    unsigned r = n % m; // 这里的r可不能省
    unsigned k = 1;
    while (p > 1)
    {
        if ((p & 1)!=0)
        {
            k = (k * r) % m; // 直接取模
        }
        r = (r * r) % m; // 同上
        p /= 2;
    }
    return (r * k) % m; // 还是同上
}

上面的代码还可以优化。下面是蒙格马利极速版：

unsigned Montgomery(unsigned n,unsigned p,unsigned m)
{ //快速计算(n^e)%m的值
      unsignedk=1;
      n%=m;
     while(p!=1)
     {
         if(0!=(p&1))k=(k*n)%m;
         n=(n*n)%m;
         p>>=1;
    }
    return(n*k)%m;
}
=====================================================
怎么判断一个数是否为素数？

笨蛋的作法：
bool IsPrime(unsigned n)
{
    if (n<2)
    { //小于2的数即不是合数也不是素数
    throw 0;
    }
    for (unsigned i=2;i<n;++i)
    { //和比它小的所有的数相除，如果都除不尽，证明素数
        if (n%i==0)
        {//除尽了，则是合数
            return false;
        }
    }
    return true;
}

一个数去除以比它的一半还要大的数，一定除不尽，所以还用判断吗？？

下面是小学生的做法：
bool IsPrime(unsigned n)
{
    if (n<2)
    {//小于2的数即不是合数也不是素数
        throw 0;
    }
    for(unsigned i=2;i<n/2+1;++i)
    { // 和比它的一半小数相除，如果都除不尽，证明素数
        if ( 0 == n % i )
        { // 除尽了，合数
            return false;
        }
    }
    return true; // 都没除尽，素数
}

一个合数必然可以由两个或多个质数相乘而得到。那么如果一个数不能被比它的一半小的所有的质数整除，那么比它一半小的所有的合数也一样不可能整除它。建立一个素数表是很有用的。

下面是中学生的做法：
bool IsPrime2(unsigned n)
{
    if ( n < 2 )
    { // 小于2的数即不是合数也不是素数
        throw 0;
    }
    static unsigned aPrimeList[] = { // 素数表
        1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,
        43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 113,
        193, 241, 257, 337, 353, 401, 433, 449, 577, 593, 641,
        673, 769, 881, 929, 977, 1009, 1153, 1201, 1217, 1249,
        1297,1361, 1409, 1489, 1553, 1601, 1697, 1777, 1873,
        1889, 2017, 2081, 2113, 2129, 2161, 2273, 2417, 2593,
        2609, 2657, 2689, 2753, 2801, 2833, 2897, 3041, 3089,
        3121, 3137, 3169, 3217, 3313, 3329, 3361, 3457, 3617,
        3697, 3761, 3793, 3889, 4001, 4049, 4129, 4177, 4241,
        4273, 4289, 4337, 4481, 4513, 4561, 4657, 4673, 4721,
        4801, 4817, 4993, 5009, 5153, 5233, 5281, 5297, 5393,
        5441, 5521, 5569, 5857, 5953, 6113, 6257, 6337, 6353,
        6449, 6481, 6529, 6577, 6673, 6689, 6737, 6833, 6961,
        6977, 7057, 7121, 7297, 7393, 7457, 7489, 7537, 7649,
        7681, 7793, 7841, 7873, 7937, 8017, 8081, 8161, 8209,
        8273, 8353, 8369, 8513, 8609, 8641, 8689, 8737, 8753,
        8849, 8929, 9041, 9137, 9281, 9377, 9473, 9521, 9601,
        9649, 9697, 9857
    };
   
    const int nListNum = sizeof(aPrimeList)/sizeof(unsigned);//计算素数表里元素的个数
    for (unsigned i=2;i<nListNum;++i )
    {
        if(n/2+1<aPrimeList[i])
        {
            return true;
        }
        if(0==n%aPrimeList[i])
        {
            return false;
        }
    }
    /*由于素数表中元素个数是有限的，那么对于用素数表判断不到的数，就只有用笨蛋办法了*/
    for (unsigned i=aPrimeList[nListNum-1];i<n/2+1;i++ )
    {
        if (0==n%i)
        { // 除尽了，合数
            return false;
        }
    }
    return true;
}
    还是太糟了，我们现在要做的对于大型素数的判断，那个素数表倒顶个P用！当然，我们可以利用动态的素数表来进行优化，这就是大学生的做法了。但是动态生成素数表的策略又复杂又没有效率，所以我们还是直接跳跃到专家的做法吧：
    根据上面讲到的费马小定理，对于两个互质的素数N和P，必有：N^(P-1)%P=1
    那么我们通过这个性质来判断素数吧，当然，你会担心当P很大的时候乘方会很麻烦。不用担心！我们上面不是有个快速的幂模算法么？好好的利用蒙格马利这位大数学家为我们带来的快乐吧！

算法思路是这样的：
    对于N，从素数表中取出任意的素数对其进行费马测试，如果取了很多个素数，N仍未测试失败，那么则认为N是素数。当然，测试次数越多越准确，但一般来讲 50次就足够了。另外，预先用“小学生”的算法构造一个包括500个素数的数组，先对Q进行整除测试，将会大大提高通过率，方法如下:
bool IsPrime3(unsigned n)
{
    if ( n < 2 )
    { // 小于2的数即不是合数也不是素数
        throw 0;
    }
    static unsigned aPrimeList[] = {
        2, 3, 5, 7, 11, 17, 19, 23, 29, 31, 41,
        43, 47, 53, 59, 67, 71, 73, 79, 83, 89, 97
    };
    const int nListNum = sizeof(aPrimeList) / sizeof(unsigned);
    for (int i=0;i<nListNum;++i)
    { // 按照素数表中的数对当前素数进行判断
        if (1!=Montgomery(aPrimeList[i],n-1,n)) // 蒙格马利算法
        {
            return false;
        }
    }
    return true;
}

    OK，这就专家的作法了。
    等等，什么？好像有点怪，看一下这个数29341，它等于13 * 37 * 61，显然是一个合数，但是竟通过了测试！！哦，抱歉，我忘了在素数表中加入13，37，61这三个数，我其实是故意的，我只是想说明并费马测试并不完全可靠。
    现在我们发现了重要的一点，费马定理是素数的必要条件而非充分条件。这种不是素数，但又能通过费马测试的数字还有不少，数学上把它们称为卡尔麦克数，现在数学家们已经找到所有10 ^ 16以内的卡尔麦克数，最大的一个是9585921133193329。我们必须寻找更为有效的测试方法。数学家们通过对费马小定理的研究，并加以扩展，总结出了多种快速有效的素数测试方法，目前最快的算法是拉宾米勒测试算法，下面介绍拉宾米勒测试。
================================================================
拉宾米勒测试

    拉宾米勒测试是一个不确定的算法，只能从概率意义上判定一个数可能是素数，但并不能确保。算法 流程如下:
    1.选择T个随机数A，并且有A<N成立。
    2.找到R和M，使得N=2*R*M+1成立。
    快速得到R和M的方式：N用二进制数B来表示，令C=B-1。因为N为奇数（素数都是奇数），所以C的最低位为0，从C的最低位的0开始向高位统计，一直到遇到第一个1。这时0的个数即为R，M为B右移R位的值 。
    3.如果A^M%N=1，则通过A对于N的测试，然后进行下一个A的测试
    4.如果A^M%N!=1，那么令i由0迭代至R，进行下面的测试
    5.如果A^((2^i)*M)%N=N-1则通过A对于N的测试，否则进行下一个i的测试
    6.如果i=r，且尚未通过测试，则此A对于N的测试失败，说明N为合数。
    7.进行下一个A对N的测试，直到测试完指定个数的A

    通过验证得知，当T为素数，并且A是平均分布的随机数，那么测试有效率为1 / ( 4 ^ T )。如果T > 8那么测试失误的机率就会小于10^(-5)，这对于一般的应用是足够了。如果需要求的素数极大，或着要求更高的保障度，可以适当调高T的值。下面是代码：

bool RabbinMillerTest( unsigned n )
{
    if (n<2)
    { // 小于2的数即不是合数也不是素数
        throw 0;
    }

    const unsigned nPrimeListSize=sizeof(g_aPrimeList)/sizeof(unsigned);//求素数表元素个数
    for(int i=0;i<nPrimeListSize;++i)
    {// 按照素数表中的数对当前素数进行判断
        if (n/2+1<=g_aPrimeList[i])
        {// 如果已经小于当前素数表的数，则一定是素数
            return true;
        }
        if (0==n%g_aPrimeList[i])
        {// 余数为0则说明一定不是素数
            return false;
        }
    }
    // 找到r和m，使得n = 2^r * m + 1;
    int r = 0, m = n - 1; // ( n - 1 ) 一定是合数
    while ( 0 == ( m & 1 ) )
    {
        m >>= 1; // 右移一位
        r++; // 统计右移的次数
    }
    const unsigned nTestCnt = 8; // 表示进行测试的次数
    for ( unsigned i = 0; i < nTestCnt; ++i )
    { // 利用随机数进行测试，
        int a = g_aPrimeList[ rand() % nPrimeListSize ];
        if ( 1 != Montgomery( a, m, n ) )
        {
            int j = 0;
            int e = 1;
            for ( ; j < r; ++j )
            {
                if ( n - 1 == Montgomery( a, m * e, n ) )
                {
                    break;
                }
                e <<= 1;
            }
            if (j == r)
            {
                return false;
            }
        }
    }
    return true;
}

 

]虽然临近期末考试，做题的欲望却丝毫没有减退。算了，这三天堕落了一下，把一些OJ上含有Fibon字样的题目给割掉了。
下面就要正式的进入复习阶段了，等考试完再回来。
下面是一些题目的解题报告了，后面如果发现新的会补上来（不知道 FH算不算:-)）
1. pku3070 Fibonacci http://acm.pku.edu.cn/JudgeOnline/problem?id=3070 这个题目按照题意做就可以轻松的搞定了.问题关键是给出了我们另外一种求解Fibonacci的方法。利用矩阵乘法！
那为什么它就比递推式和通项公式要好呢,原因是因为Fibonacci的增长速度非常快,通常普通解法没有办法求得
后面的大数,所以题目通常要求取模,而如果用通项,肯定不够精确了,用递推公式又肯定会超时(比如要你求第1000000000个
Fibonacci数模10000000的最终结果),这个时候转换为求矩阵的幂次,可以在logn的时间内搞定.
我比较喜欢下面这个形式
|0 1|    |f(0)|                   |f(1)|
===>
|1 1|    |f(1)|                   |f(0)+f(1)| 后面这个也就是f(2)了
如果令 A为那个0,1矩阵,可以发现
f(n) = (A^n).a12 * f(1) (a12为右上角的元素值)
特别的 ，
A^0 = I (I是单位矩阵)
2. hdu1021 Fibonacci Again http://acm.hdu.edu.cn/showproblem.php?pid=1021 通过一个模3障眼法，其实也可以使用上述的矩阵幂次进行求解.不过看题目的数据量,貌似直接用推倒也可以搞定
3. hdu1568 http://acm.hdu.edu.cn/showproblem.php?pid=1568 输出一个fibonacci的前4个字母
这个要使用fibonacci的通项公式了
((1+√5)/2 )^n - ((1-√5)/2)^ n
--------------------------------
√5
可以发现后面的((1-√5)/2)^ n 当n很大的时候数值非常小，可以忽略，不会对前4位造成影响。
因此转变为了
((1+√5)/2 )^n
------------------
√5
这个式子随着n的增长也是增长很快，但是我们可以发现，题目要求的是前4位，因此我们发现将长度固定在某个范围
之内(即如果数大了，就除一下)，这样的话就可以取到前4位了.
4. hdu1250 Hat's Fibonacci http://acm.hdu.edu.cn/showproblem.php?pid=1250 我用的赤裸裸的高精度算的把结果存起来了,大汗.内存占了18M. Statistic里面的貌似都很少的样子.
后来我用数组去循环递推做的,内存自然就下降下来了.这个题目貌似也可以构造矩阵,只不过要在计算
里面处理好高精度.
5. hit1214 Fibonacci Coding http://acm.hit.edu.cn/ojs/show.php?Proid=1214&Contestid=0 题目大意是去除相邻的1,改变成前面一个.
我是直接模拟的,从高位开始判断,如果相邻位为1且前面一个为0,则变化,直至没有出现该情况为止
要注意的是,如果使用string,且在高位前产生进位的时候,要更新len,如果你没有直接调用size的话.
6. hit1533 Fibonacci Numbers http://acm.hit.edu.cn/ojs/show.php?Proid=1533&Contestid=0 高精度,与hdu1350类似
7. hit1864 Fibonacci http://acm.hit.edu.cn/ojs/show.php?Proid=1864&Contestid=0 这个是求第k个finbonacci数的位数,求位数,用大数就太不明智了.
这里我同样的使用了通项的, 当n比较大的时候,利用变化过的通项
((1+√5)/2 )^n
---------------
√5
那么求位数就用 log10(F) + 1
再利用log10(a*b) = log10(a) + log10(b), log10(a/b) = log10(a) - log10(b)
8. hit2060 Fibonacci Problem Again http://acm.hit.edu.cn/ojs/show.php?Proid=2060&Contestid=0 大意求第a个到第b个之间所有fibonacci的数和。
这个有点小技巧在里面，可以自己推导一下：
F(3) = F(1) + F(2)
F(4) = F(2) + F(3) = 1 * F(1) + 2 * F(2)
F(5) = F(3) + F(4) = 2 * F(1) + 3 * F(2)
F(6) = F(4) + F(5) = 3 * F(1) + 5 * F(2)
F(7) = F(5) + F(6) = 5 * F(1) + 8 * F(2)
F(8) = F(6) + F(7) = 8 * F(1) + 13 * F(2)
S(3) = 2 * F(1) + 2 * F(2)
S(4) = 3 * F(1) + 4 * F(2)
S(5) = 5 * F(1) + 7 * F(2)
S(6) = 8 * F(1) + 12 *F(2)
S(7) = 13 *F(1) + 20 *F(2)
不难发现，
S(n) = F(n + 2) - F(2)
因此题目就转换为了求 F(b + 2) - F(a + 2 - 1)
而又是大数取模，考虑使用矩阵求幂
9. hit2065 Fibonacci Number http://acm.hit.edu.cn/ojs/show.php?Proid=2065&Contestid=0 赤裸裸的Fibonacci，最基本的那种
10. hit2204 Fibonacci Numbers http://acm.hit.edu.cn/ojs/show.php?Proid=2204&Contestid=0 这个题目相当于求第k个Fibonacci的前4位和后4位，那么前面的题目
综合起来，就是前4位用通项，后4位用矩阵求模。
就不难了。要注意的是，前面求前4位的时候我用的O(n)的做法，过了。
这边如果仍然采用O(n)的会TLE，必须采用O(logn),求幂次的方法才能过.
11. hit2255 Not Fibonacci http://acm.hit.edu.cn/ojs/show.php?Proid=2255&Contestid=0 这个题目是个BT题,交了我n次,后来发现是s=0的时候负数情况没有处理,汗.
粗看一下，这个题目和hit2060一样，不过这里好像都泛化了。这里面的规律
可没有上面那题那么好找了。wy他们貌似直接推导S的公式，都是采用的什么
等比数列弄的,我和felica讨论了一下，发现了一个好的解法。通过构造分块
矩阵，可以很好的解决这个问题。
首先构造
|0 1| |f(0)|
|q p| |f(1)|
这个可以利用矩阵幂次求得任何一个F(n),可惜这里是求和,
那么令A=前面那个p,q的矩阵
则S(n) = f(0) + f(1) + ... + f(n) = f(0) + (A + A^2 + ... + A^n)*F所得矩阵的右上角的值
那么如何求得A + A^2 + ... + A^n呢
可以继续构造如下的分块矩阵，其中I是单位矩阵
|A I|
|0 I|
令R等于上面的矩阵，则
R^2 = |A^2   A*I + I|
|0     I      |
R^3 = |A^3   A^2 * I + A * I + I|
| 0     I                 |
可以发现右上角即为I + A + A^2 + ... + A^n，多个I后面给减掉就可以了
这样我们同样可以再次利用矩阵幂次求得R^n
则S(n) = ((R^(n+1).m12 - I) * F).a12 + f(0)
可以很完美的解决这个问题.
12. hnu10072 Fibonacci Number http://acm.hnu.cn:8080/online/?action=problem&type=show&id=10072&courseid=0 与hit2065一样
13. tju1208 Fibonacci Numbers http://acm.tju.edu.cn/toj/showp1208.html 与hit1533一样
14. zju2060 Fibonacci Again http://acm.zju.edu.cn/show_problem.php?pid=2060 与hdu1021一样
15. zju1828 Fibonacci Numbers http://acm.zju.edu.cn/show_problem.php?pid=1828 与hit1533一样
16. zju2672 Fibonacci Subsequence http://acm.zju.edu.cn/show_problem.php?pid=2672 又是一道bt题,时间卡的非常之紧.
这个题目我使用的dp+hash.
记d(i,j)为以第i个数作为结尾,前一个数是第j个的fibonacci串的最大长度
则 d(i,j) = max{d(j, i - j)} + 1,
注意查找data[i] - data[j]的时候利用hash查找到的应该是j小的所有当中最大下标的那个。
因此需要在每次循环结束的时候更新已存在的数的hash。

17. hzu1060 Fibonacci数列 http://acm.fzu.edu.cn/problem.php?pid=1060 与hit1533相同
18. fjnu1158 Fibonacci数列 http://acm.fjnu.edu.cn/show?problem_id=1158 与hit1533相同
19. zjut Fibonacci数 http://acm.zjut.edu.cn/ShowProblem.aspx?ShowID=1029 与hit2065一样,题目描述稍有不一样

总结一下，主要能够考察的是
1. 求解Fibonacci的某一项（这个范围一般在45之内)
2. 求解Fibonacci的某一项模K(这个一般是大数),通用解决方法是构造矩阵求幂次
3. 求解Fibonacci的前多少位 (这个一般是大数), 通用解法是使用通项公式
4. 求解Fibonacci的后多少位，这个和取模类似
5. 求解Fibonacci的前n项和，利用推导式, S(n) = F(n + 2) - F(2)即可
6. 更多的是基于Fibonacci的综合题,包括DP,构造,等等.
=======================================================
另外,gardon补充了一题:
就普通的f(0)=0,f(1)=1,f(n)=f(n-1)+f(n-2)
给数列an=k*n+b
保证k>0,b>=0
求1到n的f(an)之和
输入k、b、m，求出那个和模m的值
b, m, n数的范围都是10^9以内
k <=20
这个题目仍然可以使用hit那个BT题的做法。先构造矩阵
应该是
0 1
1 1的幂次了(主要是使得F'(n + 1) = F'(n)*矩阵),然后构造求和矩阵就可以了.
矩阵的方法。
思想来源上上次PKU的acm warmup http://acm.pku.edu.cn/JudgeOnline/problem?id=3070 {{F(n+1),F(n)},
{F(n),F(n-1)}}
=
{{1,1}      ^n
{1,0}}
我当时就想，用加法算都只能是O(n)的，换成n个矩阵的连乘不是更废劲吗？其实不然，乘法方法有比加法更好的性质，用加法只能利用前两个的值，而乘法却不同，因为乘法有结合律，可以大幅度下降算法的耗时。
令A(n)表示一个矩阵序列
A(n)=
{{F(n),F(n-1)},
{F(n-1),F(n-2)}
那么A(2)={{1,1},{1,0}}，由那个表达式得到：
A(n)=A(2)^(n-1)，A(2)是己知的2*2矩阵，现在的问题就是如何求
A(2)^n
因为方阵的乘法有结合律，所以
A(2)^n=A(2)^(n/2)*A(2)^(n/2)，不妨设n是偶数
所以求A(n)就可以化成求A(n/2)并作一次乘法，所以递归方程是：
T(n)=T(n/2)+O(1)
它的解是O(lbn)，不可思议吧。
最后我们要观察这样的矩阵方法可否推广到一般情形，我们暂且称这样的递推序列为‘c阶线性递推方程’吧，它的意思就是说F(n)由与n相距c远的前面的一些F(m)组成的线性组合，比如Fibonacci数列的最大距离就是2，所以这样的递推式需要确定一个序列至少要知道最初的c个己知初值，我们的目的就是要求这样的递归式的一般算法：
F(n)={F(n-1),F(n-2),...F(n-c)}*K
{F(1),F(2),,,F(c)}^T=C
其中K={x1,x2,...xc}^T，C={c1,c2,...cc}^T都是非0的常向量。
我们构造A(n)表示一个c*c的矩阵序列，且形式为
A(n)=
{{F(n),F(n-1),F(n-2),...F(n-c+1)},
{F(n-1),F(n-2),F(n-3),...F(n-c)},
......
{F(n-c+1),F(n-c),F(n-c-c),...F(n-2c+2)}}
为了使它有意义，n-2c+2>=1，所以n>=2c-1，所以A(n)的初始矩阵是A(2c-1)，可以按递归式直接构造它，然后我们考虑A(n)的生成矩阵会是什么样子，设为G，于是应满足：
A(n)=A(n-1)G
观察A(n-1)的第一行，它刚好就是{F(n-1),F(n-2),...F(n-c)}，由递推式，知G的第一列应为K，而第二列设为x，则要满足{F(n-1),F(n-2),...F(n-c)}x=F(n-1)，很简单，x={1,0,0,...}，同样第三列是{0,1,0,0,...}，所以右边应该是一个(c-1)*(c-1)的单位矩阵I再在最下面补上一行0，于是G可以表示为：
G={K,{I,0}^T}
于是A(n)的通项就是
A(n)=A(2c-1)*G^(n-2c+1)
比如Fibonacci数列的K={1,1}^T，所以
G={K,{I,0}^T}=
{{1,1},
{1,0}}
再举一个例子，也是以前做过的一个题，母牛生小牛的题，在一个农场上有很多母牛，第1年有一只母牛，而且每过4年，母牛会开始生小母牛，小母牛过4年后也会开始生，此后就每年生一只小母年，问第n年时有多少只母牛？
简单例出前几个数字：1，1，1，2，3，4，6，。。。会发现一个规律，就是第n年的母牛数F(n)恰等于前一年的母牛数F(n-1)，以及前三年的母牛数F(n-3)之和，其实直接得到这个递推式也不难，因为F(n-3)到F(n)刚好4年（第4年），所以F(n-3)都是一些可以生仔的牛数，它们会添加F(n-3)个小牛，同时F(n-1)是所有的牛妈妈，那当然有F(n)=F(n-1)+F(n-3)。
它的K={1,0,1},c=3，生成矩阵G=
{{1,1,0},
{0,0,1},
{1,0,0}}
所以有A(n)=A(3)*G^(n-2)
A(3)=
{{1,1,1},
{1,1,0},
{1,0,0}}
从O(2^n)到O(n)，再到O(lbn)，效率提高得也太猛烈了，不过一点要提出来，不管什么样的算法，Fiboacci数列本身数值的递增是指数型的，所以要想算很高阶的结果，也近乎不太可能，所以经常会在一个模数范围内考虑。

 

http://hi.baidu.com/xiehuixb/blog/item/b6d84a878a7c513467096ee3.html

 

本文来自CSDN博客，转载请标明出处：http://blog.csdn.net/q3498233/archive/2010/02/26/5327805.aspx