又是黑色星期五!为什么总是黑色星期五呢?

    你是否还记得,今年的4月13日是黑色星期五。短短三个月后,黑色星期五再次现身!为什么这一天老是出现呢?
    恐怕心理原因是最好的解释。人们对黑色星期五的出现记忆更深刻,给人一种黑色星期五常常出现的错觉。有趣的是,仔细算一算,你会发现13日是星期五的次数真的要多一些。
    很多人以为现在实行的历法是4年一循环,这是不对的。现在实行的历法以400年为一个循环。大家很容易忽略整百年的问题。一个很有意思的智力题就是问一个人是否可能连续5年不过生日。有个笑话说文科MM感叹她等了她男友4年,整整1460天,学理科的好友脱口而出“难道你男友是1900年的”。我们可以算一下在这400年中共有多少天:365*303+366*97=146097。这个数正好能被7整除。换句话说,现在与400年后的星期数不变,日历完全相同。在这400年里一共有4800个月,利用Zeller公式(见这里的最后一小节)可以编程统计出13日是星期几的次数最多,这对于OIer们再熟悉不过了,因为USACO有一道题就是干这种无聊的事情。下面就是程序运行后的结果:

  • Sunday     687
  • Monday     685
  • Tuesday    685
  • Wednesday  687
  • Thursday   684
  • Friday     688
  • Saturday   684


    可以看到,事实上13日是星期五的概率确实是最高的。
    另外,注意到了么,利用“400年一周期”这个结论我们可以对USACO的那个题进行扩展,出一个Friday the Thirteenth数据加强版。
    做人要厚道,转贴不注明出处者将受到黑色星期五的诅咒。

经典谬论:用复数来证明1=2

    从来没有见到过一个纯数学的东西在digg上这么受欢迎,我也转上来。

    大家肯定都见过用除以零、平方根、无穷级数等“技巧”去“证明”类似的结论,但我第一次看到这个用虚数来玩的戏法。大家看看哪里错了:

Step 1: -1/1 = 1/-1
Step 2: 两边同时开方:sqrt( -1/1 ) = sqrt( 1/-1 )
Step 3: 化简得:sqrt(-1)/sqrt(1) = sqrt(1)/sqrt(-1)
Step 4: 也就是说,i/1 = 1/i
Step 5: 那么,i / 2 = 1 / (2i)
Step 6: 两边同时加一个数:i/2 + 3/(2i) = 1/(2i) + 3/(2i)
Step 7: 同时乘以一个数:i (i/2 + 3/(2i) ) = i ( 1/(2i) + 3/(2i) )
Step 8: 展开:(i^2)/2 + (3i)/(2i) = i/(2i) + (3i)/(2i)
Step 9: 于是有:(-1)/2 + 3/2 = 1/2 + 3/2
Step 10: 这说明1=2

    如果你也被搞晕了,去这里看看吧。这对于我这种课本上根本没有讲虚数,只是道听途说知道一些东西的文科生来说尤其具有迷惑性。

