Jul 13

    Proofs from THE BOOK的第六章相当精彩,这一章循序渐进地介绍了多个无理性证明。先证明e是无理数,证明方法和高数课本上的基本相同;试图用类似的办法证明e^2也是无理数时,这一章的内容开始牛B了起来,一些巧妙的变换就让原来的办法继续适用于e^2的证明;加上一些更有趣的技巧,我们还能继续证明e^4也是无理数;当证明对除0外的所有有理数r,e^r都是无理数时,全章达到了高潮。
    这一章还提到了pi^2是无理数的证明方法。这个证明建立在Ivan Niven于1947年提出的“pi是无理数”的经典证明的基础上:仅仅是在原证明过程中加了一些微妙的变化就得到了pi^2也是无理数的结论。注意到,“pi^2是无理数”是一个比“pi是无理数”更强的结论。由于有理数的平方还是有理数,因此证到了pi^2是无理数也就说明了pi必然是无理数;但反过来却不行,因为无理数的平方不一定也是无理数,比如根号2的平方就不是无理数。

    证明过程用到了一个函数,其中n是一个任取的大于等于1的常数。可以想像,这个函数的分子部分展开后是一个关于x的整系数多项式,最低次数为n,最高次数为2n。我们将用到这个函数的两个性质:首先,当0<x<1时,显然有0 < f(x) < 1/n!;其次,函数f及其任意阶导数在x=0和x=1处都是整数。为了证明后一个结论,首先注意到当x=0时,不管是多少阶的导数,除了常数项以外其余项都是0;常数项只可能在n<=k<=2n时出现(k表示k阶导数),但此时它等于一个整系数乘以k!/n!,显然也是个整数。另外,由于f(x)=f(1-x),根据复合函数的微分法我们立即得到对任意x都成立,当然也就有

查看更多 »

Dec 3

    《数学与猜想》里引用了一段欧拉的这篇经典的研究报告,写的非常精彩。你可以从中看到一个数学家是如何进行发现、归纳、猜想和论证的。你可以看到两个完全不同的数学模型里出现了惊人的巧合,通过挖掘它们之间的内在联系,最终完成了伟大的统一。
    没扫描仪,拿相机拍的,效果非常不好,见谅了!
    另外,拜托大家不要盗链下面的图片。

  
  
  
  
  

Nov 24

    Quake III公开源码后,有人在game/code/q_math.c里发现了这样一段代码。它的作用是将一个数开平方并取倒,经测试这段代码比(float)(1.0/sqrt(x))快4倍:
float Q_rsqrt( float number )
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;  // evil floating point bit level hacking
  i  = 0x5f3759df - ( i >> 1 ); // what the fuck?
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
  // y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

  #ifndef Q3_VM
  #ifdef __linux__
    assert( !isnan(y) ); // bk010122 - FPE?
  #endif
  #endif
  return y;
}


    code/common/cm_trace.c中也出现了这样一段解释sqrt(x)的函数,与上面的代码唯一不同的就是这个函数返回的是number*y:
/*
================
SquareRootFloat
================
*/
float SquareRootFloat(float number) {
    long i;
    float x, y;
    const float f = 1.5F;

    x = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;
    i  = 0x5f3759df - ( i >> 1 );
    y  = * ( float * ) &i;
    y  = y * ( f - ( x * y * y ) );
    y  = y * ( f - ( x * y * y ) );
    return number * y;
}


    这样的代码速度肯定飞快,我就不用多说了;但算法的原理是什么呢?其实说穿了也不是很神,程序首先猜测了一个接近1/sqrt(number)的值,然后两次使用牛顿迭代法进行迭代。根号a的倒数实际上就是方程1/x^2 - a = 0的一个正实根,它的导数是-2/x^3。运用牛顿迭代公式x' = x - f(x)/f'(x),式子化简为x' = x * (1.5 - 0.5a * x^2)。迭代几次后,x的值将趋于1/sqrt(a)。
    但这段代码真正牛B的是那个神秘的0x5f3759df,因为0x5f3759df - (i >> 1)出人意料地接近根号y的倒数。人们都不知道这个神秘的常数是怎么来的,只能把它当作神来膜拜。这个富有传奇色彩的常数到底咋回事,很少有人说得清楚。这篇论文比较科学地解释了这个常数。

Nov 24

    下面这种方法可以很有效地求出根号a的近似值:首先随便猜一个近似值x,然后不断令x等于x和a/x的平均数,迭代个六七次后x的值就已经相当精确了。
    例如,我想求根号2等于多少。假如我猜测的结果为4,虽然错的离谱,但你可以看到使用牛顿迭代法后这个值很快就趋近于根号2了:

