斐波那契序列 集锦


[定理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

优化后的算法如下:

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

完美版:
========================================================
蒙格马利”快速幂模算法

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

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

下面的讲解是给尚没有做出这样动作的同学们准备的。X^Y可以看作Y个X相乘,即然有积模分解公式,那么我们就可以把Y个X相乘再取模的过程分解开来,比 如:(17^25)%29则可分解为:( ( 17 * 17 ) % 29 * ( 17 * 17 ) % 29 * ……
如果用上面的代码将这个过程优化,那么我们就得到了著名的“蒙格马利”快速幂模算法:

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

=====================================================
怎么判断一个数是否为素数?

笨蛋的作法:

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

下面是小学生的做法:

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

下面是中学生的做法:
 

    还是太糟了,我们现在要做的对于大型素数的判断,那个素数表倒顶个P用!当然,我们可以利用动态的素数表来进行优化,这就是大学生的做法了。但是动态生成素数表的策略又复杂又没有效率,所以我们还是直接跳跃到专家的做法吧:
    根据上面讲到的费马小定理,对于两个互质的素数N和P,必有:N^(P-1)%P=1
    那么我们通过这个性质来判断素数吧,当然,你会担心当P很大的时候乘方会很麻烦。不用担心!我们上面不是有个快速的幂模算法么?好好的利用蒙格马利这位大数学家为我们带来的快乐吧!

算法思路是这样的:
    对于N,从素数表中取出任意的素数对其进行费马测试,如果取了很多个素数,N仍未测试失败,那么则认为N是素数。当然,测试次数越多越准确,但一般来讲 50次就足够了。另外,预先用“小学生”的算法构造一个包括500个素数的数组,先对Q进行整除测试,将会大大提高通过率,方法如下:

    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的值。下面是代码:

 

 

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值