神奇的分形艺术(一):无限长的曲线可能围住一块有限的面积

    很多东西都是吹神了的,其中麦田圈之谜相当引人注目。上个世纪里人们时不时能听见某个农民早晨醒了到麦田地一看立马吓得屁滚尿流的故事。上面这幅图就是97年在英国Silbury山上发现的麦田圈,看上去大致上是一个雪花形状。你或许会觉得这个图形很好看。看了下面的文字后,你会发现这个图形远远不是“好看”可以概括的,它的背后还有很多东西。

    在说明什么是分形艺术前,我们先按照下面的方法构造一个图形。看下图,首先画一个线段,然后把它平分成三段,去掉中间那一段并用两条等长的线段代替。这样,原来的一条线段就变成了四条小的线段。用相同的方法把每一条小的线段的中间三分之一替换为等边三角形的两边,得到了16条更小的线段。然后继续对16条线段进行相同的操作,并无限地迭代下去。下图是这个图形前五次迭代的过程,可以看到这样的分辨率下已经不能显示出第五次迭代后图形的所有细节了。这样的图形可以用Logo语言很轻松地画出来。

    你可能注意到一个有趣的事实:整个线条的长度每一次都变成了原来的4/3。如果最初的线段长为一个单位,那么第一次操作后总长度变成了4/3,第二次操作后总长增加到16/9,第n次操作后长度为(4/3)^n。毫无疑问,操作无限进行下去,这条曲线将达到无限长。难以置信的是这条无限长的曲线却“始终只有那么大”。

        
    当把三条这样的曲线头尾相接组成一个封闭图形时,有趣的事情发生了。这个雪花一样的图形有着无限长的边界,但是它的总面积却是有限的。换句话说,无限长的曲线围住了一块有限的面积。有人可能会问为什么面积是有限的。虽然从上面的图上看结论很显然,但这里我们还是要给出一个简单的证明。三条曲线中每一条的第n次迭代前有4^(n-1)个长为(1/3)^(n-1)的线段,迭代后多出的面积为4^(n-1)个边长为(1/3)^n的等边三角形。把4^(n-1)扩大到4^n,再把所有边长为(1/3)^n的等边三角形扩大为同样边长的正方形,总面积仍是有限的,因为无穷级数Σ4^n/9^n显然收敛。这个神奇的雪花图形叫做Koch雪花,其中那条无限长的曲线就叫做Koch曲线。他是由瑞典数学家Helge von Koch最先提出来的。本文最开头提到的麦田圈图形显然是想描绘Koch雪花。

    分形这一课题提出的时间比较晚。Koch曲线于1904年提出,是最早提出的分形图形之一。我们仔细观察一下这条特别的曲线。它有一个很强的特点:你可以把它分成若干部分,每一个部分都和原来一样(只是大小不同)。这样的图形叫做“自相似”图形(self-similar),它是分形图形(fractal)最主要的特征。自相似往往都和递归、无穷之类的东西联系在一起。比如,自相似图形往往是用递归法构造出来的,可以无限地分解下去。一条Koch曲线中包含有无数大小不同的Koch曲线。你可以对这条曲线的尖端部分不断放大,但你所看到的始终和最开始一样。它的复杂性不随尺度减小而消失。另外值得一提的是,这条曲线是一条连续的,但处处不光滑(不可微)的曲线。曲线上的任何一个点都是尖点。

    分形图形有一种特殊的计算维度的方法。我们可以看到,在有限空间内就可以达到无限长的分形曲线似乎已经超越了一维的境界,但说它是二维图形又还不够。Hausdorff维度就是专门用来对付这种分形图形的。简单地说,Hausdorff维度描述分形图形中整个图形的大小与一维大小的关系。比如,正方形是一个分形图形,因为它可以分成四个一模一样的小正方形,每一个小正方形的边长都是原来的1/2。当然,你也可以把正方形分成9个边长为1/3的小正方形。事实上,一个正方形可以分割为a^2个边长为1/a的小正方形。那个指数2就是正方形的维度。矩形、三角形都是一样,给你a^2个同样的形状才能拼成一个边长为a倍的相似形,因此它们都是二维的。我们把这里的“边长”理解为一维上的长度,那个1/a则是两个相似形的相似比。如果一个自相似形包含自身N份,每一份的一维大小都是原来的1/s,则这个相似形的Hausdorff维度为log(N)/log(s)。一个立方体可以分成8份,每一份的一维长度都是原来的一半,因此立方体的维度为log(8)/log(2)=3。同样地,一个Koch曲线包含四个小Koch曲线,大小两个Koch曲线的相似比为1/3,因此Koch曲线的Hausdorff维度为log(4)/log(3)。它约等于1.26,是一个介于1和2之间的实数。

    我们常说分形图形是一门艺术。把不同大小的Koch雪花拼接起来可以得到很多美丽的图形。如果有MM看了前面的文字一句也不懂,下面这些图片或许会让你眼前一亮。

放图前留下一句话:
Matrix67原创
转贴请注明出处