(       4  + 2/   4     ) / 2 = 2.25
(    2.25  + 2/   2.25  ) / 2 = 1.56944..
( 1.56944..+ 2/1.56944..) / 2 = 1.42189..
( 1.42189..+ 2/1.42189..) / 2 = 1.41423..

....

      
    这种算法的原理很简单,我们仅仅是不断用(x,f(x))的切线来逼近方程x^2-a=0的根。根号a实际上就是x^2-a=0的一个正实根,这个函数的导数是2x。也就是说,函数上任一点(x,f(x))处的切线斜率是2x。那么,x-f(x)/(2x)就是一个比x更接近的近似值。代入f(x)=x^2-a得到x-(x^2-a)/(2x),也就是(x+a/x)/2。
    同样的方法可以用在其它的近似值计算中。Quake III的源码中有一段非常牛B的开方取倒函数。

Oct 19

    给出一个连续函数,某一点上的导数为正说明函数在这一点是上升的,换句话说函数从左边充分靠近该点时函数值总小于这个点,从右边靠近该点时函数值总大于这个点。但这并不等于说这一点左右是一个单增区间,也就是说该点左右任意小的邻域内函数都不是单调递增的。你能找出这样的函数来吗?

    昨天数学课上,我学到了一个比较牛B的东西:函数上某一点导数为正,该点邻域不一定形成单增区间。虽然左边的点都比该点低,右边的点都比该点高,但这并不能说明左边和右边各自都是单增的。这样的函数确实存在,而且并不是那种很怪的函数,仅仅是一个简单的初等函数:f(x) = x + 2x^2*sin(1/x)。由于x=0时函数没有定义,我们规定f(0)=0。按照导数的定义,函数在x=0时的导数值为
   Limit[ (f(0+Δx)-f(0))/(Δx-0), Δx->0 ]
= Limit[ f(Δx)/Δx, Δx->0 ]
= Limit[ 1 + 2Δx*sin(1/Δx) , Δx->0 ]
= 1

    这说明函数在x=0处的导数确实是正的。当x≠0时,按照求导法则可以求出f'(x) = 1 - 2*cos(1/x) + 4x*sin(1/x)。当|x|充分小时,最后一项可以忽略不计;此时只要1/x恰好等于2πn (n为整数),那么f'(x)保证是负的。这就告诉我们,x=0左右任意近的位置都存在导数为负的情况,这样不管邻域范围多小总能找到一个函数值在减小的地方。
    其实,看一下f(x)的函数图象,你会立即明白这是怎么回事。这个函数越接近原点抖动频率越快(到原点时“周期”无限小),同时振幅也越小(到原点时振幅为0,这样可以保证导数存在);但这个函数总的来说呈上升趋势。因此,这个函数才有我们前面提到的奇怪性质。

Dec 23

    我们年级有许多漂亮的MM。一班有7个左右吧,二班大概有4个,三班最多,16个,四班最可怜,一个漂亮的MM都没有,五班据说有1个。如果用一个函数“f(班级)=漂亮MM的个数”,那么我们可以把上述信息表示成:f(1)=7,f(2)=4,f(3)=16,f(4)=0,f(5)=1,等等。
    生成函数(也有叫做“母函数”的,但是我觉得母函数不太好听)是说,构造这么一个多项式函数g(x),使得x的n次方系数为f(n)。于是,上面的f函数的生成函数g(x)=7x+4x^2+16x^3+x^5+...。这就是传说中的生成函数了。关键是,这个有什么用呢?一会儿要慢慢说。我敢打赌这绝对会是我写过的最长的一篇文章。

    生成函数最绝妙的是,某些生成函数可以化简为一个很简单的函数。也就是说,不一定每个生成函数都是用一长串多项式来表示的。比如,这个函数f(n)=1 (n当然是属于自然数的),它的生成函数就应该是g(x)=1+x+x^2+x^3+x^4+...(每一项都是一,即使n=0时也有x^0系数为1,所以有常数项)。再仔细一看,这就是一个有无穷多项的等比数列求和嘛。如果-1<x<1,那么g(x)就等于1/(1-x)了。在研究生成函数时,我们都假设级数收敛,因为生成函数的x没有实际意义,我们可以任意取值。于是,我们就说,f(n)=1的生成函数是g(x)=1/(1-x)。

    我们举一个例子说明,一些具有实际意义的组合问题也可以用像这样简单的一个函数全部表示出来。
    考虑这个问题:从二班选n个MM出来有多少种选法。学过简单的排列与组合的同学都知道,答案就是C(4,n)。也就是说。从n=0开始,问题的答案分别是1,4,6,4,1,0,0,0,...(从4个MM中选出4个以上的人来方案数当然为0喽)。那么它的生成函数g(x)就应该是g(x)=1+4x+6x^2+4x^3+x^4。这不就是……二项式展开吗?于是,g(x)=(1+x)^4。
    你或许应该知道,(1+x)^k=C(k,0)x^0+C(k,1)x^1+...+C(k,k)x^k;但你或许不知道,即使k为负数和小数的时候,也有类似的结论:(1+x)^k=C(k,0)x^0+C(k,1)x^1+...+C(k,k)x^k+C(k,k+1)x^(k+1)+C(k,k+2)x^(k+2)+...(一直加到无穷;式子看着很别扭,自己写到草稿纸上吧,毕竟这里输入数学式子很麻烦)。其中,广义的组合数C(k,i)就等于k(k-1)(k-2)...(k-i+1)/i!,比如C(4,6)=4*3*2*1*0*(-1)/6!=0,再比如C(-1.4,2)=(-1.4)*(-2.4)/2!=1.68。后面这个就叫做牛顿二项式定理。当k为整数时,所有i>k时的C(k,i)中分子都要“越过”0这一项,因此后面C(k,k+1),C(k,k+2)之类的都为0了,与我们的经典二项式定理结论相同;不同的是,牛顿二项式定理中的指数k可以是任意实数。

    我们再举一个例子说明一些更复杂的生成函数。n=x1+x2+x3+...+xk有多少个非负整数解?这道题是学排列与组合的经典例题了。把每组解的每个数都加1,就变成n+k=x1+x2+x3+...+xk的正整数解的个数了。教材上或许会出现这么一个难听的名字叫“隔板法”:把n+k个东西排成一排,在n+k-1个空格中插入k-1个“隔板”。答案我们总是知道的,就是C(n+k-1,k-1)。它就等于C(n+k-1,n)。它关于n的生成函数是g(x)=1/(1-x)^k。这个生成函数是怎么来的呢?其实,它就是(1-x)的-k次方。把(1-x)^(-k)按照刚才的牛顿二项式展开,我们就得到了x^n的系数恰好是C(n+k-1,n),因为C(-k,n)*(-x)^n=[(-1)^n*C(n+k-1,n)]*[(-1)^n*x^n]=C(n+k-1,n)x^n。这里看晕了不要紧,后文有另一种方法可以推导出一模一样的公式。事实上,我们有一个纯组合数学的更简单的解释方法。因为我们刚才的几何级数1+x+x^2+x^3+x^4+...=1/(1-x),那么(1+x+x^2+x^3+x^4+...)^k就等于1/(1-x)^k。仔细想想k个(1+x+x^2+x^3+x^4+...)相乘是什么意思。(1+x+x^2+x^3+x^4+...)^k的展开式中,n次项的系数就是我们的答案,因为它的这个系数是由原式完全展开后k个指数加起来恰好等于n的项合并起来得到的。

    现在我们引用《组合数学》上暴经典的一个例题。很多书上都会有这类题。
    我们要从苹果、香蕉、橘子和梨中拿一些水果出来,要求苹果只能拿偶数个,香蕉的个数要是5的倍数,橘子最多拿4个,梨要么不拿,要么只能拿一个。问按这样的要求拿n个水果的方案数。
    结合刚才的k个(1+x+x^2+x^3+x^4+...)相乘,我们也可以算出这个问题的生成函数。