同余运算及其基本性质

    100除以7的余数是2,意思就是说把100个东西七个七个分成一组的话最后还剩2个。余数有一个严格的定义:假如被除数是a,除数是b(假设它们均为正整数),那么我们总能够找到一个小于b的自然数r和一个整数m,使得a=bm+r。这个r就是a除以b的余数,m被称作商。我们经常用mod来表示取余,a除以b余r就写成a mod b = r。
    如果两个数a和b之差能被m整除,那么我们就说a和b对模数m同余(关于m同余)。比如,100-60除以8正好除尽,我们就说100和60对于模数8同余。它的另一层含义就是说,100和60除以8的余数相同。a和b对m同余,我们记作a≡b(mod m)。比如,刚才的例子可以写成100≡60(mod 8)。你会发现这种记号到处都在用,比如和数论相关的书中就经常把a mod 3 = 1写作a≡1(mod 3)。
    之所以把同余当作一种运算,是因为同余满足运算的诸多性质。比如,同余满足等价关系。具体地说,它满足自反性(一个数永远和自己同余)、对称性(a和b同余,b和a也就同余)和传递性(a和b同余,b和c同余可以推出a和c同余)。这三个性质都是显然的。
    同余运算里还有稍微复杂一些的性质。比如,同余运算和整数加减法一样满足“等量加等量,其和不变”。小学我们就知道,等式两边可以同时加上一个相等的数。例如,a=b可以推出a+100=b+100。这样的性质在同余运算中也有:对于同一个模数m,如果a和b同余,x和y同余,那么a+x和b+y也同余。在我看来,这个结论几乎是显然的。当然,我们也可以严格证明这个定理。这个定理对减法同样有效。

    性质:如果a≡b(mod m),x≡y(mod m),则a+x≡b+y(mod m)。
    证明:条件告诉我们,可以找到p和q使得a-mp = b-mq,也存在r和s使得x-mr = y-ms。于是a-mp + x-mr = b-mq + y-ms,即a+x-m(p+r) = b+y-m(q+s),这就告诉我们a+x和b+y除以m的余数相同。

    容易想到,两个同余式对应相乘,同余式两边仍然相等:
    如果a≡b(mod m),x≡y(mod m),则ax≡by(mod m)。
    证明:条件告诉我们,a-mp = b-mq,x-mr = y-ms。于是(a-mp)(x-mr) = (b-mq)(y-ms),等式两边分别展开后必然是ax-m(…) = by-m(…)的形式,这就说明ax≡by(mod m)。

    现在你知道为什么有的题要叫你“输出答案mod xxxxx的结果”了吧,那是为了避免高精度运算,因为这里的结论告诉我们在运算过程中边算边mod和算完后再mod的结果一样。假如a是一个很大的数,令b=a mod m,那么(a * 100) mod m和(b * 100) mod m的结果是完全一样的,这相当于是在a≡b (mod m)的两边同时乘以100。这些结论其实都很显然,因为同余运算只关心余数(不关心“整的部分”),完全可以每一次运算后都只保留余数。因此,整个运算过程中参与运算的数都不超过m,避免了高精度的出现。

    在证明Fermat小定理时,我们用到了这样一个定理:
    如果ac≡bc(mod m),且c和m互质,则a≡b(mod m) (就是说同余式两边可以同时除以一个和模数互质的数)。
    证明:条件告诉我们,ac-mp = bc-mq,移项可得ac-bc = mp-mq,也就是说(a-b)c = m(p-q)。这表明,(a-b)c里需要含有因子m,但c和m互质,因此只有可能是a-b被m整除,也即a≡b(mod m)。

    可能以后还要用到更多的定理,到时候在这里更新。

Matrix67原创
转贴请注明出处

数论部分第一节:素数与素性测试

    一个数是素数(也叫质数),当且仅当它的约数只有两个——1和它本身。规定这两个约数不能相同,因此1不是素数。对素数的研究属于数论范畴,你可以看到许多数学家没事就想出一些符合某种性质的素数并称它为某某某素数。整个数论几乎就围绕着整除和素数之类的词转过去转过来。对于写代码的人来说,素数比想像中的更重要,Google一下BigPrime或者big_prime你总会发现大堆大堆用到了素数常量的程序代码。平时没事时可以记一些素数下来以备急用。我会选一些好记的素数,比如4567, 124567, 3214567, 23456789, 55566677, 1234567894987654321, 11111111111111111111111 (23个1)。我的手机号前10位是个素数。我的网站域名的ASCII码连起来(77 97 116 114 105 120 54 55 46 99 111 109)也是个素数。还有,我的某个MM的八位生日也是一个素数。每次写Hash函数之类的东西需要一个BigPrime常量时我就取她的生日,希望她能给我带来好运。偶尔我叫她素MM,没人知道是啥意思,她自己也不知道。
    素数有很多神奇的性质。我写5个在下面供大家欣赏。