g(x)=(1+x^2+x^4+...)(1+x^5+x^10+..)(1+x+x^2+x^3+x^4)(1+x)
    =[1/(1-x^2)]*[1/(1-x^5)]*[(1-x^5)/(1-x)]*(1+x) (前两个分别是公比为2和5的几何级数,
                                                     第三个嘛,(1+x+x^2+x^3+x^4)*(1-x)不就是1-x^5了吗)
    =1/(1-x)^2   (约分,把一大半都约掉了)
    =(1-x)^(-2)=C(1,0)+C(2,1)x+C(3,2)x^2+C(4,3)x^3...   (参见刚才对1/(1-x)^k的展开)
    =1+2x+3x^2+4x^3+5x^4+....


    于是,拿n个水果有n+1种方法。我们利用生成函数,完全使用代数手段得到了答案!
    如果你对1/(1-x)^k的展开还不熟悉,我们这里再介绍一个更加简单和精妙的手段来解释1/(1-x)^2=1+2x+3x^2+4x^3+5x^4+....。
    1/(1-x)=1+x+x^2+x^3+x^4+...是前面说过的。我们对这个式子等号两边同时求导数。于是,1/(1-x)^2=1+2x+3x^2+4x^3+5x^4+....。一步就得到了我们所需要的东西!不断地再求导数,我们同样可以得到刚才用复杂的牛顿二项式定理得到的那个结论(自己试试吧)。生成函数还有很多其它的处理手段,比如等式两边同时乘以、除以常数(相当于等式右边每一项乘以、除以常数),等式两边同时乘以、除以一个x(相当于等式右边的系数“移一位”),以及求微分积分等。神奇的生成函数啊。
    我们用两种方法得到了这样一个公式:1/(1-x)^n=1+C(n,1)x^1+C(n+1,2)x^2+C(n+2,3)x^3+...+C(n+k-1,k)x^k+...。这个公式非常有用,是把一个生成函数还原为数列的武器。而且还是核武器。


    接下来我们要演示如何使用生成函数求出Fibonacci数列的通项公式。
    Fibonacci数列是这样一个递推数列:f(n)=f(n-1)+f(n-2)。现在我们需要求出它的生成函数g(x)。g(x)应该是一个这样的函数:
    g(x)=x+x^2+2x^3+3x^4+5x^5+8x^6+13x^7+...
    等式两边同时乘以x,我们得到:
    x*g(x)=x^2+x^3+2x^4+3x^5+5x^6+8x^7+...
    就像我们前面说过的一样,这相当于等式右边的所有系数向右移动了一位。
    现在我们把前面的式子和后面的式子相加,我们得到:
    g(x)+x*g(x)=x+2x^2+3x^3+5x^4+8x^5+...
    把这最后一个式子和第一个式子好好对比一下。如果第一个式子的系数往左边移动一位,然后把多余的“1”去掉,就变成了最后一个式子了。由于递推函数的性质,我们神奇地得到了:g(x)+x*g(x)=g(x)/x-1。也就是说,g(x)*x^2+g(x)*x-g(x)=-x。把左边的g(x)提出来,我们有:g(x)(x^2+x-1)=-x。于是,我们得到了g(x)=x/(1-x-x^2)。
    现在的任务是要把x/(1-x-x^2)还原成通项公式。这不是我们刚才的1/(1-x)^n的形式,我们要把它变成这种形式。我们发现,1-x-x^2=[1-(1-√5)x/2]*[1-(1+√5)x/2] ((1-√5)/2和(1+√5)/2是怎么算出来的?显然它们应该是x^2-x-1=0的两个根)。那么x/(1-x-x^2)一定能表示成?/[1-(1-√5)x/2]+?/[1-(1+√5)x/2]的形式(再次抱歉,输入数学公式很麻烦,将就看吧)。这是一定可以的,因为适当的?的取值可以让两个分式通分以后分子加起来恰好为一个x。?取值应该是多少呢?假设前面一个?是c1,后面那个是c2,那么通分以后分子为c1*[1-(1+√5)x/2]+c2*[1-(1-√5)x/2],它恰好等于x。我们得到这样两个式子:常数项c1+c2=0,以及一次项-c1*(1+√5)/2-c2*(1-√5)/2=1。这两个式子足够我们解出c1和c2的准确值。你就不用解了,我用的Mathematica 5.0。解出来c1=-1/√5,c2=1/√5。你不信的话你去解吧。现在,我们把x/(1-x-x^2)变成了-(1/√5)/[1-(1-√5)x/2] + (1/√5)/[1-(1+√5)x/2]。我们已经知道了1/[1-(1-√5)x/2]的背后是以(1-√5)/2为公比的等比数列,1/[1-(1+√5)x/2]所表示的数列公比为(1+√5)/2。那么,各乘以一个常数,再相加,我们就得到了Fibonacci数列的通项公式:f(n)=-(1/√5)*[(1-√5)/2]^n + (1/√5)*[(1+√5)/2]^n。或许你会问,这么复杂的式子啊,还有根号,Fibonacci数列不都是整数吗?神奇的是,这个充满根号的式子对于任何一个自然数n得到的都是整数。熟悉用特征方程解线性递推方程的同学应该知道,以上过程实质上和找特征根求解没有区别。事实上,用上面所说的方法,我们可以求出任何一个线性齐次递推方程的通项公式。什么叫做线性齐次递推呢?就是这样的递推方程:f(n)等于多少个f(n-1)加上多少个f(n-2)加上多少个f(n-3)等等。Fibonacci数列的递推关系就是线性齐次递推关系。

    我们最后看一个例子。我们介绍硬币兑换问题:我有1分、2分和5分面值的硬币。请问凑出n分钱有多少种方法。想一下刚才的水果,我们不难得到这个问题的生成函数:g(x)=(1+x+x^2+x^3+...)(1+x^2+x^4+...)(1+x^5+x^10+..)=1/[(1-x)(1-x^2)(1-x^5)]。现在,我们需要把它变成通项公式。我们的步骤同刚才的步骤完全相同。我们把(1-x)(1-x^2)(1-x^5)展开,得到1-x-x^2+x^3-x^5+x^6+x^7-x^8。我们求出-1+x+x^2-x^3+x^5-x^6-x^7+x^8=0的解,得到了以下8个解:-1,1,1,1,-(-1)^(1/5),(-1)^(2/5),-(-1)^(3/5),(-1)^(4/5)。这个不是我解出来的,我还是用的Mathematica 5.0。不是我不想解,而是我根本不会解这个8次方程。这也是为什么信息学会涉及这些东西的原因:次数稍微一高,只好交给计算机解决了。于是,(1-x)(1-x^2)(1-x^5)=(1+x)(1-x)^3(1+(-1)^(1/5) x)()()() (省略不写了)。注意那个(1-x)^3。由于等根的出现,我们不得不把(1-x)^3所包含的(1-x)和(1-x)^2因子写进一会儿的分母里,不然会导致解不出合适的c来。你可以看到很多虚数。不过没关系,这些虚数同样参与运算,就像刚才的根式一样不会影响到最后结果的有理性。然后,我们像刚才一样求出常数满足1/(1-x)(1-x^2)(1-x^5)=c1/()+c2/(1-x)+c3/(1-x)^2+c4/(1-x)^3...+c8/()。这个解太复杂了,我用Mathematica解了几分钟,打印出了起码几十KB的式子。虽然复杂,但我确实是得到了通项公式。你有兴趣的话可以尝试用Mathematica解决一下1/[(1-x)(1-x^3)] (只有1分和3分的硬币)。解c的值时可以用SolveAlways[]函数。你可以亲眼见到,一个四五行的充满虚数的式子最后总是得到正确的整数答案。

    生成函数还有很多东西,推导Catalan数列啊,指数生成函数啊,之类的。我有空再说吧,已经5000多个字了。

    huyichen一直在问那道题。很显然,那道题目和上面的兑换硬币有些联系。事实上,很多与它类似的题目都和生成函数有关。但那个题却没有什么可以利用生成函数的地方(或许我没想到吧)。或许每个max的值有什么方法用生成函数解出来,但整个题目是不大可能用生成函数解决的。
    近来有个帖子问一道“DP天牛”题目的。那个题目也是这样,很多与它类似的题目都和DP有关,但那道题却不大可能动规。我总觉得它可以归约到装箱问题(考虑体积关系,最少要几个箱子才能把物品放完),而后者貌似属于NPC。或许我错了吧,现在没事就在研究理论的东西,很久没有想过OI题了,这方面的能力已经开始退化了。

Matrix67原创
做人要厚道,转贴请注明出处

Apr 5

    近来一些人问我,听说你们文科班的不学极限啊。我回答,嗯,而且更强的是,不学极限但学导数。这个话题非常有趣。我打算详细介绍一下我们的课本是什么样子的。
    薄,封面绿色的,上书“第三册(选修I)”。
    整本书只有两章,统计和导数。
    统计里面只有三节,抽样方法、总体分布的估计(教你写正字,画表格)、总体期望值和方差的估计。最喜剧的是,在三种抽样方法中有一种(系统抽样)被认为太难了而在文科课本中被删掉。整章内容十几页搞定,然后是些题和和你们一样的随机数表。
    导数里面有六节,分别是导数的背景、概念、怎么求多项式函数的导数、极值、最值和微积分的历史。
    有人问,极限都没有,导数怎么学?这个课程的安排绝对体现了编写人员的聪明才智:他从斜率等出发,把概念用文字表达出来,把用定义求导的那个式子当作一个很长的“记作……”的符号,并在那些书上前面从来没出现过的陌生符号(lim和→)后面加注释。然后就告诉你,x^n的导数是nx^(n-1)。又有人问,三角函数呢?没有。对数和指数函数的导数都没有。“e”在这本书上压根儿没出现过。课本的意思就是说,你只要知道多项式函数怎么求导,怎么用来求极值就行了。你不用知道这是为什么。这是一个现成的工具,用就是。