1. 素数的个数无限多(不存在最大的素数)
  证明:反证法,假设存在最大的素数P,那么我们可以构造一个新的数2 * 3 * 5 * 7 * … * P + 1(所有的素数乘起来加1)。显然这个数不能被任一素数整除(所有素数除它都余1),这说明我们找到了一个更大的素数。

2. 存在任意长的一段连续数,其中的所有数都是合数(相邻素数之间的间隔任意大)
  证明:当0<a<=n时,n!+a能被a整除。长度为n-1的数列n!+2, n!+3, n!+4, …, n!+n中,所有的数都是合数。这个结论对所有大于1的整数n都成立,而n可以取到任意大。

3. 所有大于2的素数都可以唯一地表示成两个平方数之差。
  证明:大于2的素数都是奇数。假设这个数是2n+1。由于(n+1)^2=n^2+2n+1,(n+1)^2和n^2就是我们要找的两个平方数。下面证明这个方案是唯一的。如果素数p能表示成a^2-b^2,则p=a^2-b^2=(a+b)(a-b)。由于p是素数,那么只可能a+b=p且a-b=1,这给出了a和b的唯一解。

4. 当n为大于2的整数时,2^n+1和2^n-1两个数中,如果其中一个数是素数,那么另一个数一定是合数。
  证明:2^n不能被3整除。如果它被3除余1,那么2^n-1就能被3整除;如果被3除余2,那么2^n+1就能被3整除。总之,2^n+1和2^n-1中至少有一个是合数。

5. 如果p是素数,a是小于p的正整数,那么a^(p-1) mod p = 1。
  这个证明就有点麻烦了。
    首先我们证明这样一个结论:如果p是一个素数的话,那么对任意一个小于p的正整数a,a, 2a, 3a, …, (p-1)a除以p的余数正好是一个1到p-1的排列。例如,5是素数,3, 6, 9, 12除以5的余数分别为3, 1, 4, 2,正好就是1到4这四个数。
    反证法,假如结论不成立的话,那么就是说有两个小于p的正整数m和n使得na和ma除以p的余数相同。不妨假设n>m,则p可以整除a(n-m)。但p是素数,那么a和n-m中至少有一个含有因子p。这显然是不可能的,因为a和n-m都比p小。
    用同余式表述,我们证明了:
(p-1)! ≡ a * 2a * 3a * … * (p-1)a (mod p)
    也即:
(p-1)! ≡ (p-1)! * a^(p-1) (mod p)
    两边同时除以(p-1)!,就得到了我们的最终结论:
1 ≡ a^(p-1) (mod p)

    可惜最后这个定理最初不是我证明的。这是大数学家Fermat证明的,叫做Fermat小定理(Fermat's Little Theorem)。Euler对这个定理进行了推广,叫做Euler定理。Euler一生的定理太多了,为了和其它的“Euler定理”区别开来,有些地方叫做Fermat小定理的Euler推广。Euler定理中需要用一个函数f(m),它表示小于m的正整数中有多少个数和m互素(两个数只有公约数1称为互素)。为了方便,我们通常用记号φ(m)来表示这个函数(称作Euler函数)。Euler指出,如果a和m互素,那么a^φ(m) ≡ 1 (mod m)。可以看到,当m为素数时,φ(m)就等于m-1(所有小于m的正整数都与m互素),因此它是Fermat小定理的推广。定理的证明和Fermat小定理几乎相同,只是要考虑的式子变成了所有与m互素的数的乘积:m_1 * m_2 … m_φ(m) ≡ (a * m_1)(a * m_2) … (a * m_φ(m)) (mod m)。我为什么要顺便说一下Euler定理呢?因为下面一句话可以增加我网站的PV:这个定理出现在了The Hundred Greatest Theorems里。

    谈到Fermat小定理,数学历史上有很多误解。很长一段时间里,人们都认为Fermat小定理的逆命题是正确的,并且有人亲自验证了a=2, p<300的所有情况。国外甚至流传着一种说法,认为中国在孔子时代就证明了这样的定理:如果n整除2^(n-1)-1,则n就是素数。后来某个英国学者进行考证后才发现那是因为他们翻译中国古文时出了错。1819年有人发现了Fermat小定理逆命题的第一个反例:虽然2的340次方除以341余1,但341=11*31。后来,人们又发现了561, 645, 1105等数都表明a=2时Fermat小定理的逆命题不成立。虽然这样的数不多,但不能忽视它们的存在。于是,人们把所有能整除2^(n-1)-1的合数n叫做伪素数(pseudoprime),意思就是告诉人们这个素数是假的。
    不满足2^(n-1) mod n = 1的n一定不是素数;如果满足的话则多半是素数。这样,一个比试除法效率更高的素性判断方法出现了:制作一张伪素数表,记录某个范围内的所有伪素数,那么所有满足2^(n-1) mod n = 1且不在伪素数表中的n就是素数。之所以这种方法更快,是因为我们可以使用二分法快速计算2^(n-1) mod n 的值,这在计算机的帮助下变得非常容易;在计算机中也可以用二分查找有序数列、Hash表开散列、构建Trie树等方法使得查找伪素数表效率更高。
    有人自然会关心这样一个问题:伪素数的个数到底有多少?换句话说,如果我只计算2^(n-1) mod n的值,事先不准备伪素数表,那么素性判断出错的概率有多少?研究这个问题是很有价值的,毕竟我们是OIer,不可能背一个长度上千的常量数组带上考场。统计表明,在前10亿个自然数中共有50847534个素数,而满足2^(n-1) mod n = 1的合数n有5597个。这样算下来,算法出错的可能性约为0.00011。这个概率太高了,如果想免去建立伪素数表的工作,我们需要改进素性判断的算法。

    最简单的想法就是,我们刚才只考虑了a=2的情况。对于式子a^(n-1) mod n,取不同的a可能导致不同的结果。一个合数可能在a=2时通过了测试,但a=3时的计算结果却排除了素数的可能。于是,人们扩展了伪素数的定义,称满足a^(n-1) mod n = 1的合数n叫做以a为底的伪素数(pseudoprime to base a)。前10亿个自然数中同时以2和3为底的伪素数只有1272个,这个数目不到刚才的1/4。这告诉我们如果同时验证a=2和a=3两种情况,算法出错的概率降到了0.000025。容易想到,选择用来测试的a越多,算法越准确。通常我们的做法是,随机选择若干个小于待测数的正整数作为底数a进行若干次测试,只要有一次没有通过测试就立即把
这个数扔回合数的世界。这就是Fermat素性测试。
    人们自然会想,如果考虑了所有小于n的底数a,出错的概率是否就可以降到0呢?没想到的是,居然就有这样的合数,它可以通过所有a的测试(这个说法不准确,详见我在地核楼层的回复)。Carmichael第一个发现这样极端的伪素数,他把它们称作Carmichael数。你一定会以为这样的数一定很大。错。第一个Carmichael数小得惊人,仅仅是一个三位数,561。前10亿个自然数中Carmichael数也有600个之多。Carmichael数的存在说明,我们还需要继续加强素性判断的算法。

    Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步,建立了传说中的Miller-Rabin素性测试算法。新的测试基于下面的定理:如果p是素数,x是小于p的正整数,且x^2 mod p = 1,那么要么x=1,要么x=p-1。这是显然的,因为x^2 mod p = 1相当于p能整除x^2-1,也即p能整除(x+1)(x-1)。由于p是素数,那么只可能是x-1能被p整除(此时x=1)或x+1能被p整除(此时x=p-1)。
    我们下面来演示一下上面的定理如何应用在Fermat素性测试上。前面说过341可以通过以2为底的Fermat测试,因为2^340 mod 341=1。如果341真是素数的话,那么2^170 mod 341只可能是1或340;当算得2^170 mod 341确实等于1时,我们可以继续查看2^85除以341的结果。我们发现,2^85 mod 341=32,这一结果摘掉了341头上的素数皇冠,面具后面真实的嘴脸显现了出来,想假扮素数和我的素MM交往的企图暴露了出来。
    这就是Miller-Rabin素性测试的方法。不断地提取指数n-1中的因子2,把n-1表示成d*2^r(其中d是一个奇数)。那么我们需要计算的东西就变成了a的d*2^r次方除以n的余数。于是,a^(d * 2^(r-1))要么等于1,要么等于n-1。如果a^(d * 2^(r-1))等于1,定理继续适用于a^(d * 2^(r-2)),这样不断开方开下去,直到对于某个i满足a^(d * 2^i) mod n = n-1或者最后指数中的2用完了得到的a^d mod n=1或n-1。这样,Fermat小定理加强为如下形式:
    尽可能提取因子2,把n-1表示成d*2^r,如果n是一个素数,那么或者a^d mod n=1,或者存在某个i使得a^(d*2^i) mod n=n-1 ( 0<=i<r ) (注意i可以等于0,这就把a^d mod n=n-1的情况统一到后面去了)
    Miller-Rabin素性测试同样是不确定算法,我们把可以通过以a为底的Miller-Rabin测试的合数称作以a为底的强伪素数(strong pseudoprime)。第一个以2为底的强伪素数为2047。第一个以2和3为底的强伪素数则大到1 373 653。
    Miller-Rabin算法的代码也非常简单:计算d和r的值(可以用位运算加速),然后二分计算a^d mod n的值,最后把它平方r次。程序的代码比想像中的更简单,我写一份放在下边。虽然我已经转C了,但我相信还有很多人看不懂C语言。我再写一次Pascal吧。函数IsPrime返回对于特定的底数a,n是否是能通过测试。如果函数返回False,那说明n不是素数;如果函数返回True,那么n极有可能是素数。注意这个代码的数据范围限制在longint,你很可能需要把它们改成int64或高精度计算。
function pow( a, d, n:longint ):longint;
begin
   if d=0 then exit(1)
   else if d=1 then exit(a)
   else if d and 1=0 then exit( pow( a*a mod n, d div 2, n) mod n)
   else exit( (pow( a*a mod n, d div 2, n) * a) mod n);
end;

function IsPrime( a,n:longint ):boolean;
var
   d,t:longint;
begin
   if n=2 then exit(true);
   if (n=1) or (n and 1=0) then exit(false);
   d:=n-1;
   while d and 1=0 do d:=d shr 1;
   t:=pow( a, d, n );
   while ( d<>n-1 ) and ( t<>1 ) and ( t<>n-1 ) do
   begin
      t:=(t * t)mod n;
      d:=d shl 1;
   end;
   exit( (t=n-1) or (d and 1=1) );
end;

    对于大数的素性判断,目前Miller-Rabin算法应用最广泛。一般底数仍然是随机选取,但当待测数不太大时,选择测试底数就有一些技巧了。比如,如果被测数小于4 759 123 141,那么只需要测试三个底数2, 7和61就足够了。当然,你测试的越多,正确的范围肯定也越大。如果你每次都用前7个素数(2, 3, 5, 7, 11, 13和17)进行测试,所有不超过341 550 071 728 320的数都是正确的。如果选用2, 3, 7, 61和24251作为底数,那么10^16内唯一的强伪素数为46 856 248 255 981。这样的一些结论使得Miller-Rabin算法在OI中非常实用。通常认为,Miller-Rabin素性测试的正确率可以令人接受,随机选取k个底数进行测试算法的失误率大概为4^(-k)。

    Miller-Rabin算法是一个RP算法。RP是时间复杂度的一种,主要针对判定性问题。一个算法是RP算法表明它可以在多项式的时间里完成,对于答案为否定的情形能够准确做出判断,但同时它也有可能把对的判成错的(错误概率不能超过1/2)。RP算法是基于随机化的,因此多次运行该算法可以降低错误率。还有其它的素性测试算法也是概率型的,比如Solovay-Strassen算法。另外一些素性测试算法则需要预先知道一些辅助信息(比如n-1的质因子),或者需要待测数满足一些条件(比如待测数必须是2^n-1的形式)。前几年AKS算法轰动世界,它是第一个多项式的、确定的、无需其它条件的素性判断算法。当时一篇论文发表出来,题目就叫PRIMES is in P,然后整个世界都疯了,我们班有几个MM那天还来了初潮。算法主要基于下面的事实:n是一个素数当且仅当(x-a)^n≡(x^n-a) (mod n)。注意这个x是多项式中的未知数,等式两边各是一个多项式。举个例子来说,当a=1时命题等价于如下结论:当n是素数时,杨辉三角的第n+1行除两头的1以外其它的数都能被n整除。

Matrix67原创
转贴请注明